What is the explanation of the difference caused by the order of variables in anova for the concrete example in R?

Question

I wonder the real reason behind the difference in results of the two anovas that use the same covariates. However, the results are different.

library(PASWR2)
head(HSWRESTLER); tail(HSWRESTLER)
#  age    ht    wt abs triceps subscap hwfat tanfat skfat
# 1  18 65.75 133.6   8       6    10.5 10.71   11.9  9.80
# 2  15 65.50 129.0  10       8     9.0  8.53   10.0 10.56
# ...
# 77  15 68 153.8  13       7      11 10.07   16.7 11.77
# 78  15 66 258.6  45      37      43 33.75   34.5 38.93

mod1.HSW <- lm(hwfat ~ abs + triceps + subscap, data = HSWRESTLER)
anova(mod1.HSW)
# Analysis of Variance Table
# Response: hwfat              
#            Df Sum Sq Mean Sq F value    Pr(>F)     
# abs        1 5072.8  5072.8 535.858 < 2.2e-16 ***    
# triceps    1  242.2   242.2  25.581 2.984e-06 ***    
# subscap    1    2.2     2.2   0.237    0.6278        
# Residuals 74  700.5     9.5   
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

mod2.HSW <- lm(hwfat ~ subscap + triceps + abs, data = HSWRESTLER)
anova(mod2.HSW) # ANOVA
# Analysis of Variance Table
# Response: hwfat
#           Df Sum Sq Mean Sq F value    Pr(>F)    
# subscap    1 4939.0  4939.0 521.720 < 2.2e-16 ***
# triceps    1  204.6   204.6  21.616 1.422e-05 ***
# abs        1  173.6   173.6  18.341 5.473e-05 ***
# Residuals 74  700.5     9.5                      
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

?anova help file is not explanatory enough about what is going on. The similar SOF questions does not handle the situation with a concrete example like this. The variables in the regressions seem to be continuous variables (covariates). It seems the variable order is important. But, what does that order mean? What can be inferred from the above two anovas? Any idea?

Thx etienne. I saw that post. The explanations there (lacking a concrete example) are very subtle to me. It would be nice if the problem is handled in detail. I'm sure there would be many to benefit from the answer, me first of course. — Erdogan CEVHER, Nov 23 '16 at 21:17
Yeah, I saw that one also. There, there is the case lm(y ~ V1) and lm(y ~ V1+V2) instead of lm(y~V1+V2) and lm(y~V2+V1). Also, there is no std function in explanation; i.e., there, anova, Anova, manova is used. — Erdogan CEVHER, Nov 23 '16 at 21:34
Note that this has, not surprisingly, already been asked and answered on this site http://stats.stackexchange.com/questions/212496/why-do-p-values-change-in-significance-when-changing-the-order-of-covariates-in?rq=1 — mdewey, Nov 24 '16 at 12:06
mdewey, I just saw the link you provided now. Thx. By the way, the link you provided helped me better understood the one here: http://stats.stackexchange.com/questions/20452/how-to-interpret-type-i-type-ii-and-type-iii-anova-and-manova — Erdogan CEVHER, Nov 24 '16 at 16:56

score 4 · Accepted Answer · answered Nov 23 '16 at 22:10

When you call anova on an lm model fit, under the hood you are really using ?anova.lm, which according to the documentation "gives a sequential analysis of variance table for that fit." This is a type I ANOVA, where order of the variables matter.The term abs in your second example only represents the unique portion of the regression explained given the previous two variables.

You can perform type II ANOVA using drop1(). Here order doesn't matter, thus each main effect can be understood as the individual contribution of that predictor:

> drop1(mod1.HSW)
#Single term deletions
#
#Model:
#hwfat ~ abs + triceps + subscap
#        Df Sum of Sq    RSS    AIC
#<none>               700.54 179.22
#abs      1   173.629 874.17 194.49
#triceps  1   111.837 812.38 188.77
#subscap  1     2.244 702.78 177.47

> drop1(mod2.HSW)
#Single term deletions
#
#Model:
#hwfat ~ subscap + triceps + abs
#        Df Sum of Sq    RSS    AIC
#<none>               700.54 179.22
#subscap  1     2.244 702.78 177.47
#triceps  1   111.837 812.38 188.77
#abs      1   173.629 874.17 194.49

score 2 · Answer 2 · edited May 23 '17 at 12:39

I like to understand linear regression in terms of its computation. A linear model is computed by QR factorization, giving a number of "effects" on the rotated plane. You can refer to What is the “effects” returned by aov and lm? for details.

"Effects" is the key for generating ANOVA table:

e1 <- mod1.HSW$effects[2:4]  ## we also exclude "Intercept"
#       abs     triceps     subscap 
# 71.223840   15.561899   -1.497877 

e2 <- mod2.HSW$effects[1:4]  ## we also exclude "Intercept"
#  subscap     triceps         abs 
# 70.27797    14.30491    13.17682

Note, effects are already different, due to order of terms. Two models are equivalent, so they have the same residual variance:

d <- with(mod1.HSW, c(crossprod(residual)) / df.residual)
# [1] 9.466757

ANOVA table and F-statistic can be compute as:

cbind("explained-variance" = e1 ^ 2, "F-statistic" = e1 ^ 2 / d)

#            explained-variance  F-statistic
#abs                5072.835390  535.8577814
#triceps             242.172704   25.5813797
#subscap               2.243634    0.2370014

cbind("explained-variance" = e2 ^ 2, "F-statistic" = e2 ^ 2 / d)

#            explained-variance F-statistic
#subscap              4938.9927   521.71960
#triceps               204.6305    21.61569
#abs                   173.6286    18.34087

Have a comparison with anova(mod1.HSW) and anova(mod2.HSW) yourself.

To this stage, you should realize that the difference in "effects" explains why ANOVA table is different. lm performs non-pivoted QR factorization, hence is sensitive to term orders. QR factorization for a matrix X is like a "Gram-Schmidt" process. If you have 2 columns linearly correlated, then reversing these two columns will make QR factorization different.

Now have a look at your data:

with(HSWRESTLER, plot(abs, subscap))

Isn't this a linear relationship?

What is the explanation of the difference caused by the order of variables in anova for the concrete example in R?

2 Answers2