Vorhersageintervall für lmer () -Mischeffektmodell in R

Ich möchte ein Vorhersageintervall für eine Vorhersage aus einem lmer () -Modell erhalten. Ich habe eine Diskussion darüber gefunden:

http://rstudio-pubs-static.s3.amazonaws.com/24365_2803ab8299934e888a60e7b16113f619.html

http://glmm.wikidot.com/faq

Sie scheinen jedoch die Unsicherheit der zufälligen Effekte nicht zu berücksichtigen.

Hier ist ein konkretes Beispiel. Ich rase Goldfisch. Ich habe Daten zu den letzten 100 Rennen. Ich möchte den 101. vorhersagen, wobei ich die Unsicherheit meiner RE-Schätzungen und FE-Schätzungen berücksichtige. Ich schließe einen zufälligen Abschnitt für Fische ein (es gibt 10 verschiedene Fische) und einen festen Effekt für das Gewicht (weniger schwere Fische sind schneller).

library("lme4")

fish <- as.factor(rep(letters[1:10], each=100))
race <- as.factor(rep(900:999, 10))
oz <- round(1 + rnorm(1000)/10, 3)
sec <- 9 + rep(1:10, rep(100,10))/10 + oz + rnorm(1000)/10

fishDat <- data.frame(fishID = fish, 
      raceID = race, fishWt = oz, time = sec)
head(fishDat)
plot(fishDat$fishID, fishDat$time)

lme1 <- lmer(time ~ fishWt + (1 | fishID), data=fishDat)
summary(lme1)

Nun, um das 101. Rennen vorherzusagen. Die Fische wurden gewogen und sind bereit zu gehen:

newDat <- data.frame(fishID = letters[1:10], 
    raceID = rep(1000, 10),
    fishWt = 1 + round(rnorm(10)/10, 3))
newDat$pred <- predict(lme1, newDat)
newDat

   fishID raceID fishWt     pred
1       a   1000  1.073 10.15348
2       b   1000  1.001 10.20107
3       c   1000  0.945 10.25978
4       d   1000  1.110 10.51753
5       e   1000  0.910 10.41511
6       f   1000  0.848 10.44547
7       g   1000  0.991 10.68678
8       h   1000  0.737 10.56929
9       i   1000  0.993 10.89564
10      j   1000  0.649 10.65480

Fisch D hat sich wirklich losgelassen (1.11 oz) und wird voraussichtlich gegen Fisch E und Fisch F verlieren, die beide besser waren als in der Vergangenheit. Jetzt möchte ich jedoch sagen können: "Fisch E (mit einem Gewicht von 0,91 Unzen) schlägt Fisch D (mit einem Gewicht von 1,11 Unzen) mit der Wahrscheinlichkeit p." Gibt es eine Möglichkeit, mit lme4 eine solche Aussage zu treffen? Ich möchte, dass meine Wahrscheinlichkeit p meine Unsicherheit sowohl im festen als auch im zufälligen Effekt berücksichtigt.

Vielen Dank!

Wenn man sich die predict.merModDokumentation ansieht , schlägt dies vor: "Es gibt keine Möglichkeit, Standardfehler von Vorhersagen zu berechnen, da es schwierig ist, eine effiziente Methode zu definieren, die Unsicherheiten in den Varianzparametern berücksichtigt. Wir empfehlen dies bootMerfür diese Aufgabe." wie man bootMerdas macht. Es scheint bootMerverwendet zu werden, um Bootstrap-Konfidenzintervalle für Parameterschätzungen zu erhalten, aber ich könnte mich irren.

AKTUALISIERT F:

OK, ich glaube, ich habe die falsche Frage gestellt. Ich möchte in der Lage sein zu sagen, "Fisch A, der w oz wiegt, wird eine Rennzeit haben, die (lcl, ucl) 90% der Zeit ist."

In dem Beispiel, das ich angelegt habe, hat Fisch A, der 1,0 Unzen wiegt, eine 9 + 0.1 + 1 = 10.1 secdurchschnittliche Rennzeit mit einer Standardabweichung von 0,1. Somit liegt seine beobachtete Rennzeit zwischen

x <- rnorm(mean = 10.1, sd = 0.1, n=10000)
quantile(x, c(0.05,0.50,0.95))
       5%       50%       95% 
 9.938541 10.100032 10.261243 

