Pourquoi les valeurs de p changent-elles de signification lors du changement de l'ordre des covariables dans le modèle aov?

10

J'ai un ensemble de données de 482 observations.

data=Populationfull

Je vais faire une analyse d'association de génotypes pour 3 SNP. Im essayant de construire un modèle pour mon analyse et Im utilisant l'aov (y ~ x, data = ...). Pour un trait, j'ai plusieurs effets fixes et covariables que j'ai inclus dans le modèle, comme ceci:

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

J'ai découvert que si je modifiais l'ordre des variables dans le modèle, j'obtenais différentes valeurs de p, veuillez voir ci-dessous.

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

Pourquoi est-ce que j'obtiens des valeurs de p différentes selon l'ordre dans lequel les variables / facteurs / covariables / effets fixes (?) Sont codés? Existe-t-il un moyen de "corriger" cela? Se peut-il que j'utilise le mauvais modèle? Je suis encore assez nouveau chez R, donc si vous pouvez m'aider, veuillez rester très simple pour que je puisse comprendre la réponse hehe ... Merci, j'espère que quelqu'un pourra m'aider à comprendre cela!

Rbeginner
la source
3
Veuillez fournir quelques exemples de données pour Populationfullrendre votre problème reproductible . Cela ne se produit pas avec l'exemple de la aov()page d'aide. summary(aov(yield ~ block + N + P + K, npk)); summary(aov(yield ~ K + P + block + N , npk))
MrFlick
Les valeurs p changent car tout le champ de valeurs a changé. votre première course Earnings 1 109.01 109.01 782.838 < 2e-16 ***votre deuxième course Earnings 1 46.64 46.64 334.903 < 2e-16 ***. Vos résultats ne sont pas les mêmes. Commencez par vérifier que vous n'avez pas fait plus que réorganiser les variables.
1
AUSSI: dans le deuxième modèle, vous utilisez Earn, not Earnings ... s'il y a deux variables de noms différents, cela pourrait être un problème si ce n'est pas une faute de frappe lors de la copie vers l'espace de questions SO.
Oui, les valeurs changent mais pourquoi? J'ai utilisé exactement les mêmes colonnes de la même trame de données exacte dans les deux modèles (la chose Gains vs Gains dans le deuxième modèle est juste que je l'ai mal écrit, je l'ai corrigé maintenant).
Rbeginner
1
Cela se produit parce que votre conception est déséquilibrée. Vous trouverez beaucoup d'aide à ce sujet si vous recherchez ce forum ou simplement Google "ANOVA déséquilibrée en R". Je recommanderais de regarder dans le carpackage - il implémente les ANOVA de type II et de type III, qui ne dépendent pas de l'ordre des variables, alors aovque l'ANOVA de type I.
Loris lents

Réponses:

15

Le problème vient de la façon dont il aov()effectue ses tests de signification par défaut. Il utilise ce que l'on appelle l'analyse ANOVA "Type I", dans laquelle les tests sont effectués dans l'ordre dans lequel les variables sont spécifiées dans votre modèle. Ainsi, dans le premier exemple, il détermine la variance expliquée par sexet teste sa signification, puis quelle partie de la variance restante est expliquée DMRT3et teste sa signification en termes de cette variance restante , etc. Dans le second exemple, DMRT3est seulement évaluée après Voltsec, Autosecet sex, dans cet ordre, donc il y a moins de variance restante pour DMRT3expliquer.

Si deux variables prédictives sont corrélées, alors la première entrée dans le modèle obtiendra un «crédit» complet, laissant moins de variance à être «expliquée par» la seconde, qui peut donc apparaître moins «statistiquement significative» que la première même si elle est pas, fonctionnellement. Cette question et sa réponse expliquent les différents types d'analyses ANOVA.

Une façon de contourner cela est de vous extraire des restrictions de l'ANOVA classique et d'utiliser un modèle linéaire simple, avec lm()en R, plutôt que aov(). Cela analyse efficacement tous les prédicteurs en parallèle, "corrigeant" tous les prédicteurs à la fois. Dans ce cas, deux prédicteurs corrélés pourraient finir par avoir de grandes erreurs-types de leurs coefficients de régression estimés, et leurs coefficients pourraient différer entre différents échantillons de la population, mais au moins l'ordre dans lequel vous entrez les variables dans la spécification du modèle n'aura pas d'importance.

Si votre variable de réponse est un certain type de variable de comptage, comme son nom l' Startsindique, vous ne devriez probablement pas utiliser ANOVA de toute façon car les résidus sont peu susceptibles d'être normalement distribués, comme le requiert l'interprétation de la valeur de p . Les variables de comptage sont mieux gérées avec des modèles linéaires généralisés (par exemple, glm()dans R), qui peuvent être considérés comme une généralisation de lm()pour d'autres types de structures d'erreur résiduelles.

EdM
la source