Using the anova function in R with more than two models

Question

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.

Answer 1

This is a deliberate choice by the authors of the program. From the help page for anova.lm, with respect to calling it with multiple models:

Normally the F statistic is most appropriate, which compares the mean square for a row to the residual sum of squares for the largest model considered. (Emphasis added.)

When you specify only M1 and M2 in your call, M2 is "the largest model considered" and you get the result you expected. I can't speak for the authors, but I suppose with more than 2 nested models you could argue that the residual sum of squares in the largest model is the best estimate of the true residual variance and should thus be used for all comparisons.

The Type I versus Type II versus Type III issue in ANOVA, suggested in another answer, isn't responsible for this particular observation. It's important if you call anova() on a single model, however. Here's the result for the full model M3 as specified in the OP:

> anova(M3)
Analysis of Variance Table

Response: y
          Df Sum Sq Mean Sq F value   Pr(>F)   
x          1 604116  604116  25.571 0.001469 
x2         1 279861  279861  11.846 0.010811  
Residuals  7 165376   23625 

Now switch the order of entry of the x and x2 terms in the call to lm():

 M3s <- lm(formula = y ~ 1 + x2 + x, data = dat)

and call anova() on that:

> anova(M3s)
Analysis of Variance Table

Response: y
          Df Sum Sq Mean Sq F value    Pr(>F)    
x2         1 816194  816194 34.5477 0.0006131
x          1  67783   67783  2.8691 0.1341244    
Residuals  7 165376   23625

Compare against the standard summary of M3:

> summary(M3)
## some lines omitted from output
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1101.70266  119.17661   9.244 3.58e-05 
x             11.42016    6.74218   1.694   0.1341    
x2            -0.29281    0.08508  -3.442   0.0108   

You'll note that the p-values in the summary() output for the predictors are those reported by anova() when the predictor was listed last in the call to lm().

Answer 2

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 ;).