Inconsistent results for ANCOVA with interaction term in R

56 Views Asked by At

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 and cyl.
  • When comparing Model 2 and Model 3, the F scores are the same for vs and the interaction term, but they differ for cyl.
  • When comparing Model 1 and Model 3, the F scores are the same for cyl and the interaction term, but they differ for vs.

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!

0

There are 0 best solutions below