# R: comparing three lm models using anova gives unexpected F-value

First, I generated some data, fitted three models, and compared them with each other using the anova function:

# Generate data
set.seed(100)
intercept <- rnorm(100) + 1
a <- rnorm(100) * 2.5
b <- rnorm(100) * .5
c <- rnorm(100) * .1
y <- intercept + a + b + c

# Fit models
fit0 <- lm(y ~ 1)
fit1 <- lm(y ~ 1 + a)
fit2 <- lm(y ~ 1 + a + b + c)

# Compare models
comp <- stats::anova(fit0, fit1, fit2)


When I checked the comp model comparison object it gave me an unexpected F-value for the comparison between fit0 and fit1. While I expected that it would give me the same F-value as the overall model fit F-test that you get when calling summary(fit1) (because the overall model fit test also compares the model with the null model), it didn't.

summary(fit1) gave me:

Residual standard error: 1.173 on 98 degrees of freedom
Multiple R-squared:  0.6856,    Adjusted R-squared:  0.6824
F-statistic: 213.7 on 1 and 98 DF,  p-value: < 2.2e-16


and comp gave me:

  Res.Df    RSS Df Sum of Sq      F    Pr(>F)
1     99 428.62
2     98 134.75  1   293.874 281.07 < 2.2e-16 ***
3     96 100.37  2    34.378  16.44 7.245e-07 ***


When I looked at the anova help file, ?anova.lm, I found out why this is the case:

Optionally the table can include test statistics. 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.

So apparently it uses the residual sum of squares for the largest model (in this case fit2). However, I cannot find why it does that, so my questions are:

1. Why does the anova function use the residual sum of squared for the largest model (instead of the current models')?
2. Would you recommend only comparing two models at a time with the anova function? So, as follows:

comp1 <- stats::anova(fit0, fit1)
comp2 <- stats::anova(fit1, fit2)


## 2 Answers

1. Using the largest model as the yardstick for comparison is a choice. This choice makes sense if you have reason to think that the largest model considered is the right yardstick for comparison. If the largest model considered is not a good yardstick, for whatever reason, then this choice doesn't make sense, and you should do something else.

The largest model might be a bad choice for any number of reasons. For example, it might be very close to saturated, or it might use information to which no workable model could have access.

1. It is not necessary to compare only two models at a time, but this might make sense if you don't have a large model that is a good standard for comparison.

Pairwise comparison of models might yield a different result from the simultaneous comparison of all models. For an example of this, which prompted this post, see Faraway, "Extending the Linear Model with R" (2ed), Chapter 6, Exercise 3 (d).

To see intuitively that the pairwise comparison can yield a different answer from a simultaneous comparison of 3 models, let's call our models m1, m2, m3, in order from smallest to largest. One of the simultaneous comparisons will be a comparison of m1 to m2 with m3 as the yardstick. Such an F-statistic will never appear in the 3 possible pairwise comparisons.

There are many different ways to compare different models and it is a bit arbitrary. For example, think about the type I/II/III sums of squares.

What the stats::anova() command does is compute type I sums. The total variance is split up into 3 components from residual and the cascading models. Then those parts are compared.

$$\underbrace{428.62}_{total} = \underbrace{100.37}_{residual} + \underbrace{34.378}_{b+c} + \underbrace{293.874}_{a}$$

With car::Anova() you could also compute the type II and type III sums. In your case, if you want to have something resemble, summary(fit), then you need to do it manually.

Why type I/II/III sums instead of manual pairwise comparisons? A disadvantage of manual pairwise comparisons might be a lower power. If the residual sum of squares contains partly variance from a model then it influences the comparison. For example consider

Analysis of Variance Table

Model 1: y ~ 1
Model 2: y ~ 1 + b + c
Res.Df    RSS Df Sum of Sq     F Pr(>F)
1     99 428.62
2     97 421.27  2    7.3571 0.847 0.4318


This has a low significance because the residuals contains a large variance due to the 'a' factor. We regard the 421.27 RSS as random noise. But, that random variance can be for a large part deterministic when we know 'a'. So if you know 'a' it might be better to include it

Analysis of Variance Table

Model 1: y ~ 1 + a
Model 2: y ~ 1 + a + b + c
Res.Df    RSS Df Sum of Sq     F    Pr(>F)
1     98 134.75
2     96 100.37  2    34.378 16.44 7.245e-07


On the other hand, if the additional variables have no effect, then they will reduce the power of the test.