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


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.

  • $\begingroup$ 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 ? $\endgroup$ – Stéphane Laurent Aug 22 '13 at 21:50
  • $\begingroup$ Here is a very similar question, with a somehow different answer $\endgroup$ – Arne Jonas Warnke Jan 8 '16 at 15:19

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.

  • 2
    $\begingroup$ 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 $\endgroup$ – nomad545 Aug 22 '13 at 23:57
  • 2
    $\begingroup$ 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. $\endgroup$ – ndoogan May 16 '14 at 16:19
  • 1
    $\begingroup$ @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? $\endgroup$ – user1357015 Apr 7 '15 at 22:45
  • 1
    $\begingroup$ @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) $ $\endgroup$ – probabilityislogic Apr 9 '15 at 15:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.