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

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!

## migrated from stackoverflow.comMay 14 '16 at 7:45

This question came from our site for professional and enthusiast programmers.

• Please provide some sample data for Populationfull to make your problem reproducible. This does not happen with the example from the aov() help page. summary(aov(yield ~ block + N + P + K, npk)); summary(aov(yield ~ K + P + block + N , npk)) – MrFlick May 13 '16 at 14:38
• The p values are changing because the entire field of values has changed. your first run Earnings 1 109.01 109.01 782.838 < 2e-16 *** your second runEarnings 1 46.64 46.64 334.903 < 2e-16 ***. Your results are not the same. Begin by checking to see that you have not done more than reorder variables. – bethanyP May 13 '16 at 20:41
• ALSO: in the second model, you use Earn, not Earnings...if there are two variables of different names it could be an issue if that is not a typo in copying to the S.O questions space. – bethanyP May 13 '16 at 20:49
• Yes, the values change but why? I have used the exact same columns from the exact same data frame in both models (the Earnings vs Earn thing in the second model is just that I wrote it wrong, I have corrected it now). – Rbeginner May 15 '16 at 17:56
• This is happening because you have an unbalanced design. You will find a lot of help about this if you search this forum or just Google "unbalanced ANOVA in R". I'd recommend looking into the car package- it implements Type II and Type III ANOVA, which don't depend on the order of variables, whereas aov does Type I ANOVA. – Slow loris May 15 '16 at 18:26

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