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