Extract prediction band from lme fit
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
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)
Related videos on Youtube
ECII
Updated on June 18, 2020Comments
-
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)
with
predict
andlevel=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 almost 11 yearsgood 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 almost 11 yearsYes. With the lower and upper CIs.
-
agstudy almost 11 yearsI 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 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 almost 11 yearsin fact , @ECII did you try the hard way ...I mean calculating the band by yourself..?
-
ECII almost 11 yearsWhat do you mean "by myself"?
-
Roland almost 11 yearsRead the r-sig-mixed models FAQ.
-
-
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 almost 11 yearsIt 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 almost 11 yearsI 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 almost 11 yearsWould any of you mind explaining a bit what happens in the
newdat
? Is this code applicable for any number of regressors? -
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.