theta.hat.sp <- matrix(rnorm(1000 * J, mean = Schools$y, sd = Schools$sigma), nrow = 1000, byrow = TRUE) df <- data.frame(y = Schools$y, median = apply(theta.hat.sp, 2, median), lower = apply(theta.hat.sp, 2, function(x) quantile(x, .025)), upper = apply(theta.hat.sp, 2, function(x) quantile(x, .975))) separating <- ggplot(df, aes(x = y, y = median)) + geom_linerange(aes(ymin = lower, ymax = upper