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


37

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:

Szacunki FishID

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

Interwał prognozy dla ryb

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?


1
predictIntervalobejmuje błąd / niepewność zarówno dla stałych, jak i losowych warunków efektu. W dotplotwidzisz tylko niepewność ze względu na losowy części przepowiedni zasadniczo niepewności wokół szacunków ryb określonych przechwytuje. Jeśli twój model ma dużo niepewności w stałym parametrze, fishWta ten parametr steruje większością przewidywanej wartości, to niepewność wokół każdego konkretnego przechwytu ryby jest banalna i nie zobaczysz dużej różnicy w szerokości interwałów. Powinniśmy to wyjaśnić w predictIntervalwynikach.
jknowles

Odpowiedzi:


18

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.


3
To wygląda niesamowicie! Czytanie przez winiet teraz. Dzięki!
możliwe

Przedziały prognozowania nie są dokładnie takie, jak się spodziewałem. Zobacz aktualizację 3 powyżej.
możliwe

Nie predictInterval()lubi zagnieżdżonych efektów losowych? Na przykład przy użyciu msleepzestawu danych z ggplot2pakietu: mod <- lmer(sleep_total ~ bodywt + (1|vore/order), data=msleep); predInt <- predictInterval(merMod=mod, newdata=msleep) Zwraca błąd:Error in '[.data.frame'(newdata, , j) : undefined columns selected
możliwe

Założę się, że nie lubi efektów zagnieżdżonych. Nie sądzę, żebyśmy mieli jakieś testy w naszym pakiecie testowym. Zgłoszę problem na GitHub, aby to sprawdzić. Polecam również wypróbowanie devtools::install_github("jknowles/merTools")najpierw wersji deweloperskiej z GitHub .
jknowles

2
Jako aktualizacja, najnowsza wersja programistyczna merTools pozwala na zagnieżdżone efekty. Wkrótce zostanie przekazany do CRAN.
jknowles

15

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

Otrzymuję nieoczekiwane wyniki. Gdybym ustawił fishWt = 1 dla wszystkich ryb w newDat, oczekiwałbym średniego / mediany czasu dla ryby A wynoszącej ~ 10,1, ryby B ~ 10,2, ..., ryby J ~ 11,0 (ponieważ ich czas w danych treningowych wynosi zdefiniowane jako:) sec <- 9 + rep(1:10, rep(100,10))/10 + oz + rnorm(1000)/10. Gdy używam predict(), czasy przewidywania dla ryb A, E i J wynoszą 10,09, 10,49 i 10,99, zgodnie z oczekiwaniami. Jednak mediana czasów dla opisanej metody bootMer wynosi: 10,52, 10,59 i 10,50. Spodziewałbym się więcej zgody?
możliwe

Używanie use.u=TRUEjak w: bb <- bootMer(lme1,nsim=200,FUN=predFun,seed=101,use.u=TRUE)wydaje się dać mi to, czego chcę. Dzięki!
możliwe

OK, to staje się trochę trudne. Musisz spojrzeć na use.uargument do bootMer. Pytanie brzmi: kiedy mówisz „niepewność co do efektu stałego i efektu losowego”, co rozumiesz przez „efekt losowy”? Czy masz na myśli niepewność co do wariancji efektów losowych lub trybów warunkowych (tj. Efektów specyficznych dla ryb)? Możesz użyć use.u=TRUE, ale nie sądzę, żeby to
zrobiło,

Jeśli użyję use.u=TRUE, to „wartości u [zostań] ustalone na ich wartości szacunkowe”. Interpretuję to jako znaczenie, niezależnie od tego, jaki jest nasz szacunkowy punkt losowy dla Ryby A, jest to traktowane jako Boska Uczciwa Prawda, jeśli wolisz. bootMerzakłada, że ​​nie ma błędu w naszym oszacowaniu punktu RE. Jeśli używam use.u=FALSE, czy bootMerw ogóle bierze pod uwagę oszacowania punktu RE? Wydaje się, że bootMerwyniki przy użyciu use.u=FALSEsą równoważne (lub asymptotycznie równoważne) z użyciem re.form=NAw predict()instrukcji. Czy to prawda?
możliwe

1
Myślę, że nie jest zaimplementowany ATM, ale możesz wyodrębnić wariancje warunkowe trybów warunkowych / BLUP poprzez c(attr(ranef(lme1,condVar=TRUE)[[1]],"postVar"))(wszystkie są identyczne w tym przykładzie), a następnie przetestować te wartości.
Ben Bolker
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.