[R-sig-ME] prediction interval for prediction from mixed effects model

[Rasmussen, Paul W - DNR]

Hi,

I have a question about the description of computing prediction intervals at the glmm.wikidot.com/faq site.  I have copied the relevant sections of code below, with the subsections that I'm concerned about underlined.  The code example for nlme uses the variance of the fixed effects (predvar) and adds to it the residual variance (fm1$sigma^2).  The code example for lme4 uses the variance of the fixed effects (pvar1) and adds to it the variance among subjects (VarCorr(fm1)$Subject[1]).  It seems to me that the nlme example is correct and that the lme4 example is incorrect - shouldn't the lme4 example add the residual variance instead of the variance among subjects?  I would think that the variance among subjects would only be appropriate when predicting for a subject not previously observed.


lme
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)

lme4

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 <- predict(fm1,newdat)
## or 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

Thank you,
Paul Rasmussen
WI DNR





	[[alternative HTML version deleted]]