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

Question

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

Answer 1

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.

Answer 2

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

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