Extract prediction band from lme fit

23,344

Solution 1

Warning: Read this thread on r-sig-mixed models before doing this. Be very careful when you interpret the resulting prediction band.

From r-sig-mixed models FAQ adjusted to your example:

set.seed(42)
x <- rep(0:100,10)
y <- 15 + 2*rnorm(1010,10,4)*x + rnorm(1010,20,100)
id<-rep(1:10,each=101)

dtfr <- data.frame(x=x ,y=y, id=id)

library(nlme)

model.mx <- lme(y~x,random=~1+x|id,data=dtfr)

#create data.frame with new values for predictors
#more than one predictor is possible
new.dat <- data.frame(x=0:100)
#predict response
new.dat$pred <- predict(model.mx, newdata=new.dat,level=0)

#create design matrix
Designmat <- model.matrix(eval(eval(model.mx$call$fixed)[-2]), new.dat[-ncol(new.dat)])

#compute standard error for predictions
predvar <- diag(Designmat %*% model.mx$varFix %*% t(Designmat))
new.dat$SE <- sqrt(predvar) 
new.dat$SE2 <- sqrt(predvar+model.mx$sigma^2)

library(ggplot2) 
p1 <- ggplot(new.dat,aes(x=x,y=pred)) + 
geom_line() +
geom_ribbon(aes(ymin=pred-2*SE2,ymax=pred+2*SE2),alpha=0.2,fill="red") +
geom_ribbon(aes(ymin=pred-2*SE,ymax=pred+2*SE),alpha=0.2,fill="blue") +
geom_point(data=dtfr,aes(x=x,y=y)) +
scale_y_continuous("y")
p1

plot

Solution 2

Sorry for coming back to such an old topic, but this might address a comment here:

it would be nice if some package could provide this functionality

This functionality is included in the ggeffects-package, when you use type = "re" (which will then include the random effect variances, not only residual variances, which is - however - the same in this particular example).

library(nlme)
library(ggeffects)

x  <- rep(seq(0, 100, by = 1), 10)
y  <- 15 + 2 * rnorm(1010, 10, 4) * x + rnorm(1010, 20, 100)
id <- NULL
for (i in 1:10) {
  id <- c(id, rep(i, 101))
}
dtfr <- data.frame(x = x, y = y, id = id)
m <- lme(y ~ x,
         random =  ~ 1 + x | id,
         data = dtfr,
         na.action = na.omit)

ggpredict(m, "x") %>% plot(rawdata = T, dot.alpha = 0.2)

ggpredict(m, "x", type = "re") %>% plot(rawdata = T, dot.alpha = 0.2)

Created on 2019-06-18 by the reprex package (v0.3.0)

Share:
23,344

Related videos on Youtube

ECII
Author by

ECII

Updated on June 18, 2020

Comments

  • ECII
    ECII over 3 years

    I have following model

    x  <- rep(seq(0, 100, by=1), 10)
    y  <- 15 + 2*rnorm(1010, 10, 4)*x + rnorm(1010, 20, 100)
    id <- NULL
    for(i in 1:10){ id <- c(id, rep(i,101)) }
    dtfr <- data.frame(x=x,y=y, id=id)
    library(nlme)
    with(dtfr, summary(     lme(y~x, random=~1+x|id, na.action=na.omit)))
    model.mx <- with(dtfr, (lme(y~x, random=~1+x|id, na.action=na.omit)))
    pd       <- predict( model.mx, newdata=data.frame(x=0:100), level=0)
    with(dtfr, plot(x, y))
    lines(0:100, predict(model.mx, newdata=data.frame(x=0:100), level=0), col="darkred", lwd=7)
    

    enter image description here

    with predict and level=0 i can plot the mean population response. How can I extract and plot the 95% confidence intervals / prediction bands from the nlme object for the whole population?

    • agstudy
      agstudy almost 11 years
      good question! If in understand you try to have an equivalent of this curve(predict(model.lm, data.frame(x=x),interval ='confidence'),add=T) where model.lm e.g is lm(y~x)
    • ECII
      ECII almost 11 years
      Yes. With the lower and upper CIs.
    • agstudy
      agstudy almost 11 years
      I think even for a lm it is a chore to get it. there is intervals .lme function but it don't give the band confidence just one point.
    • ECII
      ECII almost 11 years
      intervals gets the CIs of the estimates/coefficients of the fits. What I need the CIs of y for any given x.
    • agstudy
      agstudy almost 11 years
      in fact , @ECII did you try the hard way ...I mean calculating the band by yourself..?
    • ECII
      ECII almost 11 years
      What do you mean "by myself"?
    • Roland
      Roland almost 11 years
  • Ben Bolker
    Ben Bolker almost 11 years
    +1 for saving me the trouble of doing it or feeling guilty for not doing it. Do note (as commented in the FAQ) that this only accounts for the uncertainty of the fixed effects conditional on the estimates of the random-effect variances and BLUPs/conditional modes
  • Dieter Menne
    Dieter Menne almost 11 years
    It would be good to have a cross-reference from you FAQ to here. I remember having re-invented that wheel quite a few time.
  • Roland
    Roland almost 11 years
    I don't think Prof. Bates would include this in predict.lme, but it would be nice if some package could provide this functionality (of course with clear warnings concerning the statistical limitations in the help file).
  • ECII
    ECII almost 11 years
    Would any of you mind explaining a bit what happens in the newdat ? Is this code applicable for any number of regressors?
  • Dieter Menne
    Dieter Menne almost 11 years
    @Roland: I fear that Douglas Bates would switch into his critical mode and ask "What for? Are you sure your non-statistical customers correctly read the graph?". "Conditional on the estimates" is a damned ugly concept to explain to clinical researchers.