Escribir la ecuación matemática para un modelo de efectos mixtos multinivel

La pregunta de CV

Estoy tratando de dar (a) representaciones matemáticas detalladas y concisas de un modelo de efectos mixtos. Estoy usando el lme4paquete en R. ¿Cuál es la representación matemática correcta para mi modelo?

---

Los datos, la pregunta científica y el código R

Mi conjunto de datos consta de especies en diferentes regiones. Estoy probando si la prevalencia de una especie cambia en el tiempo que lleva a una extinción (las extinciones no son necesariamente permanentes; puede recolonizarse), o después de una colonización.

lmer(prevalence ~ time + time:type + (1 + time + type:time | reg) + (1 + time + type:time | reg:spp))

• La prevalencia es la proporción de estratos ocupados por una especie en una región año
• El tiempo es una variable continua que indica el tiempo de extinción o colonización; siempre es positivo
• El tipo es una variable categórica con dos niveles. Estos dos niveles son "-" y "+". Cuando type es -, es una colonización (nivel predeterminado). Cuando el tipo es +, es una extinción.
• Reg es una variable categórica con nueve niveles, que indica la región
• Spp es una variable categórica; El número de niveles varía según las regiones, y varía entre 48 niveles y 144 niveles.

En palabras: la variable de respuesta es la prevalencia (proporción de estratos ocupados). Los efectos fijos incluyeron 1) e intercepción, 2) tiempo desde el evento y 3) la interacción entre el tiempo hasta el evento y el tipo de evento (colonización o extinción). Cada uno de estos 3 efectos fijos varía aleatoriamente entre las regiones. Dentro de una región, cada uno de los efectos varía aleatoriamente entre especies.

Estoy tratando de descubrir cómo escribir la ecuación matemática para el modelo. Creo que entiendo lo que está sucediendo en el código R (aunque estoy seguro de que tengo algunos vacíos de conocimiento, y espero que escribir la expresión matemática formal mejore mi comprensión).

He buscado bastante en la web y en estos foros. Encontré toneladas de información útil, para estar seguro (y tal vez vincularé a algunos de estos en una edición de esta pregunta). Sin embargo, no pude encontrar esa "Piedra de Rosetta" del código R traducida a las matemáticas (estoy más cómodo con el código) que realmente me ayudaría a confirmar que tengo estas ecuaciones correctas. De hecho, sé que ya hay algunas lagunas, pero llegaremos a eso.

---

Mi intento

La forma básica de un modelo de efectos mixtos, en notación matricial es (a mi entender):

$$Y = X \beta + Z \gamma + \epsilon$$

