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.

## 2d shrikage for random slope, random intercept model
library ( lme4 )
library ( ggplot2 )
library ( grid )
library ( mvtnorm )
library ( reshape2 )
library ( ellipse )
# fit models
m 0 <- lm ( Reaction ~ Days * Subject , data = sleepstudy )
m 1 <- lmer ( Reaction ~ Days + ( Days | Subject ), sleepstudy )
# extract fixed effect estimates
intercepts <- grepl ( "Days" , names ( coef ( m 0 ))) == FALSE
m 0 _ intercepts <- coef ( m 0 )[ intercepts ]
m 0 _ intercepts [ -1 ] <- m 0 _ intercepts [ -1 ] + m 0 _ intercepts [ 1 ]
slopes <- ! intercepts
m 0 _ slopes <- coef ( m 0 )[ slopes ]
m 0 _ slopes [ -1 ] <- m 0 _ slopes [ -1 ] + m 0 _ slopes [ 1 ]
# extract random effect estimates
m 1 _ intercepts <- coef ( m 1 ) $ Subject [, 1 ]
m 1 _ slopes <- coef ( m 1 ) $ Subject [, 2 ]
d <- data.frame ( interceptF = m 0 _ intercepts ,
slopeF = m 0 _ slopes ,
interceptR = m 1 _ intercepts ,
slopeR = m 1 _ slopes
)
# 95% bivariate normal ellipse for random effects
df_ell <- data.frame ( ellipse ( VarCorr ( m 1 ) $ Subject , centre = fixef ( m 1 )))
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 ( m 1 ),
sigma = VarCorr ( m 1 ) $ 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