# 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. • Do you have some knowledge about the way parameters are estimated with lmer()? It seems that you postulate that$\sigma^2_\alpha$is estimated by the empirical variance of the estimated random effects$\hat\alpha_s$. The description of your model is not clear (perharps$y_i$should be$y_{is}$). Is it a balanced design ? Aug 22, 2013 at 21:50 • Here is a very similar question, with a somehow different answer Jan 8, 2016 at 15:19 ## 1 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. • OK, got it! So, the sum of the squared SEs of the REs -- 1/48 * sum((se.ranef(mlm1)$State)^2) -- is 0.004557198. The variance of the point estimates of the REs (obtained, as above, using var(ranef(mlm1)$State)) is 0.001800869. The sum is 0.006358067, which is the variance reported using summary() on the lmer() model above, to 4 or 5 digits at least. Many thanks @probability Aug 22, 2013 at 23:57 • For those looking to this answer and the comment for help, note that nomad545 has also made use of the arm R package for the se.ranef() function. May 16, 2014 at 16:19 • @probabilityislogic: Can you provide some more detail how that equation was calculated? Specifically, how was the second equality achieved? Also, shoudn't there be a hat on the alpha after the first equality? Apr 7, 2015 at 22:45 • @user1357015 - one way to see this is to look at the gradient of the (marginal) log likelihood after integrating out the random effects. That is, differentiate the likelihood$ Y\sim Normal (1_n\alpha_0,\Sigma) $where$\Sigma=I_n\sigma^2_e+\sigma^2_{\alpha} ZZ^T$is the "unconditional" variance of Y. If you do this (plus using some manipulations) you get the above equality. The second equality follows because$ E (\alpha_s)=0$(under the model) meaning$ var (\alpha_s)=E (\alpha_s^2) \$ Apr 9, 2015 at 15:32
• I may be mistaken, but in the post by @probabilityislogic, it seems that "first term on the rhs" should be "second term on the rhs," and that "second term" should be "first term." The comment from nomad545 lends some support to this interpretation. Jun 7, 2020 at 13:43