Interwał prognozy dla modelu efektów mieszanych lmer () w R

Chcę uzyskać przedział przewidywania wokół prognozy z modelu lmer (). Znalazłem trochę dyskusji na ten temat:

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

http://glmm.wikidot.com/faq

ale wydaje się, że nie uwzględniają niepewności losowych efektów.

Oto konkretny przykład. Ścigam się złotą rybką. Mam dane dotyczące ostatnich 100 wyścigów. Chcę przewidzieć 101., biorąc pod uwagę niepewność moich oszacowań RE i oszacowań FE. Włączam losowe przechwytywanie ryb (jest 10 różnych ryb) i ustalony efekt dla wagi (mniej ciężkie ryby są szybsze).

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)

Teraz, aby przewidzieć 101. wyścig. Ryby zostały zważone i są gotowe do wypłynięcia:

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

Ryba D naprawdę puściła się (1,11 uncji) i przewiduje się, że przegra z Ryba E i Ryba F, które były lepsze niż w przeszłości. Jednak teraz chcę móc powiedzieć: „Ryba E (o wadze 0,91 uncji) pokona rybę D (o wadze 1,11 uncji) z prawdopodobieństwem p”. Czy istnieje sposób na wykonanie takiego oświadczenia przy użyciu lme4? Chcę, aby moje prawdopodobieństwo p uwzględniało moją niepewność zarówno dla efektu ustalonego, jak i efektu losowego.

Dzięki!

PS patrząc na predict.merModdokumentację, sugeruje: „Nie ma możliwości obliczenia standardowych błędów prognoz, ponieważ trudno jest zdefiniować skuteczną metodę uwzględniającą niepewność w parametrach wariancji; zalecamy bootMerdo tego zadania”, ale na szczęście, nie widzę jak bootMertego dokonać. Wygląda na to, bootMerże zostanie wykorzystany do uzyskania przedziałów ufności ładowania początkowego dla oszacowań parametrów, ale mogę się mylić.

ZAKTUALIZOWANY P:

OK, myślę, że zadawałem złe pytanie. Chcę móc powiedzieć: „Ryba A, ważąca w oz, będzie miała czas wyścigu, który wynosi (lcl, ucl) w 90% przypadków”.

W przedstawionym przeze mnie przykładzie Ryba A, ważąca 1,0 uncja, będzie miała 9 + 0.1 + 1 = 10.1 secśredni czas wyścigu ze standardowym odchyleniem 0,1. Tak więc jego obserwowany czas wyścigu będzie pomiędzy

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% czasu. Chcę funkcji przewidywania, która próbuje dać mi tę odpowiedź. Ustawienie wszystkich fishWt = 1.0IN newDat, ponowne uruchomienie SIM, używając (jako sugerowane przez Ben Bolker poniżej)

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

daje

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

Wydaje się, że tak naprawdę koncentruje się wokół średniej populacji? Jakby nie uwzględniał efektu FishID? Pomyślałem, że może to problem z wielkością próby, ale kiedy podniosłem liczbę obserwowanych ras od 100 do 10000, nadal otrzymuję podobne wyniki.

Domyślnie odnotuję bootMerzastosowania use.u=FALSE. Z drugiej strony, używając

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

daje

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

Ten przedział jest zbyt wąski i wydaje się być przedziałem ufności dla średniego czasu Ryby A. Chcę przedziału ufności dla obserwowanego czasu wyścigu Ryb A, a nie jego średniego czasu wyścigu. Jak mogę to zdobyć?

AKTUALIZACJA 2, PRAWIE:

Myślałem, że znalazłem to, czego szukałem w Gelman i Hill (2007) , strona 273. Potrzebuję wykorzystać armpakiet.

library("arm")

W przypadku ryby 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 

Dla wszystkich ryb:

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 

Właściwie to prawdopodobnie nie jest dokładnie to, czego chcę. Biorę tylko pod uwagę ogólną niepewność modelu. W sytuacji, gdy mam, powiedzmy, 5 zaobserwowanych wyścigów dla Ryby K i 1000 obserwowanych wyścigów dla Ryby L, myślę, że niepewność związana z moją prognozą dla Ryb K powinna być znacznie większa niż niepewność związana z moją prognozą dla Ryb L.

Przyjrzymy się bliżej Gelmanowi i Hillowi 2007. Wydaje mi się, że mogę w końcu przejść na BŁĘDY (lub Stan).

AKTUALIZACJA 3:

Być może źle sobie wyobrażam. Użycie predictInterval()funkcji podanej przez Jareda Knowlesa w poniższej odpowiedzi daje przedziały, które nie są dokładnie takie, jakich bym się spodziewał ...

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))

Dodałem dwie nowe ryby. Ryba K, dla której zaobserwowaliśmy 995 ras, i Ryba L, dla których zaobserwowaliśmy 5 ras. Obserwowaliśmy 100 wyścigów dla Fish AJ. Pasuję tak samo lmer()jak poprzednio. Patrząc na dotplot()z latticepakietu:

Domyślnie dotplot()porządkuje losowe efekty według ich oszacowania punktowego. Szacunek dla ryby L znajduje się w górnej linii i ma bardzo szeroki przedział ufności. Ryba K znajduje się na trzeciej linii i ma bardzo wąski przedział ufności. To ma dla mnie sens. Mamy wiele danych na temat Fish K, ale nie ma wielu danych na temat Fish L, więc jesteśmy bardziej pewni naszego szacunku na temat prawdziwej prędkości pływania Fish K. Teraz sądzę, że doprowadziłoby to do wąskiego przedziału prognoz dla ryby K i szerokiego przedziału prognoz dla ryby L podczas używania 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()

Wszystkie przedziały prognozowania wydają się mieć identyczną szerokość. Dlaczego nasze prognozy dotyczące Fish K nie są węższe od pozostałych? Dlaczego nasza prognoza dla Fish L nie jest szersza niż inne?

To pytanie i doskonała wymiana była impulsem do stworzenia predictIntervalfunkcji w merToolspakiecie. bootMerjest droga, ale w przypadku niektórych problemów obliczeniowo nie jest możliwe wygenerowanie przerwań całego modelu (w przypadkach, gdy model jest duży).

W takich przypadkach predictIntervaljest przeznaczony do wykorzystania arm::simfunkcji do generowania rozkładów parametrów w modelu, a następnie do wykorzystania tych rozkładów do wygenerowania symulowanych wartości odpowiedzi podanej newdataprzez użytkownika. Jest prosty w użyciu - wszystko, co musisz zrobić, to:

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

Możesz określić cały szereg innych wartości, w predictIntervaltym ustawienie interwału dla przedziałów predykcji, wybranie, czy zgłosić średnią lub medianę rozkładu oraz wybranie, czy dołączyć resztową wariancję z modelu.

Nie jest to pełny przedział predykcji, ponieważ zmienność thetaparametrów w lmerobiekcie nie jest uwzględniona, ale wszystkie inne zmiany są rejestrowane za pomocą tej metody, co daje całkiem przyzwoite przybliżenie.

Zrób to, bootMergenerując zestaw prognoz dla każdej replikacji parametrycznego ładowania początkowego:

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

Dane wyjściowe bootMerznajdują się w niezbyt przezroczystym "boot"obiekcie, ale możemy uzyskać surowe prognozy z $tkomponentu.

Ile czasu Fish E pokonuje Fish D?

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

Czasy ryb E są w kolumnie 5, czasy ryb D są w kolumnie 4, więc musimy tylko znać proporcję, że kolumna 5 jest mniejsza niż kolumna 4:

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