Why do p-values change in significance when changing the order of covariates in the aov model?

Question

I have a data set of 482 observations.

data=Populationfull

Im going to do a genotype association analysis for 3 SNPs. Im trying to build a model for my analysis and Im using the aov(y~x,data=...). For one trait I have several fixed effects and covariates that I have included in the model, like so:

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

I discovered that if i changed the order of the variables in the model i got different p-values, please see below.

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

Why do I get different p-values depending on in which order the variables/factors/covariates/fixedeffects(?) are coded? Is there a way to "correct" for it? Can it be that Im using the wrong model? I am still quite new at R so if you can help me with this please keep it really simple so I can understand the answer hehe... Thank you, hopefully someone can help me understand this!

Answer 1

The problem comes from the way that aov() does its default significance testing. It uses what is called "Type I" ANOVA analysis, in which testing is done in the order that the variables are specified in your model. So in the first example, it determines how much variance is explained by sex and tests its significance, then what portion of the remaining variance is explained by DMRT3 and tests its significance in terms of that remaining variance, and so forth. In the second example, DMRT3 is only evaluated after Voltsec, Autosec, and sex, in that order, so there is less variance remaining for DMRT3 to explain.

If two predictor variables are correlated then the first one entered into the model will get full "credit," leaving less variance to be "explained by" the second one, which thus may appear less "statistically significant" than the first even if it is not, functionally. This question and its answer explain the different Types of ANOVA analyses.

One way to get around this is to extract yourself from the strictures of classical ANOVA and use a simple linear model, with lm() in R, rather than aov(). This effectively analyzes all predictors in parallel, "correcting for" all predictors at once. In that case, two correlated predictors might end up having large standard errors of their estimated regression coefficients, and their coefficients might differ among different samples from the population, but at least the order you enter the variables into the model specification won't matter.

If your response variable is some type of count variable, as its name Starts suggests, then you probably shouldn't be using ANOVA anyway as residuals are unlikely to be normally distributed, as the p-value interpretation requires. Count variables are better handled with generalized linear models (e.g., glm() in R), which can be thought of as a generalization of lm() for other types of residual error structures.