summarize extent of pooling or shrinkage in multilevel models estimated with lmer()

Question

I am using lmer() in the "lme4" package to estimate multilevel models. The models include random intercepts for each group in my data. To fix ideas, here is a toy example:

library(lme4)
data(iris)
foo <- lmer(Sepal.Length ~ Sepal.Width + (1 | Species), 
            data = iris)

I would like to summarize the extent to which models like this shrink the estimates of the intercepts toward the grand mean of all intercepts, relative to the estimates that I would get from a simpler model, estimated with lm(), that includes dummy variables for each group. How may I do this?

In their book, Gelman and Hill (2007, 477-80) refer to this summary statistic as a "pooling factor" and they note that others sometimes speak of a related "shrinkage factor." In their notation, the intercepts are $\theta_k = \hat{\theta}_k + \epsilon_k$ for $k = 1, \ldots, K$. They suggest estimating a summary of the extent to which the variance of the residuals $\epsilon_k$ is reduced by the pooling of the multilevel model: $$ \DeclareMathOperator*{\V}{V} L = 1 - \frac{\V_\limits{k=1}^KE(\epsilon_k)}{E\left(\V_\limits{k=1}^K \epsilon_k\right) }. $$ They give instructions for computing this quantity in BUGS. But is there a relatively simple way to do it in R?

Perhaps the numerator in the equation above corresponds to sigma(foo)^2, but I'm not sure of that. And I don't have good ideas about how to compute the denominator. Can this information be extracted from objects created by lmer()?

I've looked through CrossValidated and haven't found any posts on this point.

Answer 1

Fleshing out Dimitris' comment, you can look at this by considering the estimates you get from lmer and lm. Using your lmer model, we can ask for the estimated intercepts and slopes with the coef() function. The intercept listed in coef() is based on the overall (fixed effect/grand mean) intercept plus/minus each group's random effect deviation off the fixed intercept:

coef(foo)
$Species
           (Intercept) Sepal.Width
setosa        2.277601   0.7971543
versicolor    3.726677   0.7971543
virginica     4.214224   0.7971543

Now you can compare these to the intercepts (means) you get from an OLS model with dummy variables for each of the Species:

summary(foo_fe)

Call:
lm(formula = Sepal.Length ~ -1 + Sepal.Width + as.factor(Species), 
    data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.30711 -0.25713 -0.05325  0.19542  1.41253 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
Sepal.Width                    0.8036     0.1063   7.557 4.19e-12 ***
as.factor(Species)setosa       2.2514     0.3698   6.089 9.57e-09 ***
as.factor(Species)versicolor   3.7101     0.3010  12.326  < 2e-16 ***
as.factor(Species)virginica    4.1982     0.3223  13.027  < 2e-16 ***

So the shrinkage going on with this data is not extensive likely because each of the groups has the same (rather large) number of observations and also because the level 2 (Species) intercept variance is much higher than the level 1 (within Species) variance:

 Groups   Name        Variance Std.Dev.
 Species  (Intercept) 1.0198   1.010   
 Residual             0.1918   0.438   
Number of obs: 150, groups:  Species, 3