Notación para modelos multinivel

Question

La fórmula que hay que especificar para entrenar un modelo multinivel (utilizando lmer de lme4 R biblioteca) siempre me atrapa. He leído innumerables libros de texto y tutoriales, pero nunca lo he entendido bien.

He aquí un ejemplo de este tutorial que me gustaría ver formulado en una ecuación. Estamos intentando modelar la frecuencia de la voz en función del sexo (las mujeres tienen un tono de voz más agudo que los hombres en general) y de la actitud de la persona (si ha contestado de manera educada o informal) en diferentes escenarios. Además, como puede verse en el subject columna, cada persona fue sometida a mediciones varias veces.

> head(politeness, n=20)
   subject gender scenario attitude frequency
1       F1      F        1      pol     213.3
2       F1      F        1      inf     204.5
3       F1      F        2      pol     285.1
4       F1      F        2      inf     259.7
5       F1      F        3      pol     203.9
6       F1      F        3      inf     286.9
7       F1      F        4      pol     250.8
8       F1      F        4      inf     276.8
9       F1      F        5      pol     231.9
10      F1      F        5      inf     252.4
11      F1      F        6      pol     181.2
12      F1      F        6      inf     230.7
13      F1      F        7      inf     216.5
14      F1      F        7      pol     154.8
15      F3      F        1      pol     229.7
16      F3      F        1      inf     237.3
17      F3      F        2      pol     236.8
18      F3      F        2      inf     251.0
19      F3      F        3      pol     267.0
20      F3      F        3      inf     266.0

subject , gender y attitude son factores (con informal y female considerados como niveles de base para attitude y gender en las ecuaciones siguientes). Ahora, una idea es entrenar un modelo con diferentes interceptos para cada subject y scenario :

politeness.model=lmer(frequency ~ attitude + gender + 
 (1|subject) + (1|scenario), data=politeness)

Si mi comprensión de la notación es correcta, esto corresponde a:

$y_i=a^1_{j[i]}+a^2_{k[i]}+\beta\cdot$ attitude$_{\text{pol}_i} + \gamma\cdot$ gender$_{\text{male}_i}$

donde $i$ indica $i^{th}$ punto de datos, $j[i]$ denota el nivel de grupo para subject y $k[i]$ denota el nivel de grupo para scenario para $i^{th}$ punto de datos. attitude$_\text{pol}$ y gender$_\text{male}$ son indicadores binarios.

Para introducir pendientes aleatorias para la actitud, podemos escribir:

politeness.model = lmer(frequency ~ attitude + gender + 
 (1+attitude|subject) + (1+attitude|scenario), data=politeness)

De nuevo, si mi comprensión es clara, esto corresponde a:

$y_i = a^1_{j[i]} + a^2_{k[i]} + (\beta^1_{j[i]} + \beta^2_{k[i]})\cdot$ attitude$_{\text{pol}_i} + \gamma\cdot$ gender$_{\text{male}_i}$

Ahora bien, ¿qué ecuación R corresponde?

politeness.null = lmer(frequency ~ gender +
 (1+attitude|subject) +  (1+attitude|scenario), data=politeness)

Answer 1

Yo escribiría

~ attitude + gender + (1|subject) + (1|scenario)

como

$$ y_i \sim \beta_0 + \beta_1 \cdot I(\textrm{attitude}=\textrm{pol}) + \beta_2 I(\textrm{gender}=\textrm{male}) + b_{1,j[i]} + b_{2,k[i]} + \epsilon_i \\ b_1 \sim N(0,\sigma^2_1) \\ b_2 \sim N(0,\sigma^2_2) \\ \epsilon \sim N(0,\sigma^2_r) $$ donde $\beta$ indica un coeficiente de efecto fijo, $b$ indica una variable aleatoria, $I$ es una función indicadora (esto es básicamente lo mismo que lo que has dicho antes, sólo que con una notación ligeramente diferente).

~ attitude + gender + (1+attitude|subject) + (1+attitude|scenario)

añade variación entre sujetos en respuesta a attitude y scenario (podríamos escribir de forma equivalente la parte de efectos aleatorios como (attitude|subject) + (attitude|scenario) es decir, dejando implícito el intercepto; esto es cuestión de gustos). Ahora

$$ y_i \sim \beta_0 + \beta_1 \cdot I(\textrm{attitude}=\textrm{pol}) + \beta_2 I(\textrm{gender}=\textrm{male}) + \\ b_{1,j[i]} + b_{3,j[i]} I(\textrm{attitude}=\textrm{pol}) + b_{2,k[i]} + b_{4,k[i]} I(\textrm{attitude}=\textrm{pol}) + \epsilon_i \\ \{b_1,b_3\} \sim \textrm{MVN}({\mathbf 0},\Sigma_1) \\ \{b_2,b_4\} \sim \textrm{MVN}({\mathbf 0},\Sigma_2) \\ \epsilon \sim N(0,\sigma^2_r) $$ donde $\Sigma_1$ y $\Sigma_2$ son matrices de varianza-covarianza no estructuradas, es decir, son simétricas y (semi)definidas positivas, pero no tienen ninguna otra restricción: $$ \Sigma_1 = \left( \begin{array}{cc} \sigma^2_1 & \sigma_{13} \\ \sigma_{13} & \sigma^2_3 \end{array} \right) $$ y análogamente para $\Sigma_2$ .

Puede resultar instructivo agrupar los términos de la siguiente manera: $$ y_i \sim (\beta_0 + b_{1,j[i]} + b_{2,k[i]}) + \\ ( \beta_1 + b_{3,j[i]} + b_{4,k[i]}) \cdot I(\textrm{attitude}=\textrm{pol}) + \beta_2 I(\textrm{gender}=\textrm{male}) + \epsilon_i $$ para que pueda ver qué efectos aleatorios afectan al intercepto y cuáles afectan a la respuesta a la actitud.

Ahora bien, si se deja de lado el efecto fijo attitude (es decir, establecer $\beta_1=0$ o suelte el attitude de la fórmula) se puede ver (sin reescribirlo todo) que, como se supone que los efectos aleatorios tienen media cero, estaremos suponiendo que la respuesta media a la actitud en todos los sujetos y escenarios será exactamente cero, aunque sigue habiendo variaciones entre sujetos y escenarios. No diré esto nunca tiene sentido desde un punto de vista estadístico, pero rara vez lo tiene. Hay discusiones sobre este tema en la lista de correo r-sig-mixed-models@r-project.org de vez en cuando ... (o puede ser discutido en StackExchange en algún lugar - si no, sería una buena pregunta de seguimiento SE ...)