ggplot with posterior distribution plotted over geom_smooth

Question

I'd like to recreate this plot from this blog post on Posterior predicted distribution for linear regression in JAGS) using ggplot?

Knowing all the extras available for ggplot, what methods are there to go about this?

Answer 1

Here's a bare-bones example of what I've come up with using density to add a geom_path.

library(ggplot2)
#mydat <- read.csv("HtWt30.csv")

mydat <- structure(list(male = c(0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 
                                 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 
                                 0L, 0L, 1L, 0L), 
                        height = c(64, 62.3, 67.9, 64.2, 64.8, 57.5, 
                                   65.6, 70.2, 63.9, 71.1, 66.5, 68.1, 62.9, 75.1, 64.6, 69.2, 68.1, 
                                   72.6, 63.2, 64.1, 64.1, 71.5, 76, 69.7, 73.3, 61.7, 66.4, 65.7, 
                                   68.3, 66.9), 
                        weight = c(136.4, 215.1, 173.6, 117.3, 123.3, 96.5, 
                                   178.3, 191.1, 158, 193.9, 127.1, 147.9, 119, 204.4, 143.4, 124.4, 
                                   140.9, 164.7, 139.8, 110.2, 134.1, 193.6, 180, 155, 188.2, 187.4, 
                                   139.2, 147.9, 178.6, 111.1)), 
                   .Names = c("male", "height", "weight"), class = "data.frame", row.names = c(NA, -30L))


# smooth plot
g_smooth <- ggplot(mydat, aes(x = height, y = weight)) + geom_smooth()

# fake posterior at a height = 60
p60 <- data.frame(x = 60, y = rnorm(1000, mean = 145, sd = 10))
# density kernel 
d60 <- density(p60$y)
# calculate scaling factor so that density covers 1/20 of full x range
density_scaling <- ((max(mydat$height) - min(mydat$height)) / 20) / max(d60$y)
# convert to points
d60points <- data.frame(y = d60$x, x = 60 + d60$y * density_scaling)
# add path to plot
g_smooth <- g_smooth + geom_path(data = d60points, aes(x = x, y = y))

# fake posterior at a height = 70
p70 <-  data.frame(x = 60, y = rnorm(1000, mean = 165, sd = 10))
# density kernel 
d70 <- density(p70$y)
# calculate scaling factor so that density covers 1/20 of full x range
density_scaling <- ((max(mydat$height) - min(mydat$height)) / 20) / max(d70$y)
# convert to points
d70points <- data.frame(y = d70$x, x = 70 + d70$y * density_scaling)
# add path to plot
g_smooth <- g_smooth + geom_path(data = d70points, aes(x = x, y = y))

g_smooth