Visualizing bivariate shrinkage
[This article was first published on Ecology in silico, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Inspired by this post about visualizing shrinkage on Coppelia, and this thread about visualizing mixed models on Stack Exchange, I started thinking about how to visualize shrinkage in more than one dimension. One might find themselves in this situation with a varying slope, varying intercept hierarichical (mixed effects) model, a model with two varying intercepts, etc. Then I found a great bivariate shrinkage plot in Doug Bates’ lme4 book, Figure 3.8.
Here is a snippet to reproduce a similar bivariate shrinkage plot in ggplot2, adding a color coded probability density surface and contours for the estimated multivariate normal distribution of random effects, using the same sleep study data that Bates used.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 | ## 2d shrikage for random slope, random intercept model
library(lme4)
library(ggplot2)
library(grid)
library(mvtnorm)
library(reshape2)
library(ellipse)
# fit models
m0 <- lm(Reaction ~ Days * Subject, data=sleepstudy)
m1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
# extract fixed effect estimates
intercepts <- grepl("Days", names(coef(m0))) == FALSE
m0_intercepts <- coef(m0)[intercepts]
m0_intercepts[-1] <- m0_intercepts[-1] + m0_intercepts[1]
slopes <- !intercepts
m0_slopes <- coef(m0)[slopes]
m0_slopes[-1] <- m0_slopes[-1] + m0_slopes[1]
# extract random effect estimates
m1_intercepts <- coef(m1)$Subject[, 1]
m1_slopes <- coef(m1)$Subject[, 2]
d <- data.frame(interceptF = m0_intercepts,
slopeF = m0_slopes,
interceptR = m1_intercepts,
slopeR = m1_slopes
)
# 95% bivariate normal ellipse for random effects
df_ell <- data.frame(ellipse(VarCorr(m1)$Subject, centre=fixef(m1)))
names(df_ell) <- c("intercept", "slope")
# bivariate normal density surface
lo <- 200 # length.out of x and y grid
xvals <- seq(from = min(df_ell$intercept) - .1 * abs(min(df_ell$intercept)),
to = max(df_ell$intercept) + .05 * abs(max(df_ell$intercept)),
length.out = lo)
yvals <- seq(from = min(df_ell$slope) - .4 * abs(min(df_ell$slope)),
to = max(df_ell$slope) + .1 * abs(max(df_ell$slope)),
length.out = lo)
z <- matrix(0, lo, lo)
for (i in 1:lo) {
x = xvals
y = yvals[i]
z[,i] = dmvnorm(cbind(x,y),
mean = fixef(m1),
sigma = VarCorr(m1)$Subject)
}
mv_ranef <- melt(z)
names(mv_ranef) <- c("x", "y", "z")
mv_ranef$x <- xvals[mv_ranef$x]
mv_ranef$y <- yvals[mv_ranef$y]
p <- ggplot(d) +
geom_raster(aes(x=x, y=y, fill=z), data=mv_ranef) +
scale_fill_gradient(low="white", high="red") +
guides(fill=FALSE) +
geom_path(data=df_ell, aes(x=intercept, y=slope), size=.5) +
geom_contour(aes(x=x, y=y, z=z), data=mv_ranef, size=.1, color="black") +
geom_segment(aes(x=interceptF, y=slopeF,
xend=interceptR, yend=slopeR),
arrow = arrow(length = unit(0.3, "cm")),
alpha=.7) +
xlab("Estimated intercepts") +
ylab("Estimated slopes") +
theme(legend.position="none") +
ggtitle("Bivariate shrinkage plot") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p
|
To leave a comment for the author, please follow the link and comment on their blog: Ecology in silico.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.