[R-sig-ME] lme and prediction intervals

D Chaws cat.dev.urandom at gmail.com
Thu Feb 18 18:25:32 CET 2010


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




More information about the R-sig-mixed-models mailing list