Hello,

I am trying to adapt some code for confidence intervals on lme predictions that I found on the mixed models wiki FAQ page (http://glmm.wikidot.com/faq):


library(nlme)
fm1 <- lme(distance ~ age*Sex, random = ~ 1 + age | Subject, data = Orthodont)

plot(Orthodont)
newdat <- expand.grid(age=c(8,10,12,14), Sex=c("Male","Female"))
newdat$pred <- predict(fm1, newdat, level = 0)

Designmat <- model.matrix(eval(eval(fm1$call$fixed)[-2]), newdat[-3])
predvar <- diag(Designmat %*% fm1$varFix %*% t(Designmat))
newdat$SE <- sqrt(predvar)
newdat$SE2 <- sqrt(predvar+fm1$sigma^2)

library(ggplot2)
pd <- position_dodge(width=0.4)
g0 <- ggplot(newdat,aes(x=age,y=pred,colour=Sex))+
   geom_point(position=pd)
g0 + geom_linerange(aes(ymin=pred-2*SE,ymax=pred+2*SE), position=pd)

## prediction intervals
g0 + geom_linerange(aes(ymin=pred-2*SE2,ymax=pred+2*SE2), position=pd)


However, I am confused about sigma, which when I compared with VarCorr (fm1), appears to be the standard deviation of the residual, which is the within-group variance (I think):

newdat$SE2 <- sqrt(predvar+fm1$sigma^2)

> fm1$sigma
[1] 1.31004

> VarCorr(fm1)
Subject = pdLogChol(1 + age)
                    Variance          StdDev         Corr
(Intercept)     5.78643478     2.4055009     (Intr)
age               0.03252449     0.1803455     -0.668
Residual        1.71620366     1.3100396


My questions are: why use the residual rather than the intercept?  The lme4 example (code below) uses the Intercept variance instead, why are they different?


library(lme4)
library(ggplot2) # Plotting
data("Orthodont",package="MEMSS")
fm1 <- lmer(
    formula = distance ~ age*Sex + (age|Subject)
    , data = Orthodont
)
newdat <- expand.grid(
    age=c(8,10,12,14)
    , Sex=c("Male","Female")
    , distance = 0
)
mm <- model.matrix(terms(fm1),newdat)
newdat$distance <- mm %*% fixef(fm1)
pvar1 <- diag(mm %*% tcrossprod(vcov(fm1),mm))
tvar1 <- pvar1+VarCorr(fm1)$Subject[1]  ## must be adapted for more complex models
newdat <- data.frame(
    newdat
    , plo = newdat$distance-2*sqrt(pvar1)
    , phi = newdat$distance+2*sqrt(pvar1)
    , tlo = newdat$distance-2*sqrt(tvar1)
    , thi = newdat$distance+2*sqrt(tvar1)
)
#plot confidence
g0 <- ggplot(newdat, aes(x=age, y=distance, colour=Sex))+geom_point()
g0 + geom_errorbar(aes(ymin = plo, ymax = phi))+
    opts(title="CI based on fixed-effects uncertainty ONLY")
#plot prediction
g0 + geom_errorbar(aes(ymin = tlo, ymax = thi))+
    opts(title="CI based on FE uncertainty + RE variance")

This is the key part:

tvar1 <- pvar1+VarCorr(fm1)$Subject[1]

> VarCorr(fm1)$Subject[1]
[1] 5.786427

$Subject
                (Intercept)         age
(Intercept)   5.7864269        -0.28962636
age            -0.2896264         0.03252433

attr(,"stddev")
                    (Intercept)         age
                    2.405499        0.180345

attr(,"correlation")
                        (Intercept)       age
(Intercept)          1.000000       -0.667619
age                   -0.667619      1.000000

attr(,"sc")
[1] 1.31004

Ultimately, I would like to adapt the lme procedure for my model, which looks like this:

model<-lme(y~B*C*D, random=~1|A/C/D)

Which variance term should I use to find confidence intervals based on fixed effect uncertainty + random error variance?

Many thanks for any advice

Regards,
Elizabeth

This email may be confidential and subject to legal privilege, it may
not reflect the views of the University of Canterbury, and it is not
guaranteed to be virus free. If you are not an intended recipient,
please notify the sender immediately and erase all copies of the message
and any attachments.

Please refer to http://www.canterbury.ac.nz/emaildisclaimer for more
information.

	[[alternative HTML version deleted]]

