R - Analysis of Covariance

R - Modern Analysis of Covariance (ANCOVA)

Analysis of Covariance (ANCOVA) is a statistical method used to evaluate the effect of a continuous predictor variable on a continuous response variable, while also controlling for the effect of one or more categorical variables (called covariates). It’s essentially a blend of Analysis of Variance (ANOVA) and linear regression.

For instance, you might want to model how a car’s horsepower (hp) affects its mileage (mpg). However, you suspect that the transmission type (automatic or manual) also plays a significant role. ANCOVA allows you to build a model that includes both hp as a continuous predictor and am (transmission type) as a categorical covariate, and even to test if the relationship between hp and mpg is different for automatic and manual cars (this is called an interaction effect).

Modern Setup: The Tidyverse

While base R is fully capable, modern R workflows often leverage the tidyverse, a collection of packages designed for data science that share a common design philosophy. We’ll use dplyr for data manipulation and broom for cleaning up model outputs.

# Install necessary packages if you haven't already
# install.packages(c("tidyverse", "broom"))

library(tidyverse)
library(broom)

Step 1: Data Preparation

We’ll use the built-in mtcars dataset. Our first step is to select the variables of interest and ensure our categorical variable (am) is properly formatted as a factor. This gives it meaningful labels and ensures statistical functions treat it correctly.

# Load and prepare the data using dplyr pipes
input_data <- mtcars %>%
  select(mpg, hp, am) %>%
  mutate(am = factor(am, labels = c("Automatic", "Manual")))

# Display the first few rows of our prepared data
head(input_data)

When we execute the above code, it produces the following result:

                   mpg  hp        am
Mazda RX4         21.0 110    Manual
Mazda RX4 Wag     21.0 110    Manual
Datsun 710        22.8  93    Manual
Hornet 4 Drive    21.4 110 Automatic
Hornet Sportabout 18.7 175 Automatic
Valiant           18.1 105 Automatic

Step 2: ANCOVA Modeling

We will build two models to test our hypothesis:

1. Interaction Model (mpg ~ hp * am): This model tests if the relationship between hp and mpg is different for each transmission type. The * automatically includes hp, am, and their interaction term (hp:am).
2. Parallel Slopes Model (mpg ~ hp + am): This model assumes the effect of hp on mpg is the same for both transmission types, but allows for a different baseline mpg for each.

Model 1: With Interaction

Let’s build and examine the model that includes an interaction term.

# Create the ANCOVA model with an interaction term
model_interaction <- aov(mpg ~ hp * am, data = input_data)

# Use broom::tidy() for a clean summary table
summary_interaction <- tidy(model_interaction)
print(summary_interaction)

This produces a clean, tibble-based output:

# A tibble: 4 × 6
  term        df    sumsq   meansq statistic      p.value
  <chr>    <dbl>    <dbl>    <dbl>     <dbl>        <dbl>
1 hp           1   678.    678.       77.4   0.00000000150
2 am           1   202.    202.       23.1   0.0000475
3 hp:am        1     0.005   0.005     0.001  0.981
4 Residuals   28   245.      8.77      NA           NA

Interpretation: The p.value for hp and am are both very small (<< 0.05), indicating they are significant predictors of mpg. However, the interaction term hp:am has a very high p-value (0.981), suggesting that the effect of horsepower on mileage does not significantly differ between automatic and manual cars.

Model 2: Without Interaction (Parallel Slopes)

Given the insignificant interaction, a simpler model without it is likely better.

# Create the model without an interaction term
model_parallel <- aov(mpg ~ hp + am, data = input_data)

# Get the tidy summary
summary_parallel <- tidy(model_parallel)
print(summary_parallel)

This provides the summary for the simpler model:

# A tibble: 3 × 6
  term        df   sumsq  meansq statistic      p.value
  <chr>    <dbl>   <dbl>   <dbl>     <dbl>        <dbl>
1 hp           1   678.   678.       80.2   0.000000000763
2 am           1   202.   202.       23.9   0.0000346
3 Residuals   29   245.     8.46      NA           NA

Interpretation: Both hp and am remain highly significant predictors.

Step 3: Comparing Models

To formally decide which model is better, we can use the anova() function to compare them. A non-significant p-value here means the simpler model (without the interaction) is sufficient.

# Compare the two models
anova(model_interaction, model_parallel)

The output confirms our suspicion:

Analysis of Variance Table

Model 1: mpg ~ hp * am
Model 2: mpg ~ hp + am
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     28 245.43
2     29 245.44 -1 -0.0052515 0.0006 0.9806

With a Pr(>F) of 0.9806, we confirm that adding the interaction term does not significantly improve the model. We should proceed with the simpler model_parallel.

Best Practices & Next Steps

• Check Assumptions: For any linear model, you must check assumptions like normality of residuals and homogeneity of variance. You can do this with plot(model_parallel) or more elegantly with ggplot2.
• Practical Interpretation: The final model mpg ~ hp + am suggests that for a given horsepower, manual transmission cars have a different average mpg than automatic ones, and this difference is consistent across all horsepower values.
• Further Resources: For more advanced modeling, explore the tidymodels framework, and for post-hoc analysis, the emmeans package is the industry standard.