Estratto previsione band di LME adatta

Question

Ho seguente modello

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) 

con predict e level=0 posso tracciare la risposta della popolazione media. Come posso estrarre e tracciare gli intervalli di confidenza al 95%/le bande di predizione dall'oggetto nlme per l'intera popolazione?

Answer 1

Avviso: Leggere this thread on r-sig-mixed models prima di procedere. Fai molta attenzione quando interpreti la banda di predizione risultante.

Da r-sig-mixed models FAQ regolata per il tuo esempio:

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