[R] lme - Confidence interval for predicted values
Douglas Bates
bates at stat.wisc.edu
Tue Apr 27 15:32:13 CEST 2004
"Urs Wiedemann" <wiedemann at fmp-berlin.de> writes:
> After having fit (using lme) a mixed effects model with a single
> random effect, I like to estimate the confidence interval for the
> predicted mean expectations.
> To my knowledge this is usually done (for fixed effects models) by
> calculating:
>
> cibandwidth <- sqrt(diag(Xnew %*% solve(t(X) %*% X) t(Xnew))) *
> qt((1-level)/2, DF)
>
> The CI is then the predicted value +/- cibandwidth. This is what
> predict.lm provides with ci.fit if I am not wrong.
>
> So for lme there is a very nice post on the S-news list:
> http://www.biostat.wustl.edu/archives/html/s-news/2003-09/msg00021.html
>
> Hopefully I quote the message correctly:
> >var(y.hat) = X2 %*% inverse(t(X)%*%inverse(Sig)%*%X) %*% t(X2)
> >
> >The hard part for lme is deciding what goes in X and what is Sig:
> >If you want a confidence interval on y for a particular subject
> >included in X, then that subject is part of X; otherwise,
> >it must be included in Sig.
I think the answer to this question will soon become easier. Later
today (if I can escape administrative duties, otherwise tomorrow) I
will make available an alpha-level release candidate of a completely
rewritten lme4 package. And I do mean completely rewritten. This is
the fourth and, I swear, the last time I have designed and implemented
from scratch methods for parameter estimation in linear mixed models.
The new version of lme (that is, from the not-yet-released lme4_0.6-1
or later) allows the usual optional arguments for model-fitting
functions in R, including "x=TRUE" which will cause the slot x in the
fitted model object to contain a list of model matrices. The last
element of this list is the X that you want.
> My question is now where to obtain "Sig" from lme. Probably this is
> obvious and I simply overlooked it.
The new version of the lme4 package has a method for the "vcov"
generic. This method returns this Sig.
More information about the R-help
mailing list