¿Cómo elegir la estructura de efectos aleatorios y fijos en los modelos lineales mixtos?

Question

Considere los siguientes datos de un diseño de dos vías dentro de los sujetos:

df <- "http://personality-project.org/r/datasets/R.appendix4.data"
df <- read.table(df,header=T)
head(df)

Observation Subject Task Valence Recall
1           1     Jim Free     Neg      8
2           2     Jim Free     Neu      9
3           3     Jim Free     Pos      5
4           4     Jim Cued     Neg      7
5           5     Jim Cued     Neu      9
6           6     Jim Cued     Pos     10

Me gustaría analizarlo utilizando modelos lineales mixtos. Teniendo en cuenta todos los posibles efectos fijos y aleatorios, hay múltiples modelos posibles:

# different fixed effects with random-intercept
a0 <- lmer(Recall~1 + (1|Subject), REML=F,df)
a1 <- lmer(Recall~Task + (1|Subject), REML=F,df)
a2 <- lmer(Recall~Valence + (1|Subject), REML=F,df)
a3 <- lmer(Recall~Task+Valence + (1|Subject), REML=F,df)
a4 <- lmer(Recall~Task*Valence + (1|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope
b0 <- lmer(Recall~1 + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b1 <- lmer(Recall~Task + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b2 <- lmer(Recall~Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b3 <- lmer(Recall~Task+Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b4 <- lmer(Recall~Task*Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope including variance-covariance matrix
c0 <- lmer(Recall~1 + (1 + Valence + Task|Subject), REML=F,df)
c1 <- lmer(Recall~Task + (1 + Valence + Task|Subject), REML=F,df)
c2 <- lmer(Recall~Valence + (1 + Valence + Task|Subject), REML=F,df)
c3 <- lmer(Recall~Task+Valence + (1 + Valence + Task|Subject), REML=F,df)
c4 <- lmer(Recall~Task*Valence + (1 + Valence + Task|Subject), REML=F,df)

1. ¿Cuál es la forma recomendada para seleccionar el modelo que mejor se ajusta en este contexto? Cuando se utilizan pruebas de relación de logaritmo-verosimilitud, ¿cuál es el procedimiento recomendado? ¿Generar modelos hacia arriba (del modelo nulo al más complejo) o hacia abajo (del modelo más complejo al nulo)? ¿Inclusión o exclusión por pasos? ¿O se recomienda poner todos los modelos en una prueba de razón de logaritmos y seleccionar el modelo con el valor p más bajo? ¿Cómo se comparan los modelos que no están anidados?

2. ¿Es recomendable encontrar primero la estructura de efectos fijos adecuada y luego la estructura de efectos aleatorios adecuada o al revés (he encontrado referencias para ambas opciones...)?

3. ¿Cuál es la forma recomendada de comunicar los resultados? Informar del valor p de la prueba de la razón de verosimilitud logarítmica que compara el modelo mixto completo (con el efecto en cuestión) con el modelo reducido (sin el efecto en cuestión). ¿O es mejor utilizar la prueba de relación de logaritmo-verosimilitud para encontrar el modelo que mejor se ajuste y luego utilizar lmerTest para informar de los valores p de los efectos en el modelo que mejor se ajuste?

Answer 1

No estoy seguro de que haya una respuesta canónica a esto, pero lo intentaré.

¿Cuál es la forma recomendada de seleccionar el modelo que mejor se ajusta en este contexto? Cuando se utilizan pruebas de relación de logaritmo-verosimilitud, ¿cuál es el procedimiento recomendado? ¿Generar modelos hacia arriba (del modelo nulo al más complejo) o hacia abajo (del modelo más complejo al nulo)? ¿Inclusión o exclusión por pasos? ¿O se recomienda poner todos los modelos en una prueba de relación de logaritmo-verosimilitud y seleccionar el modelo con el valor p más bajo? ¿Cómo se comparan los modelos que no están anidados?

Depende de cuáles sean tus objetivos.

• En general, debe ser muy , muy cuidado con la selección del modelo (véase, por ejemplo, el esta respuesta o este puesto , o simplemente buscar en Google "Harrell stepwise"...).
• Si le interesa que sus valores p sean significativos (es decir, si está haciendo pruebas de hipótesis confirmatorias), debe no hacer la selección del modelo. Sin embargo, : no me queda tan claro si los procedimientos de selección de modelos son tan malos si está haciendo la selección del modelo en partes no focales del mismo Por ejemplo, hacer la selección del modelo en los efectos aleatorios si su interés principal es la inferencia en los efectos fijos.
• No existe tal cosa como "poner todos los modelos en una prueba de razón de verosimilitud" - la prueba de razón de verosimilitud es un procedimiento por pares. Si Si quisiera hacer una selección de modelos (por ejemplo) sobre los efectos aleatorios, probablemente recomendaría un enfoque "todo a la vez" utilizando criterios de información como en este ejemplo -- que al menos evita algunos de los problemas de los enfoques por pasos (pero no de la selección de modelos en general).
• Barr et al. 2013 "Mantener el máximo" Revista de Memoria y Lenguaje (doi:10.1016/j.jml.2012.11.001) recomendaría utilizar el modelo maximalista (únicamente).
• Shravan Vasishth no está de acuerdo argumentando que estos modelos no van a tener suficiente potencia y, por lo tanto, van a ser problemáticos a menos que el conjunto de datos sea muy grande (y la relación señal-ruido sea alta)
• Otro planteamiento razonablemente defendible es el de encajar un gran pero razonable modelo y luego, si el ajuste es singular, eliminar términos hasta que no sea más
• Con algunas salvedades (enumeradas en el PREGUNTAS FRECUENTES DE LA GLMM ), puede utilizar los criterios de información para comparar modelos no anidados con diferentes efectos aleatorios (aunque Brian Ripley no está de acuerdo: véase al final de la página 6 aquí )

¿Es recomendable encontrar primero la estructura de efectos fijos adecuada y luego la estructura de efectos aleatorios adecuada o al revés (he encontrado referencias para ambas opciones...)?

No creo que nadie lo sepa. Véase la respuesta anterior sobre la selección de modelos en general. Si pudieras definir tus objetivos con suficiente claridad (cosa que poca gente hace), la pregunta podría ser respondida. Si tiene referencias para ambas opciones, sería útil editar su pregunta para incluirlas... (Por si sirve de algo, este ejemplo (ya citado anteriormente) utiliza los criterios de información para seleccionar la parte de efectos aleatorios, y luego evita la selección en la parte de efectos fijos del modelo.

¿Cuál es la forma recomendada de comunicar los resultados? Informar del valor p de la prueba de la razón de verosimilitud logarítmica comparando el modelo mixto completo (con el efecto en cuestión) con el modelo reducido (sin el efecto en cuestión). ¿O es mejor utilizar la prueba de relación de logaritmo-verosimilitud para encontrar el modelo que mejor se ajuste y luego utilizar lmerTest para informar de los valores p de los efectos en el modelo que mejor se ajuste?

Esta es (por desgracia) otra pregunta difícil. Si se informa de la efectos marginales según informa lmerTest Hay que preocuparse por la marginalidad (por ejemplo, si las estimaciones de los efectos principales de A y B tienen sentido cuando hay una A -por- B interacción en el modelo); esto es una enorme lata de gusanos, pero es Algo así como mitigado si se utiliza contrasts="sum" según la recomendación de afex::mixed() . Los diseños equilibrados también ayudan un poco. Si realmente quieres empapelar todas estas grietas, creo que recomendaría afex::mixed que le da una salida similar a lmerTest pero trata de abordar estas cuestiones.

Answer 2

Actualización de mayo de 2017 : Resulta que mucho de lo que he escrito aquí es un poco mal . A lo largo del post se realizan algunas actualizaciones.

---

Estoy muy de acuerdo con lo que ya ha dicho Ben Bolker (gracias por el saludo a afex::mixed() ), pero permítanme añadir algunas reflexiones más generales y específicas sobre esta cuestión.

Enfoque en los efectos fijos frente a los aleatorios y cómo comunicar los resultados

Para el tipo de investigación experimental que se representa en el conjunto de datos de ejemplo de Jonathan Baron que utiliza la pregunta importante suele ser si un manipulado tiene un efecto global. Por ejemplo, ¿encontramos un efecto principal global o una interacción de Task ? Un punto importante es que en esos conjuntos de datos normalmente todos los factores están bajo completo control experimental y asignados al azar. En consecuencia, el interés suele centrarse en los efectos fijos.
Por el contrario, los componentes de los efectos aleatorios pueden verse como parámetros "molestos" que capturan la varianza sistemática (es decir, las diferencias interindividuales en el tamaño del efecto) que no son necesariamente importantes para la cuestión principal. Desde este punto de vista, la sugerencia de utilizar la estructura de efectos aleatorios máximos, tal y como defienden Barr et al., es algo natural. Es fácil imaginar que una manipulación experimental no afecta a todos los individuos exactamente de la misma manera y queremos controlar esto. Por otra parte, el número de factores o niveles no suele ser demasiado grande, por lo que el peligro de sobreajuste parece comparativamente pequeño.

En consecuencia, seguiría la sugerencia de Barr et al. y especificaría una estructura de efectos aleatorios máximos e informaría de las pruebas de los efectos fijos como mis resultados principales. Para probar los efectos fijos también sugeriría utilizar afex::mixed() ya que informa de las pruebas de efectos o factores (en lugar de las pruebas de parámetros) y calcula esas pruebas de una manera algo sensata (por ejemplo, utiliza la misma estructura de efectos aleatorios para todos los modelos en los que se elimina un solo efecto, utiliza contrastes de suma a cero, ofrece diferentes métodos para calcular p -valores, ...).

¿Qué pasa con los datos del ejemplo?

El problema con los datos del ejemplo que has dado es que para este conjunto de datos la estructura de efectos aleatorios máximos conduce a un modelo sobresaturado ya que sólo hay un punto de datos por celda del diseño:

> with(df, table(Valence, Subject, Task))
, , Task = Cued

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

, , Task = Free

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

En consecuencia, lmer se ahoga en la estructura de efectos aleatorios máximos:

> lmer(Recall~Task*Valence + (Valence*Task|Subject), df)
Error: number of observations (=30) <= number of random effects (=30) for term
(Valence * Task | Subject); the random-effects parameters and the residual variance
(or scale parameter) are probably unidentifiable

Desgraciadamente, hasta donde yo sé, no hay una forma acordada de abordar este problema. Pero permítanme esbozar y discutir algunas:

1. Una primera solución podría ser eliminar la pendiente aleatoria más alta y probar los efectos para este modelo:

   require(afex)
   mixed(Recall~Task*Valence + (Valence+Task|Subject), df)
           Effect    F ndf  ddf F.scaling p.value
   1         Task 6.56   1 4.00      1.00     .06
   2      Valence 0.80   2 3.00      0.75     .53
   3 Task:Valence 0.42   2 8.00      1.00     .67

   Sin embargo, esta solución es un poco ad hoc y no está demasiado motivada.

   Actualización de mayo de 2017: Este es el enfoque que estoy apoyando actualmente. Véase esta entrada del blog y el borrador del capítulo del que soy coautor En la sección "Estructuras de efectos aleatorios para diseños ANOVA tradicionales".

2. Una solución alternativa (y que podría considerarse como defendida por la discusión de Barr et al.) podría ser eliminar siempre las pendientes aleatorias para el efecto más pequeño. Sin embargo, esto tiene dos problemas: (1) Qué estructura de efectos aleatorios utilizamos para averiguar cuál es el efecto más pequeño y (2) R es reacio a eliminar un efecto de orden inferior, como un efecto principal, si están presentes efectos de orden superior, como una interacción de este efecto (véase aquí ). En consecuencia, habría que configurar esta estructura de efectos aleatorios a mano y pasar la matriz del modelo así construida a la llamada lmer.

3. Una tercera solución podría ser utilizar una parametrización alternativa de la parte de efectos aleatorios, a saber, la que corresponde al modelo RM-ANOVA para estos datos. Desgraciadamente (?), lmer no permite "varianzas negativas", por lo que esta parametrización no se corresponde exactamente con el RM-ANOVA para todos los conjuntos de datos , ver discusión aquí y en otros lugares (por ejemplo aquí y aquí ). El "lmer-ANOVA" para estos datos sería:

   > mixed(Recall~Task*Valence + (1|Subject) + (1|Task:Subject) + (1|Valence:Subject), df)
           Effect    F ndf  ddf F.scaling p.value
   1         Task 7.35   1 4.00      1.00     .05
   2      Valence 1.46   2 8.00      1.00     .29
   3 Task:Valence 0.29   2 8.00      1.00     .76

Teniendo en cuenta todos estos problemas yo simplemente no usaría lmer para ajustar conjuntos de datos para los que sólo hay un punto de datos por celda del diseño, a menos que se disponga de una solución más consensuada para el problema de la estructura de efectos aleatorios máximos.

4. En su lugar, me gustaría También se podría seguir utilizando el ANOVA clásico. Utilizando uno de los wrappers para car::Anova() en afex los resultados serían:

   > aov4(Recall~Task*Valence + (Valence*Task|Subject), df)
           Effect         df  MSE      F  ges   p
   1      Valence 1.44, 5.75 4.67   1.46  .02 .29
   2         Task       1, 4 4.08 7.35 +  .07 .05
   3 Valence:Task 1.63, 6.52 2.96   0.29 .003 .71

   Tenga en cuenta que afex ahora también permite devolver el modelo ajustado con aov que se puede pasar a lsmeans para las pruebas post-hoc (pero para la prueba de los efectos los reportados por car::Anova siguen siendo más razonables):

   > require(lsmeans)
   > m <- aov4(Recall~Task*Valence + (Valence*Task|Subject), df, return = "aov")
   > lsmeans(m, ~Task+Valence)
    Task Valence lsmean       SE   df lower.CL upper.CL
    Cued Neg       11.8 1.852026 5.52  7.17157 16.42843
    Free Neg       10.2 1.852026 5.52  5.57157 14.82843
    Cued Neu       13.0 1.852026 5.52  8.37157 17.62843
    Free Neu       11.2 1.852026 5.52  6.57157 15.82843
    Cued Pos       13.6 1.852026 5.52  8.97157 18.22843
    Free Pos       11.0 1.852026 5.52  6.37157 15.62843
   
   Confidence level used: 0.95