Por que os valores de p mudam de significado ao alterar a ordem das covariáveis ​​no modelo aov?

10

Eu tenho um conjunto de dados de 482 observações.

data=Populationfull

Vou fazer uma análise de associação de genótipo para 3 SNPs. Estou tentando construir um modelo para minha análise e estou usando o aov (y ~ x, data = ...). Para uma característica, tenho vários efeitos fixos e covariáveis ​​que incluí no modelo, como:

Starts <- aov(Starts~Sex+DMRT3+Birthyear+Country+Earnings+Voltsec+Autosec, data=Populationfull)

summary(Starts)
                Df Sum Sq Mean Sq F value   Pr(>F)    
Sex              3  17.90    5.97  42.844  < 2e-16 ***
DMRT3            2   1.14    0.57   4.110    0.017 *  
Birthyear        9   5.59    0.62   4.461 1.26e-05 ***
Country          1  11.28   11.28  81.005  < 2e-16 ***
Earnings         1 109.01  109.01 782.838  < 2e-16 ***
Voltsec          1  12.27   12.27  88.086  < 2e-16 ***
Autosec          1   8.97    8.97  64.443 8.27e-15 ***
Residuals      463  64.48    0.14                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Descobri que se eu alterei a ordem das variáveis ​​no modelo, obtive valores p diferentes, veja abaixo.

Starts2 <- aov(Starts~Voltsec+Autosec+Sex+DMRT3+Birthyear+Country+Earnings, data=Populationfull)

summary(Starts2)

                Df Sum Sq Mean Sq F value   Pr(>F)    
Voltsec   1   2.18    2.18  15.627 8.92e-05 ***
Autosec   1 100.60  100.60 722.443  < 2e-16 ***
Sex              3  10.43    3.48  24.962 5.50e-15 ***
DMRT3            2   0.82    0.41   2.957  0.05294 .  
Birthyear        9   3.25    0.36   2.591  0.00638 ** 
Country          1   2.25    2.25  16.183 6.72e-05 ***
Earnings      1  46.64   46.64 334.903  < 2e-16 ***
Residuals      463  64.48    0.14                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Por que obtenho valores-p diferentes, dependendo da ordem em que as variáveis ​​/ fatores / covariáveis ​​/ efeitos fixos (?) São codificados? Existe uma maneira de "corrigir" isso? Será que estou usando o modelo errado? Eu ainda sou muito novo na R, por isso, se você puder me ajudar com isso, mantenha-o realmente simples para que eu possa entender a resposta hehe ... Obrigado, espero que alguém possa me ajudar a entender isso!

r anova Rbeginner
fonte
3
Forneça alguns dados de amostra para Populationfulltornar seu problema reproduzível . Isso não acontece com o exemplo da aov()página de ajuda. summary(aov(yield ~ block + N + P + K, npk)); summary(aov(yield ~ K + P + block + N , npk))
13136 MrFlick
Os valores p estão mudando porque todo o campo de valores mudou. sua primeira corrida Earnings 1 109.01 109.01 782.838 < 2e-16 ***sua segunda corrida Earnings 1 46.64 46.64 334.903 < 2e-16 ***. Seus resultados não são os mesmos. Comece verificando se você não fez mais do que reorganizar variáveis.
11
TAMBÉM: no segundo modelo, você usa Earn, não Earnings ... se houver duas variáveis ​​de nomes diferentes, pode ser um problema se esse não for um erro de digitação ao copiar para o espaço de perguntas do SO.
Sim, os valores mudam, mas por quê? Eu usei exatamente as mesmas colunas do mesmo quadro de dados em ambos os modelos (a coisa Earnings vs Earn no segundo modelo é apenas porque eu escrevi errado, eu o corrigi agora).
Rbeginner
11
Isso está acontecendo porque você tem um design desequilibrado. Você encontrará muita ajuda sobre isso se pesquisar neste fórum ou apenas no Google "ANOVA desequilibrado em R". Eu recomendo olhar para o carpacote - ele implementa ANOVA Tipo II e Tipo III, que não depende da ordem das variáveis, enquanto aovque ANOVA Tipo I.
Slow loris 15/05

Respostas:

15

O problema vem da maneira como aov()faz seu teste de significância padrão. Ele usa o que é chamado de análise ANOVA "Tipo I", na qual o teste é realizado na ordem em que as variáveis ​​são especificadas no seu modelo. Portanto, no primeiro exemplo, ele determina quanta variação é explicada sexe testa sua significância; em seguida, qual parte da variação restante é explicada DMRT3e testa sua significância em termos dessa variação restante , e assim por diante. No segundo exemplo, DMRT3é avaliado apenas depois de Voltsec, Autosece sex, nessa ordem, para que haja menos variação restante DMRT3a ser explicada.

Se duas variáveis ​​preditivas forem correlacionadas, a primeira inserida no modelo receberá "crédito" total, deixando menos variação a ser "explicada" pela segunda, que pode parecer menos "estatisticamente significativa" do que a primeira, mesmo que seja não, funcionalmente. Esta pergunta e sua resposta explicam os diferentes tipos de análises ANOVA.

Uma maneira de contornar isso é extrair-se das restrições da ANOVA clássica e usar um modelo linear simples, com lm()em R, e não aov(). Isso efetivamente analisa todos os preditores em paralelo, "corrigindo" todos os preditores de uma só vez. Nesse caso, dois preditores correlacionados podem acabar tendo grandes erros padrão de seus coeficientes de regressão estimados e seus coeficientes podem diferir entre amostras diferentes da população, mas pelo menos a ordem em que você insere as variáveis ​​na especificação do modelo não importa.

Se a sua variável de resposta for algum tipo de variável de contagem, como o nome Startssugere, provavelmente você não deve usar ANOVA de qualquer maneira, pois é improvável que os resíduos sejam distribuídos normalmente, conforme a interpretação do valor p . As variáveis ​​de contagem são melhor tratadas com modelos lineares generalizados (por exemplo, glm()em R), que podem ser pensados ​​como uma generalização lm()para outros tipos de estruturas de erro residual.

EdM
fonte