Make ggplot with regression line and normal distribution overlay

user_11409494 0
Viewed 580

I am trying to make a plot to show the intuition behind logistic (or probit) regression. How would I make a plot that looks something like this in ggplot?

Wolf & Best, The Sage Handbook of Regression Analysis and Causal Inference, 2015, p. 155

(Wolf & Best, The Sage Handbook of Regression Analysis and Causal Inference, 2015, p. 155)

Actually, what I would rather even do is have one single normal distribution displayed along the y axis with mean = 0, and a specific variance, so that I can draw horizontal lines going from the linear predictor to the y axis and sideways normal distribution. Something like this:

What this is supposed to show (assuming I haven't misunderstood something) is equation. I haven't had much success so far...

library(ggplot2)

x <- seq(1, 11, 1)
y <- x*0.5

x <- x - mean(x)
y <- y - mean(y)

df <- data.frame(x, y)

# Probability density function of a normal logistic distribution 
pdfDeltaFun <- function(x) {
  prob = (exp(x)/(1 + exp(x))^2)
  return(prob)
}

# Tried switching the x and y to be able to turn the 
# distribution overlay 90 degrees with coord_flip()
ggplot(df, aes(x = y, y = x)) + 
  geom_point() + 
  geom_line() + 
  stat_function(fun = pdfDeltaFun)+ 
  coord_flip()

enter image description here

logistic-regression ggplot2 r
1 Answers

I think this comes pretty close to the first illustration you give. If this is a thing you don't need to repeat many times, it is probably best to compute the density curves prior to plotting and use a seperate dataframe to plot these.

library(ggplot2)

x <- seq(1, 11, 1)
y <- x*0.5

x <- x - mean(x)
y <- y - mean(y)

df <- data.frame(x, y)

# For every row in `df`, compute a rotated normal density centered at `y` and shifted by `x`
curves <- lapply(seq_len(NROW(df)), function(i) {
  mu <- df$y[i]
  range <- mu + c(-3, 3)
  seq <- seq(range[1], range[2], length.out = 100)
  data.frame(
    x = -1 * dnorm(seq, mean = mu) + df$x[i],
    y = seq,
    grp = i
  )
})
# Combine above densities in one data.frame
curves <- do.call(rbind, curves)


ggplot(df, aes(x, y)) +
  geom_point() +
  geom_line() +
  # The path draws the curve
  geom_path(data = curves, aes(group = grp)) +
  # The polygon does the shading. We can use `oob_squish()` to set a range.
  geom_polygon(data = curves, aes(y = scales::oob_squish(y, c(0, Inf)),group = grp))

The second illustration is pretty close to your code. I simplified your density function by the standard normal density function and added some extra paramters to stat function:

library(ggplot2)

x <- seq(1, 11, 1)
y <- x*0.5

x <- x - mean(x)
y <- y - mean(y)

df <- data.frame(x, y)

ggplot(df, aes(x, y)) +
  geom_point() +
  geom_line() +
  stat_function(fun = dnorm,
                aes(x = after_stat(-y * 4 - 5), y = after_stat(x)),
                xlim = range(df$y)) +
  # We fill with a polygon, squishing the y-range
  stat_function(fun = dnorm, geom = "polygon",
                aes(x = after_stat(-y * 4 - 5), 
                    y = after_stat(scales::oob_squish(x, c(-Inf, -1)))),
                xlim = range(df$y))

user_113748270
Related