[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
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]]
More information about the R-sig-mixed-models
mailing list