Using the anova function in R with more than two models My question is about the anova function from the stats package in R. The function produces unexpected results (at least unexpected for me) when I use it with more than two models. To illustrate my question, I first simulate some data:
set.seed(363)
n <- 10
x <- runif(n, min = 0, max = 80)
y <- 1000 + 15*x - 0.3*x^2 + rnorm(n, mean = 0, sd = 100)
dat <- data.frame(y = y, x = x, x2 = x^2)

I then fit three linear models to these data:
M1 <- lm(y ~ 1, data = dat)
M2 <- lm(y ~ 1 + x, data = dat)
M3 <- lm(y ~ 1 + x + x2, data = dat)

Next, I use the anova function to compare the fitted models by means of $F$-tests. From the documentation of the anova function I conclude that I can pass multiple (i.e., more than two) fitted models to the function. So, I do the following:
anv <- anova(M1, M2, M3)

I expected that the $F$-tests from the anova function are equivalent to the $t$-tests from the summaries of the fitted models (the comments show the results):
anv$`Pr(>F)`  # NA 0.001468653 0.010810972
summary(M2)$coefficients[2, 4]  # 0.0109458
summary(M3)$coefficients[3, 4]  # 0.01081097

From the results it can be seen that the $p$-values from the $F$- and $t$-test for model M2 are not the same (0.001468653 vs. 0.0109458), whereas for model M3 they are identical (0.01081097). I looked for the reason and found that the $F$ value for M2 is calculated using the RSS and Res.Df from M3:
anv$F[2]  # 25.5709
anv$`Sum of Sq`[2] / anv$Df[2] / (anv$RSS[3] / anv$Res.Df[3])  # 25.5709

My question is the following: Why does the anova function use the RSS and Res.Df from M3 for testing M1 against M2? Is this commonly done? I checked some books on regression and could not find information about this. I also checked a similar analysis in SPSS, using multiple blocks and the R squared change option in the linear regression module. In the SPSS output, the $p$-values from the $F$-tests are in fact identical to the $p$-values from the $t$-tests. Finally, I checked the behavior of the anova function when passing it only models M1 and M2:
anova(M1, M2)$`Pr(>F)`[2]  # 0.0109458
summary(M2)$coefficients[2, 4]  # 0.0109458

In this case the $p$-value from the $F$-test is identical to the $p$-value from the $t$-test.

I think this is either due to (I) anova() acting a little weird in selecting the reference model when comparing more than two models at a time (it tries to automatically select the smallest model as the reference but sometimes fails at that), which you can maybe check by looking at dfs, or (II) due to different sums of squares being used. The t-test in lm uses Type III while anova uses Type I by default (Interpreting output from anova() when using lm() as input).
At any rate, comparing models pairwise is usually a lot more transparent and save ;).
