Interval around variance components

Question

I am interested in (confidence) intervals or standard errors or some thing similar (sampling based?) around the estimates of variance components for the random effects in lme4::lmer models.

I am sure that I am oversimplyfying things (and I am aware that standard errors are deemed not to be the correct thing here), but I would greatly appreciate some lines of code that give some measure of the confidence in these estimates.

In the example from the help page of VarCorr I'd like some confidence aroun the three values in the Variance colum of this output:

data(Orthodont, package="nlme")
fm1 <- lmer(distance ~ age + (age|Subject), data = Orthodont)
vc <- VarCorr(fm1)
print(vc,comp=c("Variance"))

Actually, confint(fm1,which="theta_")^2 will work for the CIs on the variances, but won't make sense for the CIs on the covariance. — Ben Bolker
@BenBolker Thanks a lot. Apparently this is quite easy. Could you still help me with the interpretation? The print(vc,comp=c("Variance")) gives three values for the variance: 'Subject (Intercept)', 'Subject age', and 'Residual'. The confint(fm1,which="theta_")^2 gives four intervals '.sig01', '.sig02', '.sig03', and '.sigma'. How do they relate? Especially, is there any CI for the residual variance? — Andreas

Ben Bolker Ben Bolker · Accepted Answer · 2014-10-28T21:52:39

Preliminaries: (I'm going to take the liberty of using a different example, as there seems to be something odd about Orthodont which I will have to look into ...)

library("lme4")
fm1 <- lmer(Reaction ~ Days + (Days|Subject), data = sleepstudy)
vc <- VarCorr(fm1)
print(vc,comp=c("Variance"))

Compute confidence intervals. Using oldNames=FALSE gives us more meaningful names to work with:

cc <- confint(fm1,which="theta_",oldNames=FALSE)

print(cc,digits=3)
##                               2.5 % 97.5 %
## sd_(Intercept)|Subject       14.382 37.716
## cor_Days.(Intercept)|Subject -0.482  0.685
## sd_Days|Subject               3.801  8.753
## sigma                        22.898 28.858

Squaring cc gives us confidence intervals for the variances -- unfortunately the confidence interval for the squared correlation coefficient is rather useless (we would probably prefer the confidence interval for the covariance, which would take more work). Drop the second (correlation) row, and drop the soon-to-be-misleading sd_/cor_ tags at the beginnings of the row names:

ccsq <- cc[-2,]^2
rownames(ccsq) <- gsub("^.*_","",rownames(ccsq))
ccsq

##                         2.5 %     97.5 %
## (Intercept)|Subject 206.82812 1422.49637
## Days|Subject         14.44885   76.62172
## sigma               524.33063  832.78396

(Note the last line is the 95% CI for the residual variance, not the residual standard deviation, even though I haven't changed the name ...)

Interval around variance components

1 Answers