¿Cómo puedo predecir mejor con (g)lmer con valores faltantes?

Question

Supongamos que estoy construyendo un modelo mixto en R, y quiero usar ese modelo para predecir nuevos datos para los cuales podría no conocer el valor de todas las características. O en algunos casos, puede que no sea tanto que no conozca el valor de todas las características, sino más bien que quiero una predicción para un valor "neutral" de esas características en algún sentido significativo. ¿Existe una forma "correcta" de hacer esto?

Aquí hay un ejemplo tonto.

install.packages("lme4")
require(lme4)
mtcars2 = mtcars
mtcars2$cyl = paste0(mtcars2$cyl, " cilindros")
# Quiero asegurarme de no tratar 'cyl' numéricamente. Esta
# probablemente no sea la mejor aproximación, pero debería funcionar.

modelo = lmer(mpg ~ (1|cyl) + hp, mtcars2)

Con mi modelo tonto ya construido, puedo usarlo para hacer una predicción sobre una nueva fila de datos.

ejemploTest = data.frame(c("4 cilindros"), c(100))
names(ejemploTest) = c("cyl", "hp")

prediccionTest = predict(modelo, ejemploTest)
print(prediccionTest)

Esto produce un resultado numérico, como se esperaba por supuesto. Pero supongamos que quisiera saber qué predice el modelo si solo conociera el valor de hp y no conociera el valor de cyl.

ejemploTest$cyl = c("Nuevo valor no visto antes")
prediccionTest = predict(modelo, ejemploTest, allow.new.levels=TRUE)
print(prediccionTest)

Puedo agregar allow.new.levels=TRUE, y esto al menos me permite obtener un número (21.66101). No estoy convencido de que sea el número correcto, pero volveremos a eso en un momento.

Me parece natural que debería poder aplicar el mismo tipo de lógica si en su lugar tuviera un valor legítimo para cyl pero hp estuviera ausente. La forma natural para mí sería poner NA en lugar del valor numérico para hp de la siguiente manera.

ejemploTest$cyl = c("4 cilindros") # restaurando el valor original
ejemploTest$hp = c(NA)
prediccionTest = predict(modelo, ejemploTest, allow.new.levels=TRUE)
print(prediccionTest)

Este enfoque no produce un número. No estoy seguro si la diferencia entre los dos ejemplos radica más en la diferencia entre efectos fijos y aleatorios o en uno siendo numérico y el otro no.

Entonces, ¿hay alguna forma directa de manejar esto en R? Creo (pero no estoy seguro) que lo que podría querer es que se impute la media en mi caso numérico de hp. Eso no parece tan simple en mi caso de cyl, ya que no puedo tomar la media de algunas cadenas. Parece que lo que allow.new.levels hace efectivamente en mi caso es ignorar la parte del modelo para esa característica por completo como si produjera un valor de 0, y simplemente calcula el intercepto del efecto fijo más 100 veces la pendiente del efecto fijo. ¿Hay formas estándar de manejar este tipo de situaciones?

Answer 1

Un par de puntos:

• Comencemos con una definición general del modelo mixto lineal, es decir, sea $y_i$ el resultado para el i-ésimo sujeto ($i = 1, \ldots, n$), $X_i$ la matriz de diseño para los efectos fijos $\beta$, y $Z_i$ la matriz de diseño para los efectos aleatorios $b_i$, entonces tenemos: $$\left\{ \begin{array}{l} y_i = X_i \beta + Z_i b_i + \varepsilon_i,\\\\ b_i \sim \mathcal N(0, D), \quad \varepsilon_i \sim \mathcal N(0, \sigma^2 \mbox{I}), \end{array} \right.$$ donde $\varepsilon_i$ son los términos de error. Cuando ajustas este modelo obtienes estimaciones para los efectos fijos $\hat \beta$. También puedes obtener estimaciones de los efectos aleatorios $\hat b_i$ que se definen como los modos de la distribución condicional: $$p(b_i \mid y_i; \hat \theta) = \frac{p(y_i \mid b_i; \hat \theta) \, p(b_i; \hat \theta)}{p(y_i; \hat \theta)},$$ donde $\theta = (\beta, \mbox{vech}(D), \sigma^2)$ denota todos los parámetros del modelo.
• De este modelo puedes obtener dos tipos de predicciones, es decir, predicciones poblacionales definidas como $$\hat y_i^{pop} = X_i \hat \beta,$$ y predicciones específicas por sujeto, definidas como $$\hat y_i^{subj} = X_i \hat \beta + Z_i \hat b_i.$$
• Supongamos ahora que deseas calcular la predicción para un nuevo sujeto $j$. Si no tienes datos de resultado para este sujeto, entonces la única predicción que puedes hacer es la predicción poblacional, es decir, $$\hat y_j^{pop} = X_j \hat \beta,$$ que será el promedio para su grupo. Sin embargo, si también tienes algunos datos $y_j^{new}$ para este sujeto, entonces podrías estimar su efecto aleatorio $\hat b_j$ como el modo de $[b_j \mid y_j^{new}; \hat \theta]$, y luego usar esto en la ecuación anterior de la predicción específica por sujeto, es decir, $$\hat y_j^{subj} = X_j \hat \beta + Z_j \hat b_j.$$

Puedes encontrar información adicional al respecto en el Capítulo 3, Sección 3.4 de mis apuntes del curso de Mediciones Repetidas.