[R-sig-ME] lme and prediction intervals

[D Chaws]

And the saga continues...

I checked on the SAS website to see how they compute standard errors
for predictions in LSMEANS, just for giggles.  According to
http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/statug_mixed_sect014.htm

SE = L(X'V-1X)-L'

or in R terms

SE = Designmat %*% solve(formXtViX(extract.lme.cov2(fm1, Orthodont),
extract.lmeDesign(fm1)$X), t(Designmat))

which is, yes...you guessed it, the exact same thing as Designmat %*%
fm1$varFix %*% t(Designmat).

And of course, the SE's don't match those in LSMEANS.  Any thoughts?
I can't believe a solution to this has never come up.  I don't care
about replicating LSMEANS, I just want to be confident that the SE's I
present are accurate and meaningful.  Since they currently leave out
random-effects variance, I suspect something important is missing.

-- DC

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
>