3
$\begingroup$

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)
head(data0)

# -- 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))
[1] 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?

$\endgroup$
  • 1
    $\begingroup$ Does this answer you question? stats.stackexchange.com/questions/69882/… $\endgroup$ – Steve Sep 16 '14 at 13:42
  • $\begingroup$ @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. $\endgroup$ – giordano Sep 17 '14 at 10:59
  • 1
    $\begingroup$ 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. $\endgroup$ – Steve Sep 17 '14 at 11:24
  • 1
    $\begingroup$ @Steve Oviously, the formula is valid for REML but not EML. $\endgroup$ – giordano May 1 '16 at 17:43
2
$\begingroup$

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
$\endgroup$
  • 1
    $\begingroup$ 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 $\endgroup$ – giordano May 2 '16 at 14:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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