90% der Zeit. Ich möchte eine Vorhersagefunktion, die versucht, mir diese Antwort zu geben. Einstellen aller fishWt = 1.0in newDat, wieder laufen die SIM und mit (wie von Ben Bolker vorgeschlagen unten)

predFun <- function(fit) {
  predict(fit,newDat)
}
bb <- bootMer(lme1,nsim=1000,FUN=predFun, use.u = FALSE)
predMat <- bb$t

gibt

> quantile(predMat[,1], c(0.05,0.50,0.95))
      5%      50%      95% 
10.01362 10.55646 11.05462 

Dies scheint tatsächlich um den Bevölkerungsdurchschnitt zu zentrieren? Als würde der FishID-Effekt nicht berücksichtigt? Ich dachte, es handele sich möglicherweise um ein Problem mit der Stichprobengröße, aber als ich die Anzahl der beobachteten Rennen von 100 auf 10000 erhöhte, erhalte ich immer noch ähnliche Ergebnisse.

Ich werde standardmäßig bootMerVerwendungen vermerken use.u=FALSE. Auf der anderen Seite mit

bb <- bootMer(lme1,nsim=1000,FUN=predFun, use.u = TRUE)

gibt

> quantile(predMat[,1], c(0.05,0.50,0.95))
      5%      50%      95% 
10.09970 10.10128 10.10270 

Dieses Intervall ist zu eng und scheint ein Konfidenzintervall für die mittlere Zeit von Fisch A zu sein. Ich möchte ein Konfidenzintervall für die beobachtete Rennzeit von Fisch A, nicht für seine durchschnittliche Rennzeit. Wie kann ich das bekommen?

UPDATE 2, FAST:

Ich dachte, ich hätte in Gelman and Hill (2007) , Seite 273, gefunden, wonach ich gesucht habe. Ich muss das armPaket verwenden.

library("arm")

Für Fisch A:

x.tilde <- 1    #observed fishWt for new race
sigma.y.hat <- sigma.hat(lme1)$sigma$data        #get uncertainty estimate of our model
coef.hat <- as.matrix(coef(lme1)$fishID)[1,]    #get intercept (random) and fishWt (fixed) parameter estimates
y.tilde <- rnorm(1000, coef.hat %*% c(1, x.tilde), sigma.y.hat) #simulate
quantile (y.tilde, c(.05, .5, .95))

  5%       50%       95% 
 9.930695 10.100209 10.263551 

Für alle Fische:

x.tilde <- rep(1,10)  #assume all fish weight 1 oz
#x.tilde <- 1 + rnorm(10)/10  #alternatively, draw random weights as in original example
sigma.y.hat <- sigma.hat(lme1)$sigma$data
coef.hat <- as.matrix(coef(lme1)$fishID)
y.tilde <- matrix(rnorm(1000, coef.hat %*% matrix(c(rep(1,10), x.tilde), nrow = 2 , byrow = TRUE), sigma.y.hat), ncol = 10, byrow = TRUE)
quantile (y.tilde[,1], c(.05, .5, .95))
       5%       50%       95% 
 9.937138 10.102627 10.234616 

Eigentlich ist das wahrscheinlich nicht genau das, was ich will. Ich berücksichtige nur die allgemeine Modellunsicherheit. In einer Situation, in der ich beispielsweise 5 Rennen für Fish K und 1000 Rennen für Fish L beobachtet habe, sollte die mit meiner Vorhersage für Fish K verbundene Unsicherheit viel größer sein als die mit meiner Vorhersage für Fish L verbundene Unsicherheit.

Ich werde weiter auf Gelman und Hill 2007 eingehen. Ich habe das Gefühl, dass ich möglicherweise zu BUGS (oder Stan) wechseln muss.

AKTUALISIEREN SIE DEN 3.:

Vielleicht konzeptualisiere ich die Dinge schlecht. Die Verwendung der predictInterval()Funktion von Jared Knowles in einer der folgenden Antworten ergibt Intervalle, die nicht ganz so sind, wie ich es erwarten würde ...

library("lattice")
library("lme4")
library("ggplot2")

fish <- c(rep(letters[1:10], each = 100), rep("k", 995), rep("l", 5))
oz <- round(1 + rnorm(2000)/10, 3)
sec <- 9 + c(rep(1:10, each = 100)/10,rep(1.1, 995), rep(1.2, 5)) + oz + rnorm(2000)

