The order of variables in ANOVA matters, doesn't it?

Question

Am I correct to understand that the order in which variables are specified in a multifactorial ANOVA makes a difference but that the order does not matter when doing a multiple linear regression?

So assuming an outcome such as measured blood loss y and two categorical variables

1. adenoidectomy method a,
2. tonsillectomy method b.

The model y~a+b is different to the model y~b+a (or so my implementation in R seems to indicate).

Am I correct to understand that the term here is that ANOVA is a hierarchical model since it first attributes as much variance as it can to the first factor before trying to attribute residual variance to the second factor?

In the example above the hierarchy makes sense because I always do the adenoidectomy first before doing the tonsillectomy but what would happen if one had two variables with no inherent order?

Answer 1

This question evidently came from a study with an unbalanced two-way design, analyzed in R with the aov() function; this page provides a more recent and detailed example of this issue.

The general answer to this question, as to so many, is: "It depends." Here it depends on whether the design is balanced and, if not, which flavor of ANOVA is chosen.

First, it depends on whether the design is balanced. In the best of all possible worlds, with equal numbers of cases in all cells of a factorial design, there would be no difference due to the order of entering the factors into the model, regardless of how ANOVA is performed.* The cases at hand, evidently from a retrospective clinical cohort, seem to be from a real world where such balance was not found. So the order might matter.

Second, it depends on how the ANOVA is performed, which is a somewhat contentious issue. The types of ANOVA for unbalanced designs differ in the order of evaluating main effects and interactions. Evaluating interactions is fundamental to two-way and higher-order ANOVA, so there are disputes over the best way to proceed. See this Cross Validated page for one explanation and discussion. See the Details and the Warning for the Anova() (with a capital "A") function in the manual for the car package for a different view.

The order of factors does matter in unbalanced designs under the default aov() in R, which uses what are called type-I tests. These are sequential attributions of variance to factors in the order of entry into the model, as the present question envisioned. The order does not matter with the type-II or type-III tests provided by the Anova() function in the car package in R. These alternatives, however, have their own potential disadvantages noted in the above links.

Finally, consider the relation to multiple linear regression as with lm() in R, which is essentially the same type of model if you include interaction terms. The order of entry of variables in lm() does not matter in terms of regression coefficients and p-values reported by summary(lm()), in which a k-level categorical factor is coded as (k-1) binary dummy variables and a regression coefficient is reported for each dummy.

It is, however, possible to wrap the lm() output with anova() (lower-case "a," from the R stats package) or Anova() to summarize the influence of each factor over all of its levels, as one expects in classical ANOVA. Then the ordering of factors will matter with anova() as for aov(), and will not matter with Anova(). Similarly, the disputes over which type of ANOVA to use would return. So it's not safe to assume order-independence of factor entry with all downstream uses of lm() models.

---

*Having equal numbers of observations in all cells is sufficient but, as I understand it, not necessary for the order of factors to be irrelevant. Less demanding types of balance may allow for order-independence.

Answer 2

The term hierarchical model refers to the structure between the factors. For example, a multi-center study is hierarchical: You have the patients nested within the hospitals treating them. Each hospital treats patients with placebo and verum, but recieving each of them in either hospital A or B is slightly different due to some common effect of the hospital governing on all their patients (might even be an interaction effect with the experimental agent). So it's called hierarchical effect.

Now your ectomy methods may be hierarchical: Is it plausible that a certain tonsillectomy method is slightly different (in itself, not yet in the effect, because that's what your are going to estimate and test) depending on the adenoidectomy method used prior on the same patient? If yes, you should specify it in your model.

Your observation that y~a+b may be different from y~b+a indicates that there is something wrong. Additive effects commute, so there should not be a difference (apart from small numerical differences). It is neither plausible nor desired that the effect of the surgery methods may depend on the order in which the statistician later specifies the effects. So you probably picked the wrong approach to feed R with the data.