Uso de lmer para el modelo lineal de medidas repetidas de efectos mixtos

Question

EDIT 2: Originalmente pensé que tenía que ejecutar un ANOVA de dos factores con medidas repetidas en un factor, pero ahora creo que un modelo lineal de efectos mixtos funcionará mejor para mis datos. Creo que casi sé lo que tiene que pasar, pero todavía estoy confundido por algunos puntos.

Los experimentos que necesito analizar son los siguientes:

• Los sujetos fueron asignados a uno de varios grupos de tratamiento
• Las mediciones de cada sujeto se realizaron en varios días
• Así que:
  • El sujeto está anidado dentro del tratamiento
  • El tratamiento se cruza con el día

(cada sujeto se asigna a un solo tratamiento, y las mediciones se realizan en cada sujeto en cada día)

Mi conjunto de datos contiene la siguiente información:

• Sujeto = factor de bloqueo (factor aleatorio)
• Día = factor intra-sujeto o de medidas repetidas (factor fijo)
• Tratamiento = factor entre sujetos (factor fijo)
• Obs = variable medida (dependiente)

ACTUALIZACIÓN Vale, fui a hablar con un estadístico, pero es un usuario de SAS. Él piensa que el modelo debe ser:

Tratamiento + Día + Sujeto(Tratamiento) + Día*Sujeto(Tratamiento)

Obviamente su notación es diferente de la sintaxis de R, pero se supone que este modelo da cuenta de:

• Tratamiento (fijo)
• Día (fijo)
• la interacción Tratamiento*Día
• Sujeto anidado dentro del tratamiento (aleatorio)
• Día cruzado con "Sujeto dentro del tratamiento" (aleatorio)

Entonces, ¿es esta la sintaxis correcta a utilizar?

m4 <- lmer(Obs~Treatment*Day + (1+Treatment/Subject) + (1+Day*Treatment/Subject), mydata)

Me preocupa especialmente si la parte de "Día cruzado con el sujeto dentro del tratamiento" es correcta. ¿Alguien que esté familiarizado con SAS, o que esté seguro de entender lo que sucede en su modelo, puede comentar si mi triste intento de sintaxis en R coincide?

Aquí están mis intentos anteriores de construir un modelo y escribir la sintaxis (discutido en las respuestas y comentarios):

m1 <- lmer(Obs ~ Treatment * Day + (1 | Subject), mydata)

¿Cómo puedo tratar el hecho de que el sujeto esté anidado dentro del tratamiento? ¿Cómo se m1 difieren de:

m2 <- lmer(Obs ~ Treatment * Day + (Treatment|Subject), mydata)
m3 <- lmer(Obs ~ Treatment * Day + (Treatment:Subject), mydata)

y son m2 y m3 equivalente (y si no, por qué)?

Además, ¿tengo que utilizar nlme en lugar de lme4 si quiero especificar la estructura de correlación (como correlation = corAR1 )? Según Medidas repetidas Para un análisis de medidas repetidas con medidas repetidas en un factor, la estructura de covarianza (la naturaleza de las correlaciones entre las medidas del mismo sujeto) es importante.

Cuando intenté hacer un ANOVA de medidas repetidas, decidí utilizar un SS de tipo II; ¿sigue siendo relevante, y si es así, cómo lo especifico?

Este es un ejemplo de cómo se ven los datos:

mydata <- data.frame(
  Subject  = c(13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 
               34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 
               19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
               40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 
               29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65), 
  Day       = c(rep(c("Day1", "Day3", "Day6"), each=28)), 
  Treatment = c(rep(c("B", "A", "C", "B", "C", "A", "A", "B", "A", "C", "B", "C", 
                      "A", "A", "B", "A", "C", "B", "C", "A", "A"), each = 4)), 
  Obs       = c(6.472687, 7.017110, 6.200715, 6.613928, 6.829968, 7.387583, 7.367293, 
                8.018853, 7.527408, 6.746739, 7.296910, 6.983360, 6.816621, 6.571689, 
                5.911261, 6.954988, 7.624122, 7.669865, 7.676225, 7.263593, 7.704737, 
                7.328716, 7.295610, 5.964180, 6.880814, 6.926342, 6.926342, 7.562293, 
                6.677607, 7.023526, 6.441864, 7.020875, 7.478931, 7.495336, 7.427709, 
                7.633020, 7.382091, 7.359731, 7.285889, 7.496863, 6.632403, 6.171196, 
                6.306012, 7.253833, 7.594852, 6.915225, 7.220147, 7.298227, 7.573612, 
                7.366550, 7.560513, 7.289078, 7.287802, 7.155336, 7.394452, 7.465383, 
                6.976048, 7.222966, 6.584153, 7.013223, 7.569905, 7.459185, 7.504068, 
                7.801867, 7.598728, 7.475841, 7.511873, 7.518384, 6.618589, 5.854754, 
                6.125749, 6.962720, 7.540600, 7.379861, 7.344189, 7.362815, 7.805802, 
                7.764172, 7.789844, 7.616437, NA, NA, NA, NA))

