# Mixed model with lmer: Variance of residuals should give the same as level 1 variance?

I expected that the variance of residuals from a mixed model computed by, for example, lmer should give the same as the residual variance from the summary output.

set.seed(123)
# -- Data simulation
group  <- rep(letters[1:10], each=10)
groupx <- rep(rnorm(10,0,5), each=10)
x0     <- 1:10
y0     <- rep(3*x0,10)
y      <- y0+groupx + rnorm(100,0,3)
data0  <- data.frame(grp=group,x=x0,y=y)

# -- library(lattice)
xyplot(y ~ x|grp, data=data0)

# -- Model
f0  <- lmer(y ~ x + (1|grp), data= data0, REML=FALSE)
f0s <- summary(f0)
f0s

# -- Compare level 1 variance (Residual):
as.data.frame(VarCorr(f0))[,c("grp","vcov")]
# -- with variance of residuals:
var(resid(f0))


In my example I got:

> as.data.frame(VarCorr(f0))[,c("grp","vcov")]
grp     vcov
1      grp 20.621959
2 Residual 7.090458

> var(resid(f0))
 6.469677


I expected that var(resid(f0)) should be exactly the same as the variance of Residuals from summary. Could someone gives me a hint what is wrong in my argumentation?

UPDATE 1: According to the formula in link

$\hat{\sigma}^2_{\alpha}=E\left[\frac{1}{10}\sum_{s=1}^{10} \alpha_s^2\right]= \frac{1}{10}\sum_{s=1}^{10}\hat{ \alpha }_s^2 +\frac{1}{10}\sum_{s=1}^{48}var(\hat{ \alpha }_s)$

I tried to compute the variance of the random variable:

> 1/10 * sum((se.ranef(f0)$grp)^2) + var(ranef(f0)$grp)
(Intercept)
(Intercept)    22.83712


which seems quite different compared to 20.621959 in the summary above.

With method REML I got the following results:

f1  <- lmer(y ~  x + (1|grp), data= data0, REML=TRUE)
f1s <- summary(f1)
as.data.frame(VarCorr(f1))[,c("grp","vcov")]
1/10 * sum((se.ranef(f1)$grp)^2) + var(ranef(f1)$grp)

> f1  <- lmer(y ~  x + (1|grp), data= data0, REML=TRUE)
>     f1s <- summary(f1)
>     as.data.frame(VarCorr(f1))[,c("grp","vcov")]
grp      vcov
1      grp 22.984104
2 Residual  7.170126
>     1/10 * sum((se.ranef(f1)$grp)^2) + var(ranef(f1)$grp)
(Intercept)
(Intercept)     22.984


Obviously, the formula is correct when method REML is used.

With

var(resid(f0))*99/90


I got a SD for residuals of 7.116645 which is near to 7.09. But how to justify a degree of freedom of 90? Or does it makes no sense to speak about df in mixed models?

• Does this answer you question? stats.stackexchange.com/questions/69882/… Sep 16, 2014 at 13:42
• @Steve Thanks a lot for the link. Supposing that the statement at the link is also valid for EML (and not only for REML) it does answer my question. Sep 17, 2014 at 10:59
• I think so.. I played around a bit with different simulations and formulas and I think that's the issue. If you have a very big sample you see that the estimator and the var of residuals are pretty much the same, which makes some sense. Sep 17, 2014 at 11:24
• @Steve Oviously, the formula is valid for REML but not EML. May 1, 2016 at 17:43

You might want to have a look at Interpretation of various output of "lmer" function in R.

Basically the estimated residual variance $\hat{\sigma}$ is an estimation for a population parameter $\sigma$. Whereas the estimated variance of your residual is an estimation based on a subset of your population. You will encounter the same thing in simple linear models:

set.seed(123)
x <- rnorm(50)
y <- 1+2*x + rnorm(50, sd = 2)
mod <- lm(y~x)
summary(mod)\$sigma #1.828483
sd(resid(mod))     #1.809729

• Thank you very much for this hint. I can compute the sgima from the residuals: sqrt(sum(resid(mod)2/(50−2))) or sd(resid(mod))∗sqrt(49/48) gives both 1.828483. But I cannot reproduce the standard deviation of the random effect and the residuals in lmer (see update 1). – giordano 20 hours ago May 2, 2016 at 14:08