I have a data set of 482 observations.


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)

                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)


                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!

  • 3
    $\begingroup$ 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)) $\endgroup$
    – MrFlick
    Commented May 13, 2016 at 14:38
  • $\begingroup$ 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. $\endgroup$
    – bethanyP
    Commented May 13, 2016 at 20:41
  • 1
    $\begingroup$ 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. $\endgroup$
    – bethanyP
    Commented May 13, 2016 at 20:49
  • $\begingroup$ 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). $\endgroup$
    – Rbeginner
    Commented May 15, 2016 at 17:56
  • 1
    $\begingroup$ 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. $\endgroup$
    – SlowLoris
    Commented May 15, 2016 at 18:26

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.