Por que lme e aov retornam resultados diferentes para medidas repetidas ANOVA em R?

Estou tentando deixar de usar o ezpacote lmepara medidas repetidas ANOVA (como espero poder usar contrastes personalizados lme).

Seguindo o conselho desta postagem no blog, consegui configurar o mesmo modelo usando os dois aov(como faz ez, quando solicitado) e lme. No entanto, enquanto no exemplo dado nesse post os valores F concordam perfeitamente entre aove lme(eu verifiquei e eles concordam ), esse não é o caso dos meus dados. Embora os valores- F sejam semelhantes, eles não são os mesmos.

aovretorna um valor f de 1,33399, lmeretorna 1,36264. Estou disposto a aceitar o aovresultado como o "correto", pois também é isso que o SPSS retorna (e é isso que conta para o meu campo / supervisor).

Questões:

1. Seria ótimo se alguém pudesse explicar por que essa diferença existe e como eu posso usar lmepara fornecer resultados confiáveis. (Eu também estaria disposto a usar lmer, em vez de lmepara este tipo de coisa, se ele dá o resultado "correta". No entanto, eu não usei até agora.)

2. Depois de resolver esse problema, gostaria de executar uma análise de contraste. Especialmente, eu estaria interessado no contraste de agrupar os dois primeiros níveis de fator (ie c("MP", "MT")) e compará-lo com o terceiro nível de fator (ie "AC"). Além disso, testando o terceiro versus o quarto nível de fator (ou seja, "AC"versus "DA").

Dados:

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))

E o código:

require(nlme)

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

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

Eles são diferentes porque o modelo lme está forçando o componente de variação ida ser maior que zero. Observando a tabela anova bruta para todos os termos, vemos que o erro quadrático médio para id é menor do que para os resíduos.

> 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   

Quando computamos os componentes da variação, isso significa que a variação devido ao id será negativa. Minha memória de quadrados médios esperados é instável, mas o cálculo é algo como

(0.15052-0.16131)/3 = -0.003597.

Isso parece estranho, mas pode acontecer. O que isso significa é que as médias para cada ID estão mais próximas umas das outras do que você esperaria uma da outra, dada a quantidade de variação residual no modelo.

Por outro lado, o uso de lme força essa variação a ser maior que zero.

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

Isso relata desvios padrão, ao quadrado para obter os rendimentos da variação 9.553e-10para a variação de id e 0.1586164para a variação residual.

Agora, você deve saber que o uso aovpara medidas repetidas só é apropriado se você acredita que a correlação entre todos os pares de medidas repetidas é idêntica; isso é chamado de simetria composta. (Tecnicamente, a esfericidade é necessária, mas isso é suficiente por enquanto.) Uma razão para usar lmemais aové que ela pode lidar com diferentes tipos de estruturas de correlação.

Nesse conjunto de dados em particular, a estimativa para essa correlação é negativa; isso ajuda a explicar como o erro quadrático médio para id foi menor que o erro quadrático residual. Uma correlação negativa significa que se a primeira medida de um indivíduo estivesse abaixo da média, em média, a segunda estaria acima da média, tornando as médias totais para os indivíduos menos variáveis ​​do que seria de esperar se houvesse uma correlação zero ou uma correlação positiva.

Usar lmecom efeito aleatório é equivalente a ajustar um modelo de simetria composta, onde essa correlação é forçada a não ser negativa; podemos ajustar um modelo em que a correlação pode ser negativa usando 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 tabela ANOVA concorda com a tabela do aovajuste e do lmajuste.

Tá, e daí? Bem, se você acredita que a variação de ide a correlação entre as observações devem ser não-negativas, o lmeajuste é realmente mais apropriado do que o ajuste usando aovou lmporque sua estimativa da variação residual é um pouco melhor. No entanto, se você acredita que a correlação entre as observações poderia ser negativo, aovou lmou glsé melhor.

Você também pode estar interessado em explorar ainda mais a estrutura de correlação; para olhar para uma estrutura geral de correlação, você faria algo como

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

Aqui, limito apenas a saída à estrutura de correlação. Os valores 1 a 4 representam os quatro níveis de factor; vemos que o fator 1 e o fator 4 têm uma correlação negativa bastante forte:

> 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

Uma maneira de escolher entre esses modelos é com um teste de razão de verossimilhança; isso mostra que o modelo de efeitos aleatórios e o modelo geral de estrutura de correlação não são estatisticamente significativamente diferentes; quando isso acontece, o modelo mais simples é geralmente 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

aov()se encaixa no modelo lm()usando o mínimo de quadrados, lmese encaixa com a máxima probabilidade. Essa diferença na forma como os parâmetros do modelo linear são estimados provavelmente explica a (muito pequena) diferença em seus valores f.

Na prática, (por exemplo, para testes de hipóteses), essas estimativas são as mesmas, portanto não vejo como uma pode ser considerada 'mais credível' que a outra. Eles vêm de diferentes paradigmas de ajuste de modelo.

Para contrastes, você precisa configurar uma matriz de contraste para seus fatores. Venebles e Ripley mostram como fazer isso nas páginas 143, p.146 e p.293-294 da 4ª edição.