[R-sig-ME] lme and prediction intervals

[D Chaws]

Ok, issue solved for the most straightforward random effects cases.
Not sure about nested random effects or more complex cases.

vcov(fm1) will give you the appropriate covariance matrix that can
then be applied to the design matrix to get SEs.

Brian S. Yandell has a nice version of lsmean for his pda package that
handles prediction SEs
well, and which appear match very closely to those computed from SAS.

The link is:

http://r-forge.r-project.org/plugins/scmsvn/viewcvs.php/*checkout*/pkg/R/lsmean.R?rev=2&root=pda

-- DC

On Fri, Feb 26, 2010 at 12:32 AM, D Chaws <cat.dev.urandom at gmail.com> wrote:
> 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
>>
>