R help - Apr 2008 - prediction intervals from a mixed-effects models?

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)

Spencer Graves <spencer.graves <at> pdf.com> writes:
> 
>       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.
> 
In theory, ci from gmodels should work

library(gmodels)
ci(df3.1)


However, I get an error message. I will forward this to Greg to let him check if
I did something stupid.

Dieter

Dear Reinhold: 

      Thanks for the suggestion.  The index of Gelman and Hill for 
"prediction" for "multilevel model" mentions pp. 272-275 and
361-363.
P. 362 starts with "predicting a new unit in an existing group", which
sounds like what I want.  Now all I need to do is study that book enough 
to be able to do what it says. 

      Thanks again. 
      Spencer

Reinhold Kliegl wrote:
> Spencer,
>
> I think the Gelman & Hill (2007) book has examples that look less
> complicated to me in comparison to what you describe (i.e. simply
> sample from the estimated distributions). I have some code for
> computation of power, also following examples in this book. Perhaps, I
> am overlooking something.
>
> Reinhold
>
> On Sun, Apr 13, 2008 at 7:10 PM, Spencer Graves <spencer.graves at
pdf.com> wrote:
>   
>>       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-sig-mixed-models at r-project.org mailing list
>>  https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>     
>
>