¿Por qué es difícil incorporar la incertidumbre en los efectos aleatorios al hacer predicciones de modelos mixtos?

Hay varios hilos en R-sig-ME sobre cómo obtener intervalos de confianza para las predicciones usando lme4y nlmeen R. Por ejemplo, aquí y aquí en 2010, incluidos algunos comentarios de Dougals Bates, uno de los autores de ambos paquetes. Dudo en citarlo literalmente, por temor a que sean sacados de contexto, pero de todos modos, un comentario que hace es

"Está combinando parámetros y variables aleatorias en sus predicciones y no estoy seguro de lo que significaría evaluar la variabilidad de esas predicciones. Un bayesiano puede tener sentido, pero no puedo entenderlo". " https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q1/003447.html

Sé que el paquete Bayesian glmm MCMCglmmpuede producir intervalos creíbles para las predicciones.

Últimamente, la versión de desarrollo de lme4on github ha recibido un predictmétodo, pero va acompañada del siguiente comentario:

"@note 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 \ code {\ link {bootMer}} para esta tarea". https://github.com/lme4/lme4/blob/master/R/predict.R

Entonces, ¿por qué es difícil incorporar la incertidumbre en los efectos aleatorios al hacer predicciones de modelos mixtos en un entorno frecuentista?

No estoy seguro sobre el comentario del método de predicción, pero un problema principal está relacionado con la generación de medidas de varianza fácilmente interpretables, no con las medidas de varianza per se. Bates no está comentando en la primera cita si puede hacerlo, sino lo que significa.

Tome un modelo simple de varios niveles de un diseño de medidas repetidas de dos niveles. Digamos que tiene los siguientes datos donde cada línea es un asunto:

En lmerel modelo se podría expresar como:

y ~ x + (1|subject)

Está prediciendo el valor y de x como un efecto fijo (la diferencia entre A y B); y la intercepción de un efecto aleatorio **. Mire cuidadosamente el gráfico y observe que si bien existe una variabilidad en el efecto x para cada sujeto (pendiente de cada línea) es relativamente pequeño en comparación con la variabilidad entre sujetos (la altura de cada línea).

El modelo analiza estos dos conjuntos de variabilidad y cada uno es significativo. Puede usar los efectos aleatorios para predecir las alturas de las líneas y puede usar los efectos fijos de x para predecir las pendientes. Incluso podría usar los dos combinados para trabajar nuestros valores y individuales. Pero lo que no puede hacer es realmente decir algo significativo con respecto a su modelo cuando combina la variabilidad de las pendientes y las alturas de las líneas. Debe hablar sobre la variabilidad de sus pendientes y alturas de líneas por separado. Esa es una característica del modelo, no una responsabilidad.

Tendrá una variabilidad del efecto de x que es relativamente fácil de estimar. Se podría decir algo sobre un intervalo de confianza en torno a eso. Pero tenga en cuenta que este intervalo de confianza tendrá una pequeña relación con la predicción de cualquier valor y particular porque el valor y está influenciado por una combinación de efecto y varianza del sujeto que es diferente de la variabilidad del efecto solo.

Cuando Bates escribe cosas como las que usted ha citado, imagino que a menudo piensa en diseños de niveles múltiples mucho más complejos que esto ni siquiera aborda. Pero incluso si considera este simple ejemplo, se pregunta qué tipo de significado real se puede extraer de la combinación de todas las medidas de varianza.

** Ignoré el efecto fijo de la intercepción por simplicidad y solo lo traté como un efecto aleatorio. Podría extraer conclusiones similares de un modelo aún más simple con una intercepción aleatoria y fija solamente, pero creo que sería más difícil de transmitir. En ese caso, una vez más, el efecto fijo y el efecto aleatorio se analizan por una razón y significan cosas diferentes y volver a unir su variabilidad para los valores pronosticados hace que esa variabilidad tenga poco sentido con respecto al modelo.

