¿Por qué lme y aov devuelven resultados diferentes para ANOVA de medidas repetidas en R?

Question

Estoy tratando de pasar de usar el ez paquete a lme para ANOVA de medidas repetidas (como espero poder usar contrastes personalizados en con lme ).

Siguiendo los consejos de esta entrada del blog Pude configurar el mismo modelo utilizando ambos aov (al igual que ez cuando se solicite) y lme . Sin embargo, mientras que en el ejemplo dado en ese puesto el F -valores coinciden perfectamente entre aov y lme (Lo he comprobado y lo hacen), no es el caso de mis datos. Aunque el F -Los valores son similares, pero no son iguales.

aov devuelve un valor f de 1,3399, lme devuelve 1,36264. Estoy dispuesto a aceptar el aov resultado como el "correcto" ya que esto es también lo que devuelve el SPSS (y esto es lo que cuenta para mi campo/supervisor).

Preguntas:

1. Sería genial si alguien pudiera explicar por qué existe esta diferencia y cómo puedo utilizar lme para proporcionar resultados creíbles. (También estaría dispuesto a utilizar lmer en lugar de lme para este tipo de cosas, si da el resultado "correcto". Sin embargo, no lo he utilizado hasta ahora).

2. Después de resolver este problema, me gustaría realizar un análisis de contraste. Especialmente me interesaría el contraste de la agrupación de los dos primeros niveles del factor (es decir, c("MP", "MT") ) y compararlo con el tercer nivel del factor (es decir "AC" ). Además, la prueba del tercer frente al cuarto nivel de factor (es decir, "AC" frente a "DA" ).

Datos:

tau.base <- structure(list(id = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 
9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 
22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 
14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 
5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 
19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L), .Label = c("A18K", 
"D21C", "F25E", "G25D", "H05M", "H07A", "H08H", "H25C", "H28E", 
"H30D", "J10G", "J22J", "K20U", "M09M", "P20E", "P26G", "P28G", 
"R03C", "U21S", "W08A", "W15V", "W18R"), class = "factor"), factor = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("MP", "MT", "AC", "DA"
), class = "factor"), value = c(0.9648092876, 0.2128662077, 1, 
0.0607615485, 0.9912814024, 3.22e-08, 0.8073856412, 0.1465590332, 
0.9981672618, 1, 1, 1, 0.9794401938, 0.6102546108, 0.428651501, 
1, 0.1710644881, 1, 0.7639763913, 1, 0.5298989196, 1, 1, 0.7162733447, 
0.7871177434, 1, 1, 1, 0.8560509327, 0.3096989662, 1, 8.51e-08, 
0.3278862311, 0.0953598576, 1, 1.38e-08, 1.07e-08, 0.545290432, 
0.1305621416, 2.61e-08, 1, 0.9834051136, 0.8044114935, 0.7938839461, 
0.9910112678, 2.58e-08, 0.5762677121, 0.4750002288, 1e-08, 0.8584252623, 
1, 1, 0.6020385797, 8.51e-08, 0.7964935271, 0.2238374288, 0.263377904, 
1, 1.07e-08, 0.3160751898, 5.8e-08, 0.3460325565, 0.6842217296, 
1.01e-08, 0.9438301877, 0.5578367224, 2.18e-08, 1, 0.9161424562, 
0.2924856039, 1e-08, 0.8672987992, 0.9266688748, 0.8356425464, 
0.9988463913, 0.2960361777, 0.0285680426, 0.0969063841, 0.6947998266, 
0.0138254805, 1, 0.3494775301, 1, 2.61e-08, 1.52e-08, 0.5393467752, 
1, 0.9069223275)), .Names = c("id", "factor", "value"), class = "data.frame", row.names = c(1L, 
6L, 10L, 13L, 16L, 17L, 18L, 22L, 23L, 24L, 27L, 29L, 31L, 33L, 
42L, 43L, 44L, 45L, 54L, 56L, 58L, 61L, 64L, 69L, 73L, 76L, 79L, 
80L, 81L, 85L, 86L, 87L, 90L, 92L, 94L, 96L, 105L, 106L, 107L, 
108L, 117L, 119L, 121L, 124L, 127L, 132L, 136L, 139L, 142L, 143L, 
144L, 148L, 149L, 150L, 153L, 155L, 157L, 159L, 168L, 169L, 170L, 
171L, 180L, 182L, 184L, 187L, 190L, 195L, 199L, 202L, 205L, 206L, 
207L, 211L, 212L, 213L, 216L, 218L, 220L, 222L, 231L, 232L, 233L, 
234L, 243L, 245L, 247L, 250L))

Y el código:

require(nlme)

summary(aov(value ~ factor+Error(id/factor), data = tau.base))

anova(lme(value ~ factor, data = tau.base, random = ~1|id))

Answer 1

Son diferentes porque el modelo lme está forzando el componente de varianza de id sea mayor que cero. Observando la tabla de anova cruda para todos los términos, vemos que el error cuadrático medio para id es menor que el de los residuos.

