I need to run a Type III ANCOVA with an interaction term between a categorical and numeric variable in R. I have found three different sets of sample code to do this, and they yield identical results when I run ANCOVA without an interaction, but when I add the interaction the results diverge.
Below is some sample code for the three approaches I am comparing. First, I will show models without interactions, which yield identical results, and then I will add the interaction term, getting inconsistent results.
# Preparing data
library(MASS)
data(mtcars)
data <- mtcars
data$vs <- as.factor (data$vs)
Models without interactions
Model 1 is based on https://www.datanovia.com/en/lessons/ancova-in-r/#computation-1
library(tidyverse)
library(rstatix)
aov1 <- data %>%
anova_test(mpg ~ vs + cyl)
> get_anova_table(aov1)
ANOVA Table (type II tests)
Effect DFn DFd F p p<.05 ges
1 vs 1 29 0.226 6.38e-01 0.008
2 cyl 1 29 30.669 5.70e-06 * 0.514
Model 2 is based on https://www.r-bloggers.com/2021/07/how-to-perform-ancova-in-r/
library(car)
aov2 <- aov (mpg ~ vs + cyl, data = data)
> Anova (aov2, type ="III")
Anova Table (Type III tests)
Response: mpg
Sum Sq Df F value Pr(>F)
(Intercept) 928.17 1 87.9765 2.765e-10 ***
vs 2.38 1 0.2255 0.6384
cyl 323.56 1 30.6691 5.695e-06 ***
Residuals 305.96 29
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Model 3 is based on https://www.r-bloggers.com/2021/10/analysis-of-covariance-ancova-using-r/
library(rstatix)
aov3 <- anova_test(data = data, formula = mpg ~ vs + cyl, type = 3, detailed = TRUE)
> aov3
ANOVA Table (type III tests)
Effect SSn SSd DFn DFd F p p<.05 ges
1 (Intercept) 1397.137 305.955 1 29 132.428 2.49e-12 * 0.820
2 vs 2.379 305.955 1 29 0.226 6.38e-01 0.008
3 cyl 323.564 305.955 1 29 30.669 5.70e-06 * 0.514
No concerns for now. All the F scores are consistent across the three models.
When I add an interaction between vs and cyl, however, the results diverge.
Models with interactions
Model 1:
library(tidyverse)
library(rstatix)
aov1_int <- data %>%
anova_test(mpg ~ vs * cyl)
> get_anova_table(aov1_int)
ANOVA Table (type II tests)
Effect DFn DFd F p p<.05 ges
1 vs 1 28 0.224 6.40e-01 0.008
2 cyl 1 28 30.482 6.69e-06 * 0.521
3 vs:cyl 1 28 0.823 3.72e-01 0.029
Model 2:
library(car)
aov2_int <- aov (mpg ~ vs * cyl, data = data)
> Anova (aov2_int, type ="III")
Anova Table (Type III tests)
Response: mpg
Sum Sq Df F value Pr(>F)
(Intercept) 540.09 1 50.8803 9.219e-08 ***
vs 5.68 1 0.5351 0.4705467
cyl 167.06 1 15.7381 0.0004589 ***
vs:cyl 8.74 1 0.8233 0.3719539
Residuals 297.22 28
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Model 3:
library(rstatix)
aov3_int <- anova_test(data = data, formula = mpg ~ vs * cyl, type = 3, detailed = TRUE)
> aov3_int
ANOVA Table (type III tests)
Effect SSn SSd DFn DFd F p p<.05 ges
1 (Intercept) 1405.720 297.216 1 28 132.430 3.98e-12 * 0.825
2 vs 5.680 297.216 1 28 0.535 4.71e-01 0.019
3 cyl 322.975 297.216 1 28 30.427 6.78e-06 * 0.521
4 vs:cyl 8.739 297.216 1 28 0.823 3.72e-01 0.029
Summary of results
- When comparing Model 1 and 2, the F scores are the same for the interaction term but they differ for
vs
andcyl
. - When comparing Model 2 and Model 3, the F scores are the same for
vs
and the interaction term, but they differ forcyl
. - When comparing Model 1 and Model 3, the F scores are the same for
cyl
and the interaction term, but they differ forvs
.
I would be grateful if anyone could help me make sense of these inconsistent results and suggest which results I should trust and report.
Many thanks everyone for your attention and support!