I am learning mixed effect model from here.

I built a simplified model and trying to extract the $\gamma$ vector from linear mixed effect model


then calculate variance of $\gamma$. (It is a toy model and we only have random effect on intercept, so $\gamma$ is a vector not a matrix)

Why it is not equal to model summary? (model summary shows DID intercept variance is $1.148$, but variance of $\gamma$ is $1.1026$). What I am missing here?

> library(lme4)
> dat <-read.csv("http://stats.idre.ucla.edu/stat/data/hdp.csv")
> dat=dat[1:100,c("Age","Sex", "DID","mobility")]
> m <- lmer(mobility ~ 1 + Age + Sex +  (1 | DID), data = dat)
> summary(m)

Linear mixed model fit by REML ['lmerMod']
Formula: mobility ~ 1 + Age + Sex + (1 | DID)
   Data: dat

REML criterion at convergence: 228.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3152 -0.4710 -0.0252  0.9330  1.7854 

Random effects:
 Groups   Name        Variance Std.Dev.
 DID      (Intercept) 1.1480   1.0714  
 Residual             0.4532   0.6732  
Number of obs: 100, groups:  DID, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)  3.28826    0.76484   4.299
Age          0.01235    0.01178   1.048
Sexmale     -0.09126    0.14836  -0.615

Correlation of Fixed Effects:
        (Intr) Age   
Age     -0.766       
Sexmale -0.097  0.030

> var(as.numeric(getME(m, "b")))
[1] 1.102682

2 Answers 2


What you examined by var(as.numeric(getME(m, "b"))) is the variance of the conditional mode of the random effects. The conditional mode is the best linear unbiased predictors (BLUP) of the random effects and is computed as the product of the relative covariance factor of the random effects $\Lambda$ and $u$ the conditional mode of the spherical random effects variable used during the optimisation of the profiled deviance. (ie. all.equal( var( as.numeric(getME(m, "b"))), var( as.numeric( crossprod(getME(m, "Lambdat"), getME(m, "u"))) ) ) # TRUE 1.10...); see here Eq. 10.

The key to remember here is optimisation. That is because during that optimisation what we optimised against ($\theta$) is the ratio of the residual deviation (see this excellent thread on details on this matter) against the REML/ML random effects deviation; formally put the relative Cholesky factor of the random effects in $Z$. So if we multiply the estimated standard deviation of the errors $\hat{\sigma}$ with $\theta$ we will get the variance (ie. all.equal( (sigma(m) * m@theta)^2, summary(m)$varcor$DID[1] ) #TRUE 1.14.... This how the inequality is materialised. But why are these not the same?

This goes back to the issue of BLUP and REML estimation; they are not guaranteed to be the same (Elvis' answer explains that too - yes, upvote it in case you haven't already; gung's answer here is highly relevant). Note also that exactly because we work on the standard deviation plane these difference will be amplified when examining the variances.

  • $\begingroup$ Thanks for the detailed answer on the optimization process. I need to spend more time on reading the math of the reference you provided in addition to ucla online tutorial. $\endgroup$
    – Haitao Du
    Commented Apr 25, 2017 at 4:28
  • 2
    $\begingroup$ +1 but I can't say I really understood the origin of the difference. $\endgroup$
    – amoeba
    Commented Apr 25, 2017 at 7:12
  • $\begingroup$ In my understanding, the difference is due to the fact that the variance of the random effects is given by the sum of the variance of the BLUP estimates plus the variance implied by the fact that BLUP values are just estimates. See the formula here: stats.stackexchange.com/a/68126/159259 , that is also the answer to the question reported at the top, flagged as the reply to this question as well. $\endgroup$ Commented Nov 8, 2022 at 21:45

The variance in the mixed model is not the variance of the shrunken BLUP random effects but the variance of what could be called the un-shrunken fixed effects (obtained from a dummy for each group) that is treating the effect for a group as a separate estimate and not as part of distribution for which you have an estimated variance.

Consequently the mixed higher-level variance is not the estimated between group variance of the sample, but the estimated between-group variance in the population.

This might help https://www.researchgate.net/publication/225303298_Contextual_Models_of_Urban_House_Prices_A_Comparison_of_Fixed-_and_Random-Coefficient_Models_Developed_by_Expansion

and this



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