Understanding the variance of random effects in lmer() models

Question

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.

Answer 1

This is a classic one way anova. A very short answer to your question is that the variance component is made up of two terms.

$$\hat{\sigma}^2_{\alpha}=E\left[\frac{1}{48}\sum_{s=1}^{48} \alpha_s^2\right]= \frac{1}{48}\sum_{s=1}^{48}\hat{ \alpha }_s^2 +\frac{1}{48}\sum_{s=1}^{48}var(\hat{ \alpha }_s)$$

So the term you computed is the first term on the rhs (as random effects have mean zero). The second term depends on whether REML of ML is used, and the the sum of squared standard errors of your random effects.