No, your example computation does NOT produce the desired "lme
prediction intervals". I ran your script and got exactly the same
numbers for upper and lower limits. Even without any consideration of
statistical theory, this suggests either shockingly precise statistical
estimation or a problem with your algorithm.
The theory behind such intervals is sufficiently complicated that
the 'nlme' developers have not found a need sufficient to justify the
effort required to develop the required capability. A reply by James
Rogers to a similar question a couple of years ago concluded, "It is not
easy to write such a function for the general case, but it may be
relatively easy to write your own for special cases of lme models."
(http://finzi.psych.upenn.edu/R/Rhelp02a/archive/42781.html)
To briefly outline some of the difficulties, we first should ask
if you want confidence intervals for the MEAN of future values or for
the future values themselves? And how do you want to handle the random
effects? If you want the mean of the fixed effects, averaging over
random effects, that should be fairly easy. Consider, for example, the
following modification of the example on the 'lme' help page:
fm1 <- lme(distance ~ age, data = Orthodont) # random is ~ age
summary(fm1)
<snip>
Fixed effects: distance ~ age
Value Std.Error DF t-value p-value
(Intercept) 16.761111 0.7752461 80 21.620375 0
age 0.660185 0.0712533 80 9.265334 0
Correlation:
(Intr)
age -0.848
From the "Std.Error" and "Correlation", we could
reconstruct the
covariance matrix of the parameter estimates, b.hat. Call this S.b.
Then the estimate for Ey for some new set of covariates coded in a row
vector x' is var(Ey.hat) = x' S.b x. The square root of this number is
a standard deviation, and you could add and subtract some multiple like
2 of this number from x' b.hat to get the desired confidence interval.
If I wanted to do that myself, I might read the code for
"summary.lme", after using getAnywhere("summary.lme") to get
it.
I know this doesn't solve your problem, but I hope it helps.
Spencer Graves
Michael Kubovy wrote:> Hi,
>
> Can someone please offer a procedure for going from CIs produced by
> intervals.lme() to fixed-effects response CIs.
>
> Here's a simple example:
>
> library(mlmRev)
> library(nlme)
> hsb.lme <- lme(mAch ~ minrty * sector, random = ~ 1 | cses, Hsb82)
> (intervals(hsb.lme))
> (hsb.new <- data.frame
> minrty = rep(c('No', 'Yes'), 2),
> sector = rep(c('Public', 'Catholic'), each = 2)))
> cbind(hsb.new, predict(hsb.lme, hsb.new, level = 0))
>
> Is the following correct (I know from the previous command that the
> estimate is correct)?
> cbind(hsb.new, rbind(hsb.int[[1]][1,], hsb.int[[1]][1,]+hsb.int[[1]]
> [2,], hsb.int[[1]][1,]+hsb.int[[1]][3,], hsb.int[[1]][1,]+hsb.int[[1]]
> [2,] + hsb.int[[1]][3,] + hsb.int[[1]][4,]))
>
> If so, is there an easier way to write it?
> _____________________________
> Professor Michael Kubovy
> University of Virginia
> Department of Psychology
> USPS: P.O.Box 400400 Charlottesville, VA 22904-4400
> Parcels: Room 102 Gilmer Hall
> McCormick Road Charlottesville, VA 22903
> Office: B011 +1-434-982-4729
> Lab: B019 +1-434-982-4751
> Fax: +1-434-982-4766
> WWW: http://www.people.virginia.edu/~mk9y/
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>