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
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.merMod
dokumentację, 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 bootMer
do tego zadania”, ale na szczęście, nie widzę jak bootMer
tego 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.0
IN 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ę bootMer
zastosowania 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ć arm
pakiet.
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 lattice
pakietu:
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?
predictInterval
obejmuje błąd / niepewność zarówno dla stałych, jak i losowych warunków efektu. Wdotplot
widzisz 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,fishWt
a 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ć wpredictInterval
wynikach.