[R-sig-ME] lme and prediction intervals

[Ben Bolker]

  As I believe someone else commented, these seem to be confidence
intervals (conditional on random effects estimates) and not prediction
intervals.  To get prediction intervals (I think) you would need to
incorporate the residual variance: see below.

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)

library(ggplot2)
pd <- position_dodge(width=0.4)
ggplot(newdat,aes(x=age,y=pred,colour=Sex))+
  geom_point(position=pd)+
  geom_linerange(aes(ymin=pred-2*SE,ymax=pred+2*SE),
                 position=pd)

## prediction intervals
ggplot(newdat,aes(x=age,y=pred,colour=Sex))+
  geom_point(position=pd)+
  geom_linerange(aes(ymin=pred-2*SE2,ymax=pred+2*SE2),
                 position=pd)


D Chaws wrote:
> Dr. Bates, would you be able to weigh in on this?
> 
> On Thu, Feb 18, 2010 at 12:25 PM, D Chaws <cat.dev.urandom at gmail.com> wrote:
>> Dear lme(r) users,
>>
>> I have seen this issue come up many times over the years, but haven't
>> come across an answer as of yet.
>> Take for example the growth model:
>>
>> fm1 <- lme(distance ~ age*Sex, random = ~ 1 + age | Subject, data = Orthodont)
>>
>> Distance increases with age, but more so in Males.
>>
>> I want to obtain the model predicted values of distance at each age
>> (c(8,10,12,14)) for males and female separately to explore this
>> interaction.
>>
>> So,
>>
>> newdat <- expand.grid(age=c(8,10,12,14), Sex=c("Male","Female"))
>> newdat$pred <- predict(fm1, newdat, level = 0)
>>
>> R# newdat
>>  age    Sex  pred
>> 1   8   Male 22.62
>> 2  10   Male 24.18
>> 3  12   Male 25.75
>> 4  14   Male 27.32
>> 5   8 Female 21.21
>> 6  10 Female 22.17
>> 7  12 Female 23.13
>> 8  14 Female 24.09
>>
>> Yup, males have a steeper increase with age than females.
>>
>> The question is, how to go about getting prediction intervals around
>> these predictions.  It seems reasonable to need to know
>> the precision of these predictions, and of course most journals
>> require the reporting of error bars etc...  However, predict.lme
>> doesn't
>> have a se.fit or intervals argument.
>>
>> The only answer I have found at the moment is to use the design matrix
>> and $varFix from the model.
>>
>> So,
>>
>> Designmat <- model.matrix(eval(eval(fm1$call$fixed)[-2]), newdat[-3])
>> R# Designmat
>>  (Intercept) age SexFemale age:SexFemale
>> 1           1   8         0             0
>> 2           1  10         0             0
>> 3           1  12         0             0
>> 4           1  14         0             0
>> 5           1   8         1             8
>> 6           1  10         1            10
>> 7           1  12         1            12
>> 8           1  14         1            14
>>
>> newdat$SE <- sqrt(diag(Designmat %*% fm1$varFix %*% t(Designmat)))
>> R# newdat
>>  age    Sex  pred     SE
>> 1   8   Male 22.62 0.5265
>> 2  10   Male 24.18 0.4848
>> 3  12   Male 25.75 0.5021
>> 4  14   Male 27.32 0.5730
>> 5   8 Female 21.21 0.6350
>> 6  10 Female 22.17 0.5847
>> 7  12 Female 23.13 0.6056
>> 8  14 Female 24.09 0.6910
>>
>> Are these true pointwise prediction intervals?
>>
>> Any help would be greatly appreciated.  I refuse to use SAS for this!
>>
>> -- D. Chaws
>>
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


-- 
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bolker at ufl.edu / people.biology.ufl.edu/bolker
GPG key: people.biology.ufl.edu/bolker/benbolker-publickey.asc