How can I get prediction intervals from a mixed-effects model?
Consider the following example:
library(nlme)
fm3 <- lme(distance ~ age*Sex, data = Orthodont, random = ~ 1)
df3.1 <- with(Orthodont, data.frame(age=seq(5, 20, 5),
Subject=rep(Subject[1], 4),
Sex=rep(Sex[1], 4)))
predict(fm3, df3.1, interval='prediction')
# M01 M01 M01 M01
# 22.69012 26.61199 30.53387 34.45574
# NOTE: The 'interval' argument to the 'predict' function was ignored.
# It works works for an 'lm' object, but not an 'lme' object.
One way to do this might be via mcmcsamp of the corresponding
'lmer' model:
library(lme4)
set.seed(3)
samp3r <- mcmcsamp(fm3r, n=10000)
samp3r[1:2,]
Then use library(coda) to check convergence and write a function
to simulate a single observation from each set of simulated parameters
and compute quantile(..., c(.025, .975)) for each prediction level
desired.
However, before I coded that, I thought I would ask if some easier
method might be available.
Thanks,
Spencer
p.s. RSiteSearch("lme prediction intervals") produced 3 hits including
2 from James A Rogers over 3 years ago. In one, he said, "I am not
aware of any published R function that gives you prediction intervals or
tolerance intervals for lme models."
(http://finzi.psych.upenn.edu/R/Rhelp02a/archive/42781.html) In the
other, he provided sample code for prediction or tolerance intervals of
lme variance components.
(http://finzi.psych.upenn.edu/R/Rhelp02a/archive/44675.html)
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
