Extrahiere Vorhersage-Band von lme fit

Habe ich folgende Modell

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)

Extrahiere Vorhersage-Band von lme fit

mit predict und level=0 ich kann das Grundstück die Durchschnittliche Bevölkerung Reaktion. Wie kann ich extrahieren und Darstellung der 95% Konfidenzintervalle /Vorhersage bands aus dem nlme-Objekt für die ganze Bevölkerung?

Kommentar zu dem Problem
gute Frage! Wenn in verstehen, Sie versuchen, ein äquivalent von curve(predict(model.lm, data.frame(x=x),interval ='confidence'),add=T), wobei-Modell.lm e.g lm(y~x) Kommentarautor: agstudy
Ja. Mit den unteren und oberen GUS-Staaten. Kommentarautor: ECII
Ich denke, auch für einen Film, ist es eine lästige Pflicht, um es zu bekommen. es ist Intervalle .lme - Funktion, aber es nicht geben die band Vertrauen, das nur einen Punkt. Kommentarautor: agstudy
Intervalle ruft die GUS-Schätzungen/ - Koeffizienten der passt. Was ich brauche, den GUS-y für alle x gegeben. Kommentarautor: ECII
in der Tat , @ECII haben Sie versucht, den harten Weg ...ich meine Berechnung von der band selbst..? Kommentarautor: agstudy

InformationsquelleAutor der Frage ECII | 2013-01-16

  1. 21

    Warnung: Lesen diesen thread r-sig-mixed-Modelle bevor Sie dies tun. Seien Sie sehr vorsichtig bei der Interpretation der resultierenden Vorhersage band.

    Vom r-sig-mixed-Modelle FAQ angepasst auf dein Beispiel:

    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

    Extrahiere Vorhersage-Band von lme fit

    InformationsquelleAutor der Antwort Roland

Schreibe einen Kommentar