# Understanding the variance of random effects in lmer() models

I’m having trouble understanding the output of my lmer() model. It is a simple model of an outcome variable (Support) with varying State intercepts / State random effects:

mlm1 <- lmer(Support ~ (1 | State))


The results of summary(mlm1) are:

Linear mixed model fit by REML
Formula: Support ~ (1 | State)
AIC   BIC logLik deviance REMLdev
12088 12107  -6041    12076   12082
Random effects:
Groups   Name        Variance  Std.Dev.
State    (Intercept) 0.0063695 0.079809
Residual             1.1114756 1.054265
Number of obs: 4097, groups: State, 48

Fixed effects:
Estimate Std. Error t value
(Intercept)  0.13218    0.02159   6.123


I take it that the variance of the varying state intercepts / random effects is 0.0063695. But when I extract the vector of these state random effects and calculate the variance

var(ranef(mlm1)\$State)


The result is: 0.001800869, considerably smaller than the variance reported by summary().

As far as I understand it, the model I have specified can be written:

$y_i = \alpha_0 + \alpha_s + \epsilon_i, \text{ for } i = \{1, 2, ..., 4097\}$

$\alpha_s \sim N(0, \sigma^2_\alpha), \text{ for } s = \{1, 2, ..., 48\}$

If this is correct, then the variance of the random effects ($\alpha_s$) should be $\sigma^2_\alpha$. Yet these are not actually equivalent in my lmer() fit.