# 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/… – Steve Sep 16 '14 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. – giordano Sep 17 '14 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. – Steve Sep 17 '14 at 11:24
• @Steve Oviously, the formula is valid for REML but not EML. – giordano May 1 '16 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 – giordano May 2 '16 at 14:08