Why does the Type-3 ANOVA using LRT via car::Anova() give different result than term-by-term LRT model comparison via anova() in R?

Question

With this very simple data:

> A
 [1] "a" "a" "a" "a" "b" "b" "b" "b" "b" "b" "b" "b"
> B
 [1] "x" "y" "x" "y" "x" "y" "x" "y" "x" "x" "x" "x"
> C
 [1] "l" "l" "m" "m" "l" "l" "m" "m" "l" "l" "l" "l"
> response
 [1] 14 30 15 35 50 51 30 32 51 55 53 55

I try to reproduce the Type-3 car::Anova by using step-by-step term elimination to better understand interactions and analysis of variance.

For example I want to assess the term "C"

> options(contrasts = c("contr.sum", "contr.poly"))
> m1 <- lm(response ~ A*B*C)  # the full model
> car::Anova(m1, type=3, test.statistic = "LR")
Anova Table (Type III tests)
Response: response
            Sum Sq Df F value  Pr(>F)

(Intercept)   9374  1 1802.78 1.8e-06 ***
A              716  1  137.69  0.0003 ***
B              182  1   35.00  0.0041 **
C              178  1   34.23  0.0043 **
A:B            178  1   34.23  0.0043 ** 
A:C            317  1   61.03  0.0014 ** 
B:C              8  1    1.63  0.2714

A:B:C            0  1    0.00  0.9755

Residuals       21  4

Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

For term "C" I got p-value = 0.0043 Now I going to assess term C by elimination:

> m2 <- lm(response ~ A + B + A:B)  # the term "C", eliminated all terms with C
> anova(m1, m2)
Analysis of Variance Table
Model 1: response ~ A * B * C
Model 2: response ~ A + B + A:B
  Res.Df RSS Df Sum of Sq    F Pr(>F)

1      4  21

2      8 648 -4      -627 30.1  0.003 **

Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Here p-value is = 0.003. Close but not the same. This is the simplest general linear model, no mixed effects, no transformations.

When I typed "test = LRT"

> anova(m1, m2, test="LRT")
Analysis of Variance Table
Model 1: response ~ A * B * C
Model 2: response ~ A + B + A:B
  Res.Df RSS Df Sum of Sq Pr(>Chi)

1      4  21

2      8 648 -4      -627   <2e-16 ***

Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The result is totally different.

Please note, my goal IS NOT to make any formal inference, only to UNDERSTAND the calculations behind. I understand, that elimination terms from model should be equivalent to the car::Anova() and give equal result numerically.

Let's confirm:

> drop1(m1, scope = ~A*B*C, test="F")
Single term deletions
Model:
response ~ A * B * C
       Df Sum of Sq RSS  AIC F value Pr(>F)

<none>               21 22.6

A       1       716 737 63.4  137.69 0.0003 ***
B       1       182 203 47.9   35.00 0.0041 **
C       1       178 199 47.7   34.23 0.0043 **
A:B     1       178 199 47.7   34.23 0.0043 ** 
A:C     1       317 338 54.1   61.03 0.0014 ** 
B:C     1         8  29 24.7    1.63 0.2714

A:B:C   1         0  21 20.6    0.00 0.9755

Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

It works! It agrees with car::Anova().

So how should I eliminate the terms to obtain the same result using anova()? Evidently comparing A+B+A:B vs. the full model is not enough!

Answer 1

For a test of C, Type III sum of squares compare the full model with a reduced model where just the main effect of C is removed but all interactions involving C are retained (see this excellent post for a detailled comparison between the different types of SS). In other words, Type III SS do not respect the principle of marginality.

In order to reproduce the output from car::Anova using anova, you need to fit the two models first and then feed them into anova. The problem is that with only categorical variables (i.e. factors), it's not possible to drop just the main effect of C using the formula interface when interactions with C are present in the model. But you can do it by manipulating the model matrix directly (see here). Save the model matrix of the full model and remove the column that corresponds to C, then refit the model to obtain the reduced model.

Here is how it's done:

library(car)
options(contrasts = c("contr.sum", "contr.poly"))

The data

A <- factor(rep(c("a", "b"), c(4, 8)))
B <- factor(c("x", "y", "x", "y", "x", "y", "x", "y", "x", "x", "x", "x"))
C <- factor(c("l", "l", "m", "m", "l", "l", "m", "m", "l", "l", "l", "l"))
y <- c(14, 30, 15, 35, 50, 51, 30, 32, 51, 55, 53, 55)
dat <- data.frame(y, A, B , C)
m <- lm(y~ABC, data = dat)
Anova(m, type = "III")
Response: y
            Sum Sq Df   F value    Pr(>F)

(Intercept) 9374.5  1 1802.7788 1.839e-06 ***
A            716.0  1  137.6934 0.0003017 ***
B            182.0  1   35.0011 0.0040878 ** 
C            178.0  1   34.2318 0.0042571 ** 
A:B          178.0  1   34.2318 0.0042571 ** 
A:C          317.3  1   61.0267 0.0014491 ** 
B:C            8.5  1    1.6250 0.2714108

A:B:C          0.0  1    0.0011 0.9754909

Residuals     20.8  4 

So far so good. To reproduce the output for C, fit the full model, extract the model matrix, remove column corresponding to C and refit.

m1 <- lm(y~A*B*C, data = dat)  # Full model
X <- model.matrix(m1) # Extract the model matrix of the full model
X0 <- X[, -c(4)] # Removing the column corresponding to C
m0 <- with(dat, lm(y~X0 + 0)) # Reduced model
anova(m0, m1, test = "F")
Model 1: y ~ X0 + 0
Model 2: y ~ A * B * C
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)

1      5 198.81

2      4  20.80  1    178.01 34.232 0.004257 **

You see that now we fully reproduced the output from Anova for the main effect of C using the two models and anova.