Comprobación de supuestos modelos mixtos lmer / lme en R

25

Ejecuté un diseño repetido mediante el cual probé 30 hombres y 30 mujeres en tres tareas diferentes. Quiero entender cómo el comportamiento de hombres y mujeres es diferente y cómo eso depende de la tarea. Utilicé tanto el paquete lmer como el lme4 para investigar esto, sin embargo, estoy atascado tratando de verificar las suposiciones para cualquiera de los métodos. El código que ejecuto es

lm.full <- lmer(behaviour ~ task*sex + (1|ID/task), REML=FALSE, data=dat)
lm.full2 <-lme(behaviour ~ task*sex, random = ~ 1|ID/task, method="ML", data=dat)

Verifiqué si la interacción era el mejor modelo comparándolo con el modelo más simple sin la interacción y ejecutando un anova:

lm.base1 <- lmer(behaviour ~ task+sex+(1|ID/task), REML=FALSE, data=dat)
lm.base2 <- lme(behaviour ~ task+sex, random= ~1|ID/task), method="ML", data=dat)
anova(lm.base1, lm.full)
anova(lm.base2, lm.full2)

P1: ¿Está bien usar estos predictores categóricos en un modelo mixto lineal?
P2: ¿Entiendo correctamente que está bien que la variable de resultado ("comportamiento") no necesite ser distribuida normalmente (por sexo / tareas)?
P3: ¿Cómo puedo verificar la homogeneidad de la varianza? Para un modelo lineal simple que uso plot(LM$fitted.values,rstandard(LM)). ¿Es plot(reside(lm.base1))suficiente el uso ?
P4: Para verificar la normalidad, ¿está utilizando el siguiente código?

hist((resid(lm.base1) - mean(resid(lm.base1))) / sd(resid(lm.base1)), freq = FALSE); curve(dnorm, add = TRUE)
crazjo
fuente
Una cosa que noté también es que la versión de lme4 que estaba usando no era la más reciente y, por lo tanto, la trama simple (myModel.lm) no funcionó, quizás esto sea útil para que otros lectores lo sepan ...
crazjo

Respuestas:

26

P1: Sí, como cualquier modelo de regresión.

P2: Al igual que los modelos lineales generales, su variable de resultado no necesita distribuirse normalmente como una variable univariada. Sin embargo, los modelos LME suponen que los residuos del modelo se distribuyen normalmente. Por lo tanto, una transformación o agregar pesos al modelo sería una forma de encargarse de esto (y verificar con diagramas de diagnóstico, por supuesto).

Q3: plot(myModel.lme)

P4: qqnorm(myModel.lme, ~ranef(., level=2)). Este código le permitirá realizar gráficos QQ para cada nivel de los efectos aleatorios. Los modelos LME suponen que no solo los residuos dentro del grupo están normalmente distribuidos, sino que cada nivel de los efectos aleatorios también. Varíe levelde 0, 1 a 2 para que pueda verificar la rata, la tarea y los residuos dentro del sujeto.

EDITAR: También debería agregar que si bien se supone la normalidad y que la transformación probablemente ayuda a reducir los problemas con errores no normales / efectos aleatorios, no está claro que todos los problemas se resuelvan realmente o que no se introduzca un sesgo. Si sus datos requieren una transformación, tenga cuidado con la estimación de los efectos aleatorios. Aquí hay un documento que aborda esto .

Alce
fuente
Gracias por tu respuesta. Me gustaría compartir mi conjunto de datos y script para el análisis, incluida la salida, para ver si lo que hice es correcto. ¿Es posible en el intercambio de pila? Además, creo que ejecuté el factor aleatorio incorrecto (1 | rata / tarea), ¿no debería ser simplemente (1 | rata)? Probé 60 ratas (30 de cada sexo) en tres tareas.
crazjo
99
Probé el código para Q4 recientemente y recibí un error sobre el objeto de tipo 'S4' no subconjustable. ¿Estaba ese código destinado a modelos que se ajustaban al paquete lme? ¿Qué pasa con lme4?
emudrak
Con respecto al cuarto trimestre, las personas que hacen esos gráficos deben tener en cuenta que la N para cada uno de los gráficos producidos será sustancialmente menor que el total y, por lo tanto, los gráficos serán mucho más variables. No espere que se vean tan consistentemente distribuidos normalmente como en general.
John
14

Parece bastante engañoso sobre los supuestos que rodean a los modelos de varios niveles. No existe una suposición de homogeneidad de varianza en los datos, solo que los residuos deben distribuirse aproximadamente normalmente. Y los predictores categóricos se usan en la regresión todo el tiempo (la función subyacente en R que ejecuta un ANOVA es el comando de regresión lineal).

Para obtener detalles sobre el examen de los supuestos, consulte el libro Pinheiro y Bates (p. 174, sección 4.3.1). Además, si planea usar lme4 (sobre el cual el libro no está escrito) puede replicar sus tramas usando plot con un lmermodelo ( ?plot.merMod).

Para comprobar rápidamente la normalidad, simplemente sería qqnorm(resid(myModel)).

John
fuente
Gracias por tu comentario. ¿Sugieres usar el lmer sobre el método lme4? ¿Y estoy en lo cierto al comprender que la variable de respuesta no necesita distribuirse normalmente? Leeré bien el libro de Pinheiro y Bates.
crazjo
Además, ¿está seguro de que funciona qqnorm (resid (myModel)) en un modelo mixto con múltiples factores?
crazjo
La nueva función lmer tiene más capacidades y un mayor rendimiento. ¿Has probado qqnorm? Siga los consejos al principio del libro sobre cómo leerlo.
John
La trama me había parecido inicialmente extraña, posiblemente porque en realidad no tenía la versión más reciente de lmer. Gracias por notar esto, ahora funciona según sea necesario.
crazjo
12

Con respecto a Q2:

De acuerdo con el libro de Pinheiro y Bates, puede usar el siguiente enfoque:

"La lmefunción permite modelar la heterocedasticidad del grupo dentro del error mediante un weightsargumento. Este tema se tratará en detalle en el § 5.2, pero, por ahora, es suficiente saber que la varIdentestructura de la función de variación permite diferentes variaciones para cada nivel de un factor y puede usarse para ajustarse al modelo heteroscedastic [...] "

Pinheiro y Bates, pág. 177

Si desea verificar las variaciones iguales entre sexusted, puede usar este enfoque:

plot( lm.base2, resid(., type = "p") ~ fitted(.) | sex,
  id = 0.05, adj = -0.3 )

Si las variaciones son diferentes, puede actualizar su modelo de la siguiente manera:

lm.base2u <- update( lm.base2, weights = varIdent(form = ~ 1 | sex) )
summary(lm.base2u)

Además, puede echar un vistazo al robustlmmpaquete que también utiliza un enfoque de pesaje. La tesis doctoral de Koller sobre este concepto está disponible como acceso abierto ("Estimación robusta de modelos lineales mixtos"). El resumen dice:

"Una nueva estimación de la escala, la estimación de la Escala Adaptativa de Diseño, se desarrolla con el objetivo de proporcionar una base sólida para pruebas robustas posteriores. Lo hace igualando la heterocedasticidad natural de los residuos y ajustando la ecuación de estimación robusta para la escala misma . Estas correcciones adaptativas de diseño son cruciales en entornos de muestras pequeñas, donde el número de observaciones podría ser simplemente cinco veces el número de parámetros a estimar o menos ".



No tengo suficientes puntos para comentarios. Sin embargo, veo la necesidad de aclarar algunos aspectos de la respuesta de @John arriba. Pinheiro y Bates afirman en la pág. 174:

Supuesto 1 : los errores dentro del grupo son independientes e idénticamente distribuidos normalmente, con media cero y varianza σ2, y son independientes de los efectos aleatorios.

De hecho, esta afirmación no es clara acerca de las variaciones homogéneas y no soy lo suficientemente profundo en las estadísticas para conocer todas las matemáticas detrás del concepto LME. Sin embargo, en la p. 175, §4.3.1, la sección que trata de la Asunción 1 que escriben:

En esta sección, nos concentramos en métodos para evaluar la suposición de que los errores dentro del grupo se distribuyen normalmente, están centrados en cero y tienen una varianza constante .

Además, en los siguientes ejemplos, las " variaciones constantes " son realmente importantes. Por lo tanto, se puede especular si implican variaciones homogéneas cuando escriben " distribuido de manera idéntica " en la p. 174 sin abordarlo más directamente.

ToJo
fuente
-6

Q1: Sí, ¿por qué no?

P2: Creo que el requisito es que los errores se distribuyan normalmente.

P3: se puede probar con la prueba de Leven, por ejemplo.

usuario12719
fuente