fishDat <- data.frame(fishID = fish, fishWt = oz, time = sec)
dim(fishDat)
head(fishDat)
plot(fishDat$fishID, fishDat$time)

lme1 <- lmer(time ~ fishWt + (1 | fishID), data=fishDat)
summary(lme1)
dotplot(ranef(lme1, condVar = TRUE))

Ich habe zwei neue Fische hinzugefügt. Fisch K, für den wir 995 Rennen beobachtet haben, und Fisch L, für den wir 5 Rennen beobachtet haben. Wir haben 100 Rennen für Fish AJ beobachtet. Ich passe lmer()wie vorher. Betrachten Sie die dotplot()aus dem latticePaket:

dotplot()Ordnet die zufälligen Effekte standardmäßig nach ihrer Punktschätzung neu an. Die Schätzung für Fisch L befindet sich in der obersten Zeile und hat ein sehr breites Konfidenzintervall. Fisch K befindet sich in der dritten Zeile und hat ein sehr enges Konfidenzintervall. Das macht für mich Sinn. Wir haben viele Daten zu Fish K, aber nicht viele Daten zu Fish L, daher sind wir uns unserer Einschätzung der tatsächlichen Schwimmgeschwindigkeit von Fish K sicherer. Nun würde ich annehmen, dass dies zu einem engen Vorhersageintervall für Fisch K und einem weiten Vorhersageintervall für Fisch L bei der Verwendung führen würde predictInterval(). Howeva:

newDat <- data.frame(fishID = letters[1:12],
                     fishWt = 1)

preds <- predictInterval(lme1, newdata = newDat, n.sims = 999)
preds
ggplot(aes(x=letters[1:12], y=fit, ymin=lwr, ymax=upr), data=preds) +
  geom_point() + 
  geom_linerange() +
  labs(x="Index", y="Prediction w/ 95% PI") + theme_bw()

Alle diese Vorhersageintervalle scheinen in der Breite identisch zu sein. Warum ist unsere Vorhersage für Fish K nicht enger als für die anderen? Warum ist unsere Vorhersage für Fisch L nicht breiter als andere?

Diese Frage und der ausgezeichnete Austausch waren der Anstoß für die Erstellung der predictIntervalFunktion im merToolsPaket. bootMerist der richtige Weg, aber für einige Probleme ist es rechnerisch nicht machbar, Bootstrap-Überholungen des gesamten Modells zu generieren (in Fällen, in denen das Modell groß ist).

In diesen Fällen predictIntervalwird entworfen, um die arm::simFunktionen zu verwenden, um Verteilungen von Parametern im Modell zu generieren und diese Verteilungen dann zu verwenden, um simulierte Werte der Antwort zu generieren, die newdatavom Benutzer bereitgestellt werden. Es ist einfach zu bedienen - alles, was Sie tun müssten, ist:

library(merTools)
preds <- predictInterval(lme1, newdata = newDat, n.sims = 999)

Sie können eine ganze Reihe anderer Werte angeben predictInterval, um das Intervall für die Vorhersageintervalle festzulegen , auszuwählen, ob der Mittelwert oder der Median der Verteilung gemeldet werden soll, und um festzulegen, ob die Restvarianz aus dem Modell einbezogen werden soll oder nicht.

Es handelt sich nicht um ein vollständiges Vorhersageintervall, da die Variabilität der thetaParameter im lmerObjekt nicht berücksichtigt wird. Alle anderen Variationen werden jedoch mit dieser Methode erfasst, was eine ziemlich gute Annäherung ergibt.

Erstellen Sie dazu bootMereine Reihe von Vorhersagen für jedes parametrische Bootstrap-Replikat:

predFun <- function(fit) {
    predict(fit,newDat)
}
bb <- bootMer(lme1,nsim=200,FUN=predFun,seed=101)

Die Ausgabe von bootMerist in einem nicht schrecklich transparenten "boot"Objekt, aber wir können die Rohprognosen aus der $tKomponente herausholen.

Wie oft schlägt Fish E Fish D?

predMat <- bb$t
dim(predMat) ## 200 rows (PB reps) x 10 (predictions)

Die Zeiten von Fisch E sind in Spalte 5, die Zeiten von Fisch D sind in Spalte 4, daher müssen wir nur wissen, in welchem ​​Verhältnis Spalte 5 kleiner ist als Spalte 4:

mean(predMat[,5]<predMat[,4])  ## 0.57