24 votos

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

Estoy tratando de pasar de un uso de la ez paquete a lme ANOVA de medidas repetidas (como tengo la esperanza de ser capaz de utilizar personalizado contrasta con lme).

Siguiendo los consejos de este blog que yo era capaz de establecer el mismo modelo con el tanto aov (como ez, cuando se solicita) y lme. Sin embargo, mientras que en el ejemplo dado en que después de la F-valores perfectamente de acuerdo entre aov y lme (lo he comprobado, y lo hacen), este no es el caso de mis datos. Aunque el F-valores son similares, no son el mismo.

aov devuelve un valor f de 1.3399, lme devuelve 1.36264. Estoy dispuesto a aceptar la aov resultado como la "correcta" uno como esto es también lo SPSS devuelve (y esto es lo que cuenta para mi campo/supervisor).

Preguntas:

  1. Sería genial si alguien pudiera explicar por qué esta diferencia existe y cómo puedo usar lme a proporcionar resultados fiables. (Yo también estaría dispuesto a utilizar lmer en lugar de lme para este tipo de cosas, si se le da a la "correcta". Sin embargo, no lo he utilizado hasta ahora.)

  2. Después de resolver este problema me gustaría ejecutar un análisis de contraste. Especialmente yo estaría interesado en el contraste de la agrupación de los dos primeros niveles del factor (es decir, c("MP", "MT")) y comparar con el tercer nivel de factor (es decir, "AC"). Además, las pruebas de la tercera versus el cuarto nivel de factor (es decir, "AC" versus "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))

21voto

Raptrex Puntos 115

Son diferentes debido a la lme modelo está forzando a los componentes de varianza de id ser mayor que cero. Mirando el raw tabla anova para todos los términos, vemos que el error cuadrático medio para la identificación 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 se calculan las componentes de varianza, esto significa que la varianza debida a la identificación será negativa. Mi memoria de espera media de los cuadrados de la memoria es débil, pero el cálculo es algo así como

(0.15052-0.16131)/3 = -0.003597.

Esto suena raro, pero puede suceder. Lo que significa es que los promedios para cada id están más cerca unos de otros de lo que puedes esperar de cada uno de los otros, dada la cantidad de variación residual en el modelo.

En contraste, el uso de lme de las fuerzas de esta variación debe ser 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

Este informe de desviaciones estándar, ajustando para obtener la varianza de los rendimientos 9.553e-10 para la identificación y la varianza 0.1586164 de la varianza residual.

Ahora, usted debe saber que el uso de aov para medidas repetidas, sólo es apropiado si usted cree que la correlación entre todos los pares de medidas repetidas es idéntico; esto se llama el compuesto de la simetría. (Técnicamente, esfericidad es necesario, pero esto es suficiente por ahora.) Una razón para usar lme sobre aov es que puede manejar diferentes tipos de correlación de las estructuras.

En este conjunto particular de datos, la estimación de esta correlación es negativa; esto ayuda a explicar cómo el error cuadrático medio para la identificación de la era menor que el residual de cuadrados de error. Una correlación negativa significa que si un individuo de la primera medición por debajo de la media, en promedio, su segunda sería por encima de la media, haciendo que el total de los promedios de los individuos menos variable que nos esperaría si había una correlación cero o una correlación positiva.

El uso de lme con un efecto aleatorio es equivalente a la colocación de un compuesto de simetría modelo en el que la correlación es forzado a ser no negativo; se puede ajustar a un modelo donde la correlación es permitido ser negativos en 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 ANOVA está de acuerdo con la tabla de la aov de ajuste y de la lm de ajuste.

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

Usted también puede estar interesado en la exploración de la estructura de las correlaciones más; mirar una correlación general de la estructura, tendría que hacer algo como

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

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

> 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 manera de elegir entre estos modelos es con una prueba de razón de verosimilitud; esto muestra que el modelo de efectos aleatorios y la correlación general de la estructura del modelo no son diferentes estadísticamente significativas; cuando esto ocurre, el modelo más simple es generalmente preferido.

> 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

20voto

Julien Chastang Puntos 161

aov() ajusta el modelo a través de la lm() el uso de mínimos cuadrados, lme se ajusta a través de máxima verosimilitud. Que diferencia en cómo los parámetros del modelo lineal que se estima probable que las cuentas de la (muy pequeño) la diferencia en su f-valores.

En la práctica, (por ejemplo, para la prueba de hipótesis) estas estimaciones son iguales, así que no veo cómo uno puede ser considerado "más creíble" que la otra. Vienen de diferentes ajuste del modelo de paradigmas.

Por los contrastes, que usted necesita para establecer un contraste de la matriz de sus factores. Venebles y Ripley a mostrar cómo hacer esto en la página 143, p.146 y p.293-294 de la 4ª edición.

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X