Durante mucho tiempo me he preguntado sobre la creencia aparentemente común de que existe una diferencia fundamental en los efectos fijos y aleatorios para los modelos de efectos mixtos (generalmente no lineales). Esta creencia es, por ejemplo, declarada por Bates en la siguiente respuesta

https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q1/003447.html

Bates afirma claramente que cree que existe una diferencia fundamental entre los efectos fijos y aleatorios para que no puedan combinarse. Creo que está equivocado y espero convencer a algunos lectores de un punto de vista alternativo. Tomo un enfoque frecuentista, así que lo que quiero hacer es definir una noción de probabilidad de perfil para una función de los efectos fijos y aleatorios. Para motivar la discusión, supongamos que tenemos un modelo de dos parámetros con los parámetros x y u (nada sobre efectos aleatorios hasta ahora). Sea la función de probabilidad donde suprimimos cualquier referencia a los datos. Supongamos que sea ​​cualquier función (agradable) de x y u. La probabilidad de perfil para la función viene dada porg ( x , u ) P g ( t ) $L(x,u)$$g(x,u)$$P_g(t)$$g$

$$P_g(t)=\max_{x,u} \{L(x,u)\ |\ g(x,u)=t \} \eqno(1)$$

Creo que nadie discutiría con esto. Ahora supongamos que tenemos una distribución de probabilidad previa para u. Entonces afirmaría que la probabilidad de perfil para todavía tiene sentido, pero deberíamos modificar (1) incluyendo el anterior.$p(u)$$g$

$$P_g(t)=\max_{x,u} \{L(x,u)p(u)\ |\ g(x,u)=t \} \eqno(2)$$ Tenga en cuenta que dado que es un parámetro con a prior es exactamente lo mismo que se conoce como efecto aleatorio. Entonces, ¿por qué mucha gente piensa que los parámetros de efectos aleatorios son de alguna manera diferentes? Creo que la diferencia proviene de la práctica habitual de estimación de parámetros para ellos. Lo que hace que los efectos aleatorios sean `` diferentes '' es que hay muchos de ellos en muchos modelos. Como resultado, para obtener estimaciones útiles para los efectos fijos (u otros parámetros), es necesario tratar los efectos aleatorios de una manera diferente. Lo que hacemos es integrarlos fuera del modelo. En el modelo anterior formaríamos la probabilidad donde Ahora el$u$$F(x)$ $$F(x) = \int L(x,u)p(u)du$$ $u$se fueron. Entonces, si todo lo que tenemos es , parece que no tiene sentido hablar sobre la probabilidad de perfil para alguna función .$F(x)$$g(x,u)$

Por lo tanto, para obtener información sobre la función no debemos integrar el parámetro . Pero lo que sucede en el caso donde hay muchos parámetros de efectos aleatorios. Luego afirmo que deberíamos integrarnos sobre `` la mayoría '', pero no todos, en un sentido que voy a precisar. Para motivar la construcción, deje que haya efectos aleatorios . Considere el caso especial donde la función solo depende de , y de hecho es la función más simple imaginable, . Integre sobre los efectos aleatorios para obtener $g(x,u)$$u$$n$$u=(u_1,u_2,...,u_{n-1},u_n)$$g(x,u)$$u_n$$g(x,u)=u_n$$u_1,u_2,...,u_{n-1}$

$$F(x,u_n) = \int L(x,u_1,...,u_n)p(u_1,...,u_n))du_1du_2...du_{n-1}\eqno(4)$$ como antes podemos formar el perfil de probabilidad Cómo generalizar para que tenga sentido para una función arbitraria . en cuenta que la definición de en es la misma que Para ver esta nota que para el caso simple , es lo mismo que $$P_g(t)=\max_{x,u_n} \{F(x,u_n) | u_n=t \} \eqno(3)$$ $(3)$$g(x,u)$$F(x,u_n)$$(4)$ $$F(x,s) = \lim_{\epsilon\rightarrow 0}{1\over\epsilon} \int_{\{(x,u_n) | s-\epsilon/2<g(x,u_n)<s+\epsilon/2\}} L(x,u_1,...,u_n)p(u_1,...,u_n))du_1du_2...du_n\eqno(5)$$ $g(x,u)=u_n$$(5)$ $$F(x,s)=\lim_{\epsilon\rightarrow 0}{1\over\epsilon} \int_{\{(x,u_n) | s-\epsilon/2<u_n<s+\epsilon/2\}} F(x,u_n)du_n\eqno(6)$$