> anova(lm1 <- lm(value~ factor+id, data=tau.base))

          Df  Sum Sq Mean Sq F value Pr(>F)
factor     3  0.6484 0.21614  1.3399 0.2694
id        21  3.1609 0.15052  0.9331 0.5526
Residuals 63 10.1628 0.16131   

Cuando calculamos los componentes de la varianza, esto significa que la varianza debida al id será negativa. Mi memoria de los cuadrados medios esperados es inestable, pero el cálculo es algo así como

(0.15052-0.16131)/3 = -0.003597.

Esto suena a impar pero puede ocurrir. Lo que significa es que los promedios de cada ID están más cerca unos de otros de lo que se esperaría dada la cantidad de variación residual en el modelo.

En cambio, el uso de lme obliga a que esta varianza sea mayor que cero.

> summary(lme1 <- lme(value ~ factor, data = tau.base, random = ~1|id))
...
Random effects:
 Formula: ~1 | id
        (Intercept)  Residual
StdDev: 3.09076e-05 0.3982667

Esto informa de las desviaciones estándar, elevando al cuadrado para obtener la varianza se obtiene 9.553e-10 para la varianza del id y 0.1586164 para la varianza residual.

Ahora bien, debes saber que el uso de aov para medidas repetidas sólo es apropiado si se cree que la correlación entre todos los pares de medidas repetidas es idéntica; esto se llama simetría compuesta. (Técnicamente, se requiere esfericidad, pero esto es suficiente por ahora). Una razón para utilizar lme en aov es que puede manejar diferentes tipos de estructuras de correlación.

En este conjunto de datos en particular, la estimación de esta correlación es negativa; esto ayuda a explicar cómo el error cuadrático medio de id fue menor que el error cuadrático residual. Una correlación negativa significa que si la primera medición de un individuo estuvo por debajo de la media, en promedio, su segunda estaría por encima de la media, haciendo que las medias totales de los individuos sean menos variables de lo que esperaríamos si hubiera una correlación cero o una correlación positiva.

Utilizando lme con un efecto aleatorio equivale a ajustar un modelo de simetría compuesta en el que se obliga a que esa correlación no sea negativa; podemos ajustar un modelo en el que se permite que la correlación sea negativa utilizando gls :

> anova(gls1 <- gls(value ~ factor, correlation=corCompSymm(form=~1|id),
                    data=tau.base))
Denom. DF: 84 
            numDF   F-value p-value
(Intercept)     1 199.55223  <.0001
factor          3   1.33985   0.267

Esta tabla de ANOVA coincide con la tabla del aov ajuste y de la lm encaja.

Bien, ¿y qué? Bueno, si usted cree que la variación de id y la correlación entre las observaciones debe ser no negativa, el lme es en realidad más apropiado que el ajuste utilizando aov o lm ya que su estimación de la varianza residual es ligeramente mejor. Sin embargo, si cree que la correlación entre las observaciones podría ser negativa, aov o lm o gls es mejor.

También puede interesarle explorar la estructura de correlación más a fondo; para ver una estructura de correlación general, haría algo como

gls2 <- gls(value ~ factor, correlation=corSymm(form=~unclass(factor)|id),
data=tau.base)

Aquí sólo limito la salida a la estructura de correlación. Los valores 1 a 4 representan los cuatro niveles de factor vemos que el factor 1 y el factor 4 tienen una correlación negativa bastante fuerte:

> summary(gls2)
...
Correlation Structure: General
 Formula: ~unclass(factor) | id 
 Parameter estimate(s):
 Correlation: 
  1      2      3     
2  0.049              
3 -0.127  0.208       
4 -0.400  0.146 -0.024

Una forma de elegir entre estos modelos es con una prueba de razón de verosimilitud; esto muestra que el modelo de efectos aleatorios y el modelo de estructura de correlación general no son estadísticamente significativos; cuando esto ocurre, se suele preferir el modelo más sencillo.

> anova(lme1, gls2)
     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
lme1     1  6 108.0794 122.6643 -48.03972                        
gls2     2 11 111.9787 138.7177 -44.98936 1 vs 2 6.100725  0.2965

Answer 2

aov() se ajusta al modelo a través de lm() utilizando los mínimos cuadrados, lme se ajusta por máxima verosimilitud. Esa diferencia en la forma de estimar los parámetros del modelo lineal probablemente explique la (muy pequeña) diferencia en sus valores f.

En la práctica, (por ejemplo, para la comprobación de hipótesis) estas estimaciones son las mismas, por lo que no veo cómo una podría considerarse "más creíble" que la otra. Proceden de diferentes paradigmas de ajuste de modelos.

Para los contrastes, es necesario establecer una matriz de contraste para sus factores. Venebles y Ripley muestran cómo hacerlo en las páginas 143, 146 y 293-294 de la 4ª edición.