Quiero obtener una predicción del intervalo de alrededor de una predicción de un lmer() del modelo. He encontrado alguna discusión acerca de esto:
http://rstudio-pubs-static.s3.amazonaws.com/24365_2803ab8299934e888a60e7b16113f619.html
pero parecen no tomar a la incertidumbre de los efectos aleatorios en cuenta.
He aquí un ejemplo específico. Yo soy de las carreras de pez de oro. Tengo los datos de los últimos 100 carreras. Quiero predecir el número 101, teniendo en cuenta la incertidumbre de mi RE estimaciones, y la FE de las estimaciones. Estoy incluyendo una al azar para interceptar peces (hay 10 tipos diferentes de peces), y el efecto fijo de peso (menos pesados que los peces son más rápidos).
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)
Ahora, para predecir el número 101 de la carrera. Los peces han sido pesadas y está listo para ir:
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
Los peces D realmente ha soltó (1.11 oz) y es realmente predijo a perder a los Peces E y de Pescado F, ambos de los cuales ha sido mejor que en el pasado. Sin embargo, ahora quiero ser capaz de decir, "Peces E (pesaje 0.91 oz) va a vencer a los Peces D (pesaje 1.11 oz) con una probabilidad de p". Hay una manera de hacer una declaración utilizando lme4? Quiero que mi probabilidad p para tener en cuenta mi incertidumbre, tanto en el de efectos fijos y el de efectos aleatorios.
Gracias!
PS mirando a la predict.merMod documentación, sugiere que "no Hay ninguna opción para calcular los errores estándar de las predicciones debido a que es difícil definir un método eficiente que incorpora la incertidumbre en los parámetros de varianza; recomendamos bootMer para esta tarea," pero mira por dónde, no puedo ver cómo utilizar bootMer a ello. Parece bootMer sería utilizado para obtener bootstrap intervalos de confianza para las estimaciones de los parámetros, pero puedo estar equivocado.
ACTUALIZADO P:
OK, creo que me estaba haciendo la pregunta equivocada. Quiero ser capaz de decir, "Pescado, de peso w oz, tendrá un tiempo de carrera que es (lcl, ucl), el 90% del tiempo."
En el ejemplo que he puesto, Pescado, pesaje 1.0 oz, tendrá un tiempo de carrera de 9 + 0.1 + 1 = 10.1 sec en promedio, con una desviación estándar de 0.1. Por lo tanto, su observó el tiempo de carrera será de entre
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
El 90% del tiempo. Quiero una función de predicción de que los intentos de darme la respuesta. Configuración de todos los fishWt = 1.0 en newDat, de volver a ejecutar la sim, y usando (como sugerido por Ben Bolker a continuación)
bb <- bootMer(lme1,nsim=1000,FUN=predFun, use.u = FALSE)
da
> quantile(predMat[,1], c(0.05,0.50,0.95))
5% 50% 95%
10.01362 10.55646 11.05462
Esto parece realmente estar centrada en torno a la media de la población? Como si no fuera a tomar el FishID efecto en la cuenta? Pensé que tal vez se trataba de un tamaño de la muestra problema, pero cuando me encontré con el número de observados carreras de 100 a 10000, yo todavía obtener resultados similares.
Por otro lado, el uso de
bb <- bootMer(lme1,nsim=1000,FUN=predFun, use.u = TRUE)
da
> quantile(predMat[,1], c(0.05,0.50,0.95))
5% 50% 95%
10.09970 10.10128 10.10270
Ese intervalo es demasiado estrecho, y parecería ser un intervalo de confianza para Peces de Una media hora. Quiero un intervalo de confianza para Un Pescado que se observa el tiempo de carrera, no es su promedio de carrera de tiempo. Cómo puedo obtener?
ACTUALIZAR 2, CASI:
Yo pensaba que encontré lo que estaba buscando, en Gelman y Hill (2007) , página 273. Necesidad de utilizar la arm paquete.
library("arm")
Para El Pescado:
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
Para todos los peces:
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
En realidad, esto probablemente no es exactamente lo que quiero. Sólo estoy tomando en cuenta el conjunto de incertidumbre del modelo. En una situación en la que yo tengo, digamos, 5 observa las razas de Peces K y 1000 observado las razas de Peces L, creo que la incertidumbre asociada con mi predicción para los Peces K debe ser mucho mayor que la incertidumbre asociada con mi predicción para los Peces L.
Tendrá un aspecto más dentro de Gelman y Hill, 2007. Siento que puede terminar encima de tener que cambiar a los ERRORES (o Stan).
ACTUALIZACIÓN EL 3 de:
Tal vez estoy conceptualizar las cosas mal. El uso de la predictInterval() función dada por Jared Knowles en una respuesta a continuación da a intervalos que no son lo que yo esperaría...
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))
He añadido dos nuevos peces. Peces de K, para los cuales hemos observado 995 carreras, y los Peces L, para los cuales hemos observado 5 carreras. Hemos observado 100 razas de Peces a-J. I ajuste de la misma lmer() como antes. Mirando a la dotplot() de la lattice paquete de:
De forma predeterminada, dotplot() reordena el de efectos aleatorios, por su punto de estimación. La estimación para el Pescado L está en la parte superior de la línea, y tiene un ancho de intervalo de confianza. Los peces K es en la tercera línea, y tiene un muy estrecho intervalo de confianza. Esto tiene sentido para mí. Tenemos un montón de datos sobre Peces K, pero no una gran cantidad de datos sobre Peces de L, por lo que nos sentimos más seguros en nuestro conjeturar acerca de los Peces K la verdadera velocidad de natación. Ahora, yo creo que esto llevaría a un estrecho intervalo de predicción para los Peces K, y un amplio intervalo de predicción para los Peces L cuando se utiliza 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()
Todos los intervalos de predicción parece ser idéntica en anchura. Por qué no es nuestra predicción para los Peces K más estrecho que los otros? Por qué no es nuestra predicción para los Peces de L más amplia que la de los demás?