$$X = \begin{bmatrix} 1 & \Delta t & \Delta t_{+} \\ \vdots & \vdots & \vdots \\ 1 & \Delta t_n & \Delta t_{+,n} \end{bmatrix}$$ Z= [ 1 I ( r 1 ) Δ t I ( r 1 ) Δ t + I ( r 1 ) … 1 I ( r 9 ) Δ t I ( r 9 ) Δ t + I ( r 9 ) ⋮ ⋮ ⋮ ⋱ $$\beta^{'} = \begin{bmatrix} \beta_0 & \beta_1 & \beta_2 \end{bmatrix}$$ $$Z = \begin{bmatrix} 1 I(r_1) & \Delta t I(r_1) & \Delta t_{+} I(r_1) & \dots & 1 I(r_9) & \Delta t I(r_9) & \Delta t_{+} I(r_9) \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 1 I(r_{1,n}) & \Delta t_n I(r_{1,n}) & \Delta t_{+,n} I(r_{1,n}) & \dots & 1 I(r_{9,n}) & \Delta t I(r_{9,n}) & \Delta t_{+,n} I(r_{9,n}) \\ \end{bmatrix}$$ ϵ∼N(0,Σ $$\gamma^{'} = \begin{bmatrix} \gamma_{0,1} & \gamma_{1,1} &\gamma_{2,1} & \dots & \gamma_{0,9} & \gamma_{1,9} &\gamma_{2,9} \end{bmatrix}$$ $$\epsilon \sim \mathcal{N}(0,\Sigma)$$

• Δ t Δ t $X$ es la matriz de diseño para los efectos fijos, es el tiempo después de la colonización ( ) y es el tiempo después de la extinción ( )$\Delta t$time$\Delta t_{+}$time:type
• $Z$ es la matriz de diseño para los efectos aleatorios (nivel 1?), I () es la función del indicador que da 1 si la muestra pertenece a la región designada y 0 de lo contrario, r se indexa para indicar una de las nueve regiones.
• $\beta$ y contienen parámetros$\gamma$
• $\epsilon$ es errores; No estoy completamente seguro de cómo explicar , aunque me doy cuenta de que una de estas matrices de varianza / covarianza expresará covarianzas entre pendientes e intersecciones, por ejemplo$\Sigma$

Asumiendo que las cosas hasta ahora son ~ correctas, eso significa que soy bueno en el nivel superior. Sin embargo, explicar la variación específica de la especie en los parámetros, que está anidado dentro de cada región, me dejó perplejo aún más.

Pero tomé una grieta en algo que tal vez tenga sentido ...

Cada uno de los parámetros en se deriva de una combinación lineal de predictores y parámetros específicos de la especie dentro de una región. Para cada región, hay 3 filas de, correspondientes a las 3 variables predictoras. Cada se puede expresar individualmente como$\gamma$$\gamma$

• $\gamma_{p,r} = U_{p,r} b_{p,r} + \eta_{p,r}$
  • donde es una matriz de diseño específica para la región y el predictor , es una matriz de parámetros 1 por S para la región (riqueza en la región = , por ejemplo, 48 o 144), y es una matriz de términos de error$U_{p,r}$$r$$p$$b_{p,r}$$S$$\eta_{p,r}$

Específicamente, para una región determinada, cada uno de los sería:$\gamma_{p,r}$

$$\gamma_{0,r} = U_{0,r} b_{0,r} + \eta_{0,r}$$ $$\gamma_{0,r} = \begin{bmatrix} 1 I(s_1) \dots 1 I(s_S) \end{bmatrix} + \begin{bmatrix} b_{0,1}\\ \vdots \\ b_{0,S} \end{bmatrix} + \eta_{0,r}$$ $$\gamma_{1,r} = U_{1,r} b_{1,r} + \eta_{1,r}$$ $$\gamma_{1,r} = \begin{bmatrix} \Delta t I(s_1) \dots \Delta t I(s_S) \end{bmatrix} + \begin{bmatrix} b_{1,1}\\ \vdots \\ b_{1,S} \end{bmatrix} + \eta_{1,r}$$ $$\gamma_{2,r} = U_{2,r} b_{2,r} + \eta_{2,r}$$ $$\gamma_{2,r} = \begin{bmatrix} \Delta t_+ I(s_1) \dots \Delta t_+ I(s_S) \end{bmatrix} + \begin{bmatrix} b_{2,1}\\ \vdots \\ b_{2,S} \end{bmatrix} + \eta_{2,r}$$

Eso se repetiría para cada región. Luego, , como . Aunque, quizás en lugar de , hay otra letra, como , que se usa comúnmente.$\eta \sim \mathcal{N}(0,\Sigma_{\eta})$$\epsilon$$\Sigma$$G$

---

Editar: otras preguntas y respuestas que fueron algo útiles

• Este Q / A fue agradable, pero no escribió las cosas en forma de matriz completa

Si entendí el código correctamente, ¿por qué no simplemente escribir algo como

$$y_{i} = \Big(\alpha + \nu_{j[i]}^{(\alpha)} + \eta_{k[i]}^{(\alpha)}\Big) + \Big(\beta + \nu_{j[i]}^{(\beta)} + \eta_{k[i]}^{(\beta)}\Big)T_{i} + \Big(\delta + \nu_{j[i]}^{(\delta)} + \eta_{k[i]}^{(\delta)}\Big)(T_{i} * Z_{i}) + \epsilon_i$$ con o, si la primera ecuación es demasiado larga, algo como y $$\begin{aligned} \Big[\nu_{j}^{(\alpha)}, \nu_j^{(\beta)}, \nu_j^{(\delta)}\Big] &\sim \text{Multi-Normal}(\mathbf 0, \boldsymbol \Sigma_\nu) \\ \Big[\eta_{j}^{(\alpha)}, \eta_j^{(\beta)}, \eta_j^{(\delta)}\Big] &\sim \text{Multi-Normal}(\mathbf 0, \boldsymbol \Sigma_\eta)\\ \epsilon_i & \sim \text{Normal}(0, \sigma_\epsilon) \end{aligned}$$ $$y_{i} = \alpha_{j[i],k[i]} + \beta_{j[i],k[i]}T_{i} + \delta_{j[i],k[i]}(T_i * Z_i) + \epsilon_i$$ $$\begin{aligned} \alpha_{j[i],k[i]} &= \alpha + \nu_{j}^{(\alpha)} + \eta_{k}^{(\alpha)} \\ \beta_{j[i],k[i]}&=\beta + \nu_{j}^{(\beta)} + \eta_{k}^{(\beta)}\\ \delta_{j[i],k[i]}&=\delta + \nu_{j}^{(\delta)} + \eta_{k}^{(\delta)}\\ \end{aligned}$$ con la misma estructura de covarianza ¿como anteriormente? Muestra la estructura anidada de los datos, así como qué coeficientes varían según los niveles.