[R] Help with response CIs for lme

[Spencer Graves]

      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.
>