Para una función general formamos la función definida por y calculamos la probabilidad de perfil $g(x,u)$$F(x,s)$$(5)$

$$P_g(s)=\max_{x,u} \{F(x,s) | g(x,u)=s \} \eqno(3)$$

La probabilidad de este perfil es un concepto bien definido y se sostiene por sí mismo. Sin embargo, para ser útil en la práctica, uno debe poder calcular su valor, al menos aproximadamente. Creo que para muchos modelos la función se puede aproximar suficientemente bien usando una variante de la aproximación de Laplace. Defina por Sea H el hessian del logaritmo de la función con respecto a los parámetros y .$F(x,s)$$\hat x(s),\hat u(s)$

$$\hat x(s),\hat u(s)= \max_{x,u} \{L(x,u)p(u)\ |\ g(x,u)=s\}$$ $-L(x,u)p(u)$$x$$u$

Los conjuntos de niveles de son submanifolds dimensionales de un espacio dimensional donde hay efectos fijos y efectos aleatorios. Necesitamos integrar una forma sobre este múltiple donde todo está linealizado en Esto implica un poco de geometría diferencial elemental. Suponga que Al volver a parametrizar podemos suponer que y . Entonces considera el mapa $g$$m+n-1$$n+m$$m$$n$$n$$du_1\wedge du_2\wedge\ldots\wedge du_n$$\hat x(s),\hat u(s)$$g_{x_n}(\hat x(s),\hat u(s))\ne 0$$\hat x(s)=0$$\hat u(s)=0$

$$(x_1,x_2,\ldots,x_{m-1},u_1,u_2,\ldots,u_n) \rightarrow (x_1,x_2,\ldots,x_{m-1}, {-\sum_{i=1}^{m-1}g_{x_i}x_i-\sum_{i=1}^ng_{u_i}u_i\over g_{x_m}}, u_1,u_2,\ldots,u_n)$$ donde se utiliza para denota la derivada parcial de con respecto a evaluada en el punto máximo. Este es un mapa lineal del espacio dimensional en el espacio tangente del conjunto de niveles de . Podemos usarlo para calcular la integral deseada. Primero, el retroceso de las formas 1 son simplemente ellos mismos.$g_{x_i}$$g$$x_i$$m+n-1$$g$$du_i$

El retroceso del hessiano es la forma cuadrática

$$T_{i,j} =H_{i+m,j+m}+{g_{u_i}g_{u_j}\over {g_{x_m}}^2}H_{m,m}\quad \hbox{\rm for} \ 1<=i,j<=n$$

Por lo tanto, la integral se puede calcular (o aproximar) mediante la aproximación de Laplace, que es la fórmula habitual que implica el logaritmo del determinante de , que se calcula mediante la descomposición de Cholesky. El valor de la aproximación de Laplace de la integral es dondeEs el determinante. todavía tenemos que tratar con el ancho del conjunto de niveles de como Para ordenar primero esto tiene el valor donde es el vector de derivadas parciales de $T$

$$L(\hat x(s),\hat u(s))|-T|^{1\over2}$$ $|\cdot|$$g$$\epsilon\rightarrow 0$$\epsilon/\|\nabla g(\hat x(s),\hat u(s))\|$$\nabla g(\hat x(s),\hat u(s)))$$g$ $( g_{x_1}, g_{x_2}, \ldots, g_{x_m}, g_{u_1}, g_{u_2}, \ldots, g_{u_n})$ para que el valor de probabilidad en el conjunto de nivel de se da por Esta es la aproximación correcta para usar para calcular la probabilidad de perfil.$g$ $${L(\hat x(s),\hat u(s))|-T|^{1\over2}\over \|\nabla g(\hat x(s),\hat u(s))\|}$$