0
$\begingroup$

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')

enter image description here

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.

enter image description here

Red line is default. Black line is what I want.

$\endgroup$
  • $\begingroup$ Possible duplicate of Line of best fit does not look like a good fit. Why? $\endgroup$ – Stephan Kolassa Jul 6 '18 at 20:28
  • $\begingroup$ Since both lines "cut through the kernel density" [contours], please explain how we are to determine, identify, or otherwise distinguish the black line from all others. What do you intend it to show? Perhaps the major axis of the second-order approximation? Beware, though, that this will not be the least squares regression line and it might have poor properties, depending on how you intend to use and interpret it. Perhaps you could explain further? $\endgroup$ – whuber Jul 6 '18 at 20:40
  • $\begingroup$ One other thing: unless you make sure to scale the axes properly, most angles in the plot will be deceptive. $\endgroup$ – whuber Jul 6 '18 at 20:49
  • $\begingroup$ Ideally, I'd like to somehow 'draw' lines through each of the kernel density circles, calculate their intercept and slopes, take an average of all those lines, and draw that averaged line. I would hope that this line would 'cut the kernel density in half', in a way. Based on my playing with the data for months, I think this would appropiately show 'higher' vs. 'lower' than expected samples (not statistically, but visually) $\endgroup$ – user3579613 Jul 6 '18 at 21:25
  • $\begingroup$ You might be on the wrong track, then: your black/gray line clearly is the wrong one in the tails, because at the right it is higher than just about all the data and at the left it eventually is lower than all the data. That's a pretty serious bias for something that's intended to split data roughly evenly into "higher" and "lower." Perhaps some of the plots I posted near the end of stats.stackexchange.com/a/71303/919 will help make this point more clearly. If by "follow the density" you mean to trace its "ridge line" (conditional modes), the red line already does that. $\endgroup$ – whuber Jul 6 '18 at 21:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.