have a large dataset of gene expression from ~10,000 patient samples (TCGA), and I'm plotting a predicted expression value (x) and the actual observed value (y) of a certain gene signature. For my downstream analysis, I need to draw a precise line through the plot and calculate different parameters in samples above/below the line. No matter how I draw a line through the data (geom_smooth(method = 'lm', 'glm', 'gam', or 'loess')), the line always seems imperfect - it doesn't cut through the data to my liking (red line is lm in figure). After playing around for a while, I realized that the 2d kernel density lines (geom_density2d) actually do a good job of showing the slope/trends of my data, so I manually drew a line that kind of cuts through the density lines (black line in figure).
My question: how can I automatically draw a line that cuts through the kernel density lines, as for the black line in the figure? (Rather than manually playing with different intercepts and slopes till something looks good).
The best approach I can think of is to somehow calculate intercept and slope of the longest diameter for each of the kernel lines, take an average of all those intercepts and slopes and plot that line, but that's a bit out of my league. Maybe someone here has experience with this and can help?
A more hacky approach may be getting the x,y coords of each kernel density line from ggplot_build, and going from there, but it feels too hacky (and is also out of my league).
Thanks!
Reprex:
library(MASS)
set.seed(123)
samples = 10000
r = 0.9
data = mvrnorm(n=samples, mu=c(0, 0), Sigma=matrix(c(2, r, r, 2), nrow=2))
x = data[, 1] # standard normal (mu=0, sd=1)
y = data[, 2] # standard normal (mu=0, sd=1)
#plot(x, y)
test.df = data.frame(x = x, y = y)
lm(y ~ x, test.df)
ggplot(test.df, aes(x, y)) +
geom_point(color = 'grey') +
#geom_hex(bins = 90) +
#stat_density2d(aes(fill = ..density..), geom = "raster", contour = FALSE) +
geom_density2d(color = 'red', lwd = 0.5, contour = T, h = c(2,2)) +
geom_smooth(method = "glm", se = F, lwd = 1, color = 'red') +
geom_abline(intercept = 0, slope = 0.7, lwd = 1, col = 'black')
EDIT: I'm predicting expression of Gene A based on a bunch of genes (Genes J, K, L) (this is x-axis), and then plotting that against the actual observed expression of Gene A (y-axis). I then hypothesize that expression of Gene B causes 'lower than expected' expression of Gene A. I then have a bunch of downstream steps to compare samples with higher vs. lower than expected Gene A expression. The figure below seems to suggest my hypothesis is correct.
Thus, I don't actually need the black line to be 'more statistically correct' than the off the shelf red 'lm' line, I just want to separate data into 'higher' or 'lower' than expected. As you can see in the figure below, the red line has issues towards the tails of the line. And further, the black line follows the 'kernel density' of the population. So if I can get the black line automatically, it wouldn't be totally manual - I'd just choose to follow the density of the samples.
Red line is default. Black line is what I want.
