Intervalo de predicción para el modelo de efectos mixtos lmer () en R

37

Quiero obtener un intervalo de predicción alrededor de una predicción de un modelo lmer (). He encontrado alguna discusión sobre esto:

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

http://glmm.wikidot.com/faq

pero parecen no tener en cuenta la incertidumbre de los efectos aleatorios.

Aquí hay un ejemplo específico. Estoy compitiendo con peces dorados. Tengo datos sobre las últimas 100 carreras. Quiero predecir el 101, teniendo en cuenta la incertidumbre de mis estimaciones de RE y las estimaciones de FE. Incluyo una intercepción aleatoria para peces (hay 10 peces diferentes) y un efecto fijo para el peso (los peces menos pesados ​​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 la carrera 101. Los peces han sido pesados ​​y están listos para salir:

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

Fish D realmente se ha dejado llevar (1.11 oz) y en realidad se predice que perderá ante Fish E y Fish F, quienes han sido mejores que en el pasado. Sin embargo, ahora quiero poder decir: "El pez E (con un peso de 0,91 oz) vencerá al pez D (con un peso de 1,11 oz) con una probabilidad p". ¿Hay alguna manera de hacer tal declaración usando lme4? Quiero que mi probabilidad p tenga en cuenta mi incertidumbre tanto en el efecto fijo como en el efecto aleatorio.

¡Gracias!

PD mirando la predict.merModdocumentación, sugiere "No hay opción para calcular errores estándar de predicciones porque es difícil definir un método eficiente que incorpore incertidumbre en los parámetros de varianza; recomendamos bootMerpara esta tarea", pero por Dios, no puedo ver cómo usar bootMerpara hacer esto. Parece bootMerque se usaría para obtener intervalos de confianza de arranque para las estimaciones de parámetros, pero podría estar equivocado.

Q ACTUALIZADO:

OK, creo que estaba haciendo la pregunta equivocada. Quiero poder decir: "El pescado A, que pese w oz, tendrá un tiempo de carrera que es (lcl, ucl) el 90% del tiempo".

En el ejemplo que he presentado, Fish A, con un peso de 1.0 oz, tendrá un tiempo de carrera 9 + 0.1 + 1 = 10.1 secpromedio, con una desviación estándar de 0.1. Por lo tanto, su tiempo de carrera observado será 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 

90% del tiempo Quiero una función de predicción que intente darme esa respuesta. Configurar todos fishWt = 1.0en newDat, volver a ejecutar el simulador, y el uso (como se sugiere por Ben Bolker abajo)

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

da

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

¿Esto parece centrarse realmente en el promedio de la población? ¿Como si no tuviera en cuenta el efecto FishID? Pensé que tal vez era un problema de tamaño de muestra, pero cuando superé el número de carreras observadas de 100 a 10000, todavía obtengo resultados similares.

Notaré bootMerusos use.u=FALSEpor defecto. Por otro lado, usando

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 parece ser un intervalo de confianza para el tiempo medio de Fish A. Quiero un intervalo de confianza para el tiempo de carrera observado de Fish A, no su tiempo de carrera promedio. ¿Cómo puedo conseguir eso?

ACTUALIZACIÓN 2, CASI:

Yo pensé que encontré lo que estaba buscando en Gelman y Hill (2007) , página 273. Necesidad de utilizar el armpaquete.

library("arm")

Para el pescado 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 

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. Solo estoy teniendo en cuenta la incertidumbre general del modelo. En una situación en la que tengo, digamos, 5 razas observadas para Fish K y 1000 razas observadas para Fish L, creo que la incertidumbre asociada con mi predicción para Fish K debería ser mucho mayor que la incertidumbre asociada con mi predicción para Fish L.

Examinaré más a fondo a Gelman y Hill 2007. Siento que podría terminar teniendo que cambiar a BUGS (o Stan).

ACTUALIZAR EL 3er:

Quizás estoy conceptualizando mal las cosas. El uso de la predictInterval()función dada por Jared Knowles en una respuesta a continuación da intervalos que no son exactamente lo que 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 agregado dos peces nuevos. Fish K, para quien hemos observado 995 razas, y Fish L, para quien hemos observado 5 razas. Hemos observado 100 carreras para Fish AJ. Me quedo igual lmer()que antes. Mirando dotplot()desde el latticepaquete:

FishID Estimates

Por defecto, dotplot()reordena los efectos aleatorios por su estimación puntual. La estimación para Fish L está en la línea superior y tiene un intervalo de confianza muy amplio. Fish K está en la tercera línea y tiene un intervalo de confianza muy estrecho. Esto tiene sentido para mí. Tenemos muchos datos sobre Fish K, pero no muchos sobre Fish L, por lo que tenemos más confianza en nuestras estimaciones sobre la verdadera velocidad de natación de Fish K. Ahora, creo que esto conduciría a un intervalo de predicción estrecho para Fish K, y un intervalo de predicción amplio para Fish L cuando se usa 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()

Intervalo de predicción para peces

Todos esos intervalos de predicción parecen ser idénticos en ancho. ¿Por qué nuestra predicción para Fish K no es más estrecha que las demás? ¿Por qué nuestra predicción para Fish L no es más amplia que otras?

hossibley
fuente
1
predictIntervalincluye el error / incertidumbre para los términos de efectos fijos y aleatorios. En dotplotusted está viendo solamente la incertidumbre debida a la parte aleatoria de la predicción, esencialmente, la incertidumbre en torno a la estimación de las intersecciones específicas de pescado. Si su modelo tiene mucha incertidumbre en el parámetro fijo fishWty este parámetro impulsa la mayor parte del valor predicho, entonces la incertidumbre en torno a cualquier intercepción de peces específica es trivial y no verá una gran diferencia en el ancho de los intervalos. Deberíamos dejar esto más claro en los predictIntervalresultados.
jknowles

Respuestas:

19

Esta pregunta y el excelente intercambio fueron el impulso para crear la predictIntervalfunción en el merToolspaquete. bootMeres el camino a seguir, pero para algunos problemas no es factible computacionalmente generar ajustes de arranque de todo el modelo (en casos donde el modelo es grande).

En esos casos, predictIntervalestá diseñado para usar las arm::simfunciones para generar distribuciones de parámetros en el modelo y luego usar esas distribuciones para generar valores simulados de la respuesta dada la newdataproporcionada por el usuario. Es fácil de usar, todo lo que necesita hacer es:

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

Puede especificar una gran cantidad de otros valores para predictIntervalincluir la configuración del intervalo para los intervalos de predicción, elegir si se debe informar la media o mediana de la distribución y elegir si se incluye o no la varianza residual del modelo.

No es un intervalo de predicción completo porque la variabilidad de los thetaparámetros en el lmerobjeto no está incluida, pero todas las otras variaciones se capturan a través de este método, dando una aproximación bastante decente.

jknowles
fuente
3
¡Esto se ve increíble! Leyendo la viñeta ahora. ¡Gracias!
Hossibley
Los intervalos de predicción no son exactamente lo que esperaba. Ver actualización 3 arriba.
Hossibley
¿ predictInterval()No le gustan los efectos aleatorios anidados? Por ejemplo, usando el msleepconjunto de datos del ggplot2paquete: mod <- lmer(sleep_total ~ bodywt + (1|vore/order), data=msleep); predInt <- predictInterval(merMod=mod, newdata=msleep) Devuelve un error:Error in '[.data.frame'(newdata, , j) : undefined columns selected
hossibley
Apuesto a que no le gustan los efectos anidados. No creo que hayamos probado nada en nuestro conjunto de pruebas. Presentaré un problema en GitHub para investigar eso. También recomendaría probar devtools::install_github("jknowles/merTools")primero la versión de desarrollo de GitHub .
jknowles
2
Como actualización, la última versión de desarrollo de merTools permite efectos anidados. Será empujado a CRAN en breve.
jknowles
15

Haga esto bootMergenerando un conjunto de predicciones para cada réplica de bootstrap paramétrica:

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

La salida de bootMerestá en un objeto no terriblemente transparente "boot", pero podemos obtener las predicciones en bruto del $tcomponente.

¿Cuánto tiempo gana Fish E a Fish D?

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

Los tiempos del pez E están en la columna 5, los tiempos del pez D están en la columna 4, por lo que solo necesitamos saber la proporción en que la columna 5 es menor que la columna 4:

mean(predMat[,5]<predMat[,4])  ## 0.57
Ben Bolker
fuente
Estoy obteniendo algunos resultados inesperados. Si configuro fishWt = 1 para todos los peces en newDat, esperaría que el tiempo medio / medio para Fish A sea ~ 10.1, Fish B ~ 10.2, ..., Fish J ~ 11.0 (ya que su tiempo en los datos de entrenamiento es definido como:) sec <- 9 + rep(1:10, rep(100,10))/10 + oz + rnorm(1000)/10. Cuando uso predict(), los tiempos de predicción para los peces A, E y J son 10.09, 10.49 y 10.99, como se esperaba. Sin embargo, los tiempos medios para el método bootMer que usted describe son: 10.52, 10.59 y 10.50. Hubiera esperado más acuerdo?
Hossibley
Usar use.u=TRUEcomo en: bb <- bootMer(lme1,nsim=200,FUN=predFun,seed=101,use.u=TRUE)parece darme lo que quiero. ¡Gracias!
Hossibley
OK, esto se pone un poco complicado. Tienes que mirar el use.uargumento para bootMer. La pregunta es, cuando dice "incertidumbre en el efecto fijo y el efecto aleatorio", ¿qué quiere decir con 'efecto aleatorio'? ¿Te refieres a la incertidumbre en la variación de los efectos aleatorios o en los modos condicionales (es decir, los efectos específicos de los peces)? Puedes usar use.u=TRUE, pero no creo que necesariamente haga lo que quieras ...
Ben Bolker
Si uso use.u=TRUE, entonces los "valores de u [permanecer] fijados en sus valores estimados". Interpreto esto como significado, cualquiera que sea nuestra estimación puntual de efecto aleatorio para el pez A, se toma como la verdad honesta de Dios, por así decirlo. bootMerasume que no hay error en nuestra estimación del punto RE. Si lo uso use.u=FALSE, ¿ bootMertiene en cuenta las estimaciones de puntos RE? Parece que los bootMerresultados cuando se usan use.u=FALSEson equivalentes (o asintóticamente equivalentes) al uso re.form=NAen la predict()declaración. ¿Es eso cierto?
Hossibley
1
Creo que no está implementado ATM, pero puede extraer las variaciones condicionales de los modos condicionales / BLUP a través de c(attr(ranef(lme1,condVar=TRUE)[[1]],"postVar"))(todas son idénticas en este ejemplo) y luego muestrear esos valores.
Ben Bolker