[R-sig-ME] prediction intervals from a mixed-effects models?

[Nicholas Lewin-Koh]

Hi,
I think Doug made a comment a few weeks ago about why there was as
yet not a predict method for lme4, to the effect that there are so
many different things that can be extracted from the model. And
I have searched before for code for prediction intervals with a 
similar outcome. I think the mcmcsamp function offers a really
nice way to get prediction intervals, since effectively any function
of the parameters can be applied to the mcmcsample and get a valid
prediction (credible set) interval (am I remembering my Bayesian 
statistics correctly, it is the joint distribution right?) 

I have been toying with writing a general predict method for
lmer objects, given that Doug has elaborated quite a bit on the
underlying mer structures and algorithms. I suspect that 
this won't be as easy for glmm's, or is it again just a
matter of twiddling with the joint posterior distribution 
of the parameters? 

Doug would you accept a predict method for lme4, or should
that go into a separate package. I am also thinking
of tackling splines, but that won't be next week.

Nicholas

> Message: 2
> Date: Sun, 13 Apr 2008 10:10:50 -0700
> From: Spencer Graves <spencer.graves at pdf.com>
> Subject: [R-sig-ME] prediction intervals from a mixed-effects models?
> To: r-sig-mixed-models at r-project.org, r-help at r-project.org
> Message-ID: <48023E9A.5060606 at pdf.com>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
> 
>       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)
> 
> 
>