Answer 1

Creo que su enfoque es correcto. Modelo m1 especifica un intercepto separado para cada sujeto. Modelo m2 añade una pendiente distinta para cada sujeto. Su pendiente es a través de los días ya que los sujetos sólo participan en un grupo de tratamiento. Si escribe el modelo m2 como sigue es más obvio que se modele un intercepto y una pendiente por separado para cada sujeto

m2 <- lmer(Obs ~ Treatment * Day + (1+Day|Subject), mydata)

Esto equivale a:

m2 <- lmer(Obs ~ Treatment + Day + Treatment:Day + (1+Day|Subject), mydata)

Es decir, los efectos principales del tratamiento, el día y la interacción entre ambos.

Creo que no hay que preocuparse por el anidamiento siempre que no se repitan los ID de los sujetos dentro de los grupos de tratamiento. El modelo correcto depende realmente de su pregunta de investigación. ¿Hay razones para creer que las pendientes de los sujetos varían además del efecto del tratamiento? Podría ejecutar ambos modelos y compararlos con anova(m1,m2) para ver si los datos apoyan alguna de ellas.

No estoy seguro de lo que quieres expresar con el modelo m3 ? La sintaxis de anidamiento utiliza un / Por ejemplo (1|group/subgroup) .

No creo que haya que preocuparse por la autocorrelación con un número tan pequeño de puntos temporales.

Answer 2

No me siento lo suficientemente cómodo como para comentar tu problema de errores autocorrelacionados (ni sobre las diferentes implementaciones en lme4 vs. nlme), pero puedo hablar del resto.

Su modelo m1 es un modelo de intercepción aleatoria, en el que se ha incluido la interacción a nivel cruzado entre el Tratamiento y el Día (se permite que el efecto del Día varíe entre los grupos de Tratamiento). Para permitir que el cambio a lo largo del tiempo difiera entre los participantes (es decir, para modelar explícitamente las diferencias individuales en el cambio a lo largo del tiempo), también debe permitir que el efecto del Día sea al azar . Para ello, hay que especificar:

m2 <- lmer(Obs ~ Day + Treatment + Day:Treatment + (Day | Subject), mydata)

En este modelo:

• El intercepto si la puntuación prevista para la categoría de referencia del tratamiento en el Día=0
• El coeficiente de Día es el cambio previsto a lo largo del tiempo por cada aumento de 1 unidad de días para la categoría de referencia del tratamiento
• Los coeficientes de los dos códigos ficticios para los grupos de tratamiento (creados automáticamente por R) son la diferencia prevista entre cada grupo de tratamiento restante y la categoría de referencia en el Día=0
• Los coeficientes de los dos términos de interacción son la diferencia del efecto del tiempo (Día) sobre las puntuaciones previstas entre la categoría de referencia y los demás grupos de tratamiento.

Tanto los interceptos como el efecto del Día en la puntuación son aleatorios (se permite que cada sujeto tenga una puntuación predicha diferente en el Día=0 y un cambio lineal diferente a lo largo del tiempo). La covarianza entre los interceptos y las pendientes también se modela (se permite que covarié).

Como puede ver, la interpretación de los coeficientes de las dos variables ficticias es condicional al Día=0. Le indicarán si la puntuación prevista en el Día=0 para la categoría de referencia es significativamente diferente de los dos grupos de tratamiento restantes. Por lo tanto, el lugar en el que se decida centrar la variable Día es importante. Si se centra en el Día 1, los coeficientes le indican si la puntuación prevista para la categoría de referencia en el día 1 es significativamente diferente de la puntuación prevista de los dos grupos restantes. De esta manera, se puede ver si hay diferencias preexistentes entre los grupos . Si se centra en el Día 3, los coeficientes le indican si la puntuación prevista para la categoría de referencia al tercer día es significativamente diferente de la puntuación prevista de los dos grupos restantes. De esta manera, se puede ver si hay diferencias entre los grupos al final de la intervención .

Por último, hay que tener en cuenta que los sujetos son no anidado dentro de Tratamiento. Sus tres tratamientos no son niveles aleatorios de una población de niveles a los que quiere generalizar sus resultados; más bien, como ha mencionado, sus niveles son fijos y quiere generalizar sus resultados sólo a estos niveles. (Por no mencionar que no debería utilizar la modelización multinivel si sólo tiene 3 unidades de nivel superior; véase Maas y Hox, 2005). En cambio, el tratamiento es un predictor de nivel 2, es decir, un predictor que toma un único valor a través de los Días (unidades de nivel 1) para cada sujeto. Por lo tanto, se incluye simplemente como un predictor en su modelo.

Referencia:
Maas, C. J. M., & Hox, J. J. (2005). Tamaños de muestra suficientes para la modelización multinivel. Metodología: Revista Europea de Métodos de Investigación para las Ciencias Sociales y del Comportamiento , 1 , 86-92.