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.
Best Answer
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.