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@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. > ______________________________________________ R-help@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.