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

Rasmussen, Paul W - DNR PaulW.Rasmussen at wisconsin.gov
Wed Jan 14 21:22:25 CET 2015


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.

fm1 <- lme(distance ~ age*Sex, random = ~ 1 + age | Subject, data = 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) # Plotting
fm1 <- lmer(
    formula = distance ~ age*Sex + (age|Subject)
    , data = Orthodont
newdat <- expand.grid(
    , 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

	[[alternative HTML version deleted]]

More information about the R-sig-mixed-models mailing list