# Why is the estimated variance of a random effect not equal to the variance of its BLUP estimates? [duplicate]

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

$$y=X\beta+Z\gamma+\epsilon$$

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=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

• I got to give to you - this is a messy little question. (+1) – usεr11852 Apr 24 '17 at 22:37
• One way to think about this is the difference between the among-group variance (=Vg) and the observed difference among groups (=Vg+Ve/n); I'd expect the variance of the BLUPs to be intermediate between these two values. – Ben Bolker Apr 25 '17 at 12:09
• – gung - Reinstate Monica Apr 25 '17 at 20:52
• @gung: That was my edit about! :D – usεr11852 Apr 25 '17 at 20:54
• @usεr11852, sorry I missed it. I think you edited while I was writing the comment. – gung - Reinstate Monica Apr 25 '17 at 21:09

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.

• 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. – Haitao Du Apr 25 '17 at 4:28
• +1 but I can't say I really understood the origin of the difference. – amoeba Apr 25 '17 at 7:12

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.

and this

https://www.researchgate.net/publication/252146040_Do_multilevel_models_ever_give_different_results