Activity: ANCOVA

Covariates

Hands-on activity: analysis of covariance, using cricket chirping data, real fruitfly longevity data (Partridge), and a heterogeneous-slopes sea urchin example.

Author

Bill Perry

Worksheet: ANCOVA — Analysis of Covariance

How to use this worksheet

Work through each part in order, at your own pace. Type every line of code yourself into a plain R script — do not copy-paste. Blocks marked ▶ Run this are code you should type and execute. Blocks marked ✏️ Your turn ask you to write, modify, or answer something. Boxes marked 🚀 If you finish early are optional bonus material that goes a bit further than what we covered in lecture.

This worksheet uses simulated cricket chirping data, the real Partridge fruitfly-longevity dataset, and simulated sea urchin data (for a heterogeneous-slopes example).

▶ Run this in your Script:

library(lmtest)
library(interactions)
library(janitor)
library(car)
library(emmeans)
library(broom)
library(patchwork)
library(tidyverse)

options(scipen = 999)

What is ANCOVA? ANCOVA (Analysis of Covariance) combines regression and ANOVA to:

Compare group means while adjusting for a continuous covariate
Increase statistical power by reducing residual error
Control for confounding variables

When to use ANCOVA:

Response variable: Continuous
Predictor variable: Categorical (factor/groups)
Covariate: Continuous variable that affects the response

Key assumptions of ANCOVA:

Independence of observations
Normality of residuals
Homogeneity of variances across groups
Linearity between response and covariate within each group
Homogeneity of slopes (most critical!) — regression slopes must not be significantly different across all groups

Important

Critical first step: Always test for homogeneity of slopes before proceeding with ANCOVA. If slopes differ significantly between groups, standard ANCOVA is inappropriate.

Part 1 · Cricket Chirping Analysis

We want to compare chirping rate of two cricket species — Oecanthus exclamationis and Oecanthus niveus. But we measured rates at different temperatures, and there’s a relationship between pulse rate and temperature. ANCOVA lets us adjust for the temperature effect to get a more powerful test!

▶ Run this — create the simulated cricket data and visualize it:

set.seed(456)
n <- 40
spp <- rep(c("O. exclamationis", "O. niveus"), each = n/2)
temp <- c(rnorm(n/2, mean = 22, sd = 2), rnorm(n/2, mean = 24, sd = 2))
chirp_rate <- 40 + 2.5 * (temp - 23) + ifelse(spp == "O. exclamationis", 10, 0) + rnorm(n, sd = 3)
c_df <- data.frame(spp = spp, temp = temp, chirp_rate = chirp_rate)

head(c_df)

ggplot(c_df, aes(x = temp, y = chirp_rate, color = spp)) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE)

Step 1: Test homogeneity of slopes

This is the most critical assumption! We test if the regression slopes are equal across all groups.

options(contrasts = c("contr.sum", "contr.poly")) # for true type 3 anovas
# options(contrasts = c("contr.treatment", "contr.poly")) # original version
# see below for interpretation...

lm_int_model <- lm(chirp_rate ~ temp * spp, data = c_df)
Anova(lm_int_model, type = 3)

Interpretation: If p > 0.05, slopes are homogeneous — proceed with ANCOVA. If p < 0.05, slopes differ and standard ANCOVA is inappropriate.

The contrasts setting changes the specific hypothesis being tested for your main effects:

With contr.treatment, the test for temp tests the effect of temperature only for the reference species (O. exclamationis comes first alphabetically) — is mean chirp_rate equal to zero for the reference species when temp is 0?
With contr.sum (“sum coding” or “deviation coding”), the test for temp tests the average effect of temperature across both species — is average mean chirp_rate (averaged across both species) equal to zero when temp is 0?

Important

The intercept is hard to interpret in both models because temp = 0 is meaningless. If you centered the temp variable (temp_c = temp - mean(temp)) and re-ran the model, the intercept would then test the chirp_rate at the average temperature — much more interpretable.

Step 2: Fit the ANCOVA model

Since slopes are homogeneous (p > 0.05), fit the ANCOVA model without the interaction term.

ancova_model <- lm(chirp_rate ~ temp + spp, data = c_df)
summary(ancova_model)

Interpretation:

spp1 is the clue that contr.sum is active — it represents the coefficient for the first species (O. exclamationis) relative to the grand mean.
(Intercept): the average chirp_rate (across both species) when temp is 0.
temp: the common slope — for every 1-degree increase in temp, chirp_rate increases by that many units.
spp1: the deviation from the grand mean for O. exclamationis; the effect for O. niveus is the negative of that value.

Anova(ancova_model, type = 2)

Both Type II (technically more appropriate here) and Type III end up doing the exact same calculation:

The test for temp gets the Sum of Squares for temp after accounting for spp.
The test for spp gets the Sum of Squares for spp after accounting for temp.

✏️ Your turn: Based on the Anova() output, is there a significant species effect on chirping rate after adjusting for temperature? ________________________

Step 3: Check model assumptions

▶ Run this:

par(mfrow = c(2, 2))
plot(ancova_model, main = "ANCOVA Diagnostic Plots")
par(mfrow = c(1, 1))

shapiro.test(ancova_model$residuals)

leveneTest(chirp_rate ~ spp, data = c_df)

# Breusch-Pagan (BP) test for homoscedasticity
lmtest::bptest(ancova_model)

Step 4: Calculate adjusted means

ANCOVA compares adjusted means — what each group’s mean would be at the overall mean of the covariate.

c_emmeans <- emmeans(ancova_model, "spp")

cricket_adj_means_df <- as.data.frame(c_emmeans)
cricket_adj_means_df

Step 5: Pairwise comparisons

pairs(c_emmeans, adjust = "sidak")

Step 6: Visualize results

plot(c_emmeans, comparisons = TRUE)

c_emmeans_df <- as.data.frame(c_emmeans)
ggplot(c_emmeans_df, aes(x = spp, y = emmean, fill = spp)) +
  geom_bar(stat = "identity", width = 0.7) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  labs(title = "Adjusted Mean Chirping Rate by Species",
       subtitle = "Means adjusted for temperature",
       x = "Species",
       y = "Adjusted Chirping Rate") +
  theme_minimal() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1))

Tip

🚀 If you finish early: Refit ancova_model with options(contrasts = c("contr.treatment", "contr.poly")) instead of contr.sum, and compare the coefficient table. Does the Anova() F-test for spp change? Does the intercept’s meaning change?

# Write your code here:

Part 2 · Partridge Longevity Analysis

We’ll analyze the effect of mating strategy on male fruitfly longevity, using thorax length as a covariate — this is real data from Partridge & Farquhar (1981).

▶ Run this — load and visualize the data:

p_df <- read.csv("data/partridge.csv") %>% clean_names() %>%
  rename(treat = treatmen)

p_df$treat <- factor(p_df$treat,
                      levels = 1:5,
                      labels = c("No females",
                                 "One virgin female daily",
                                 "Eight virgin females daily",
                                 "One inseminated female daily",
                                 "Eight inseminated females daily"))

head(p_df)

ggplot(p_df, aes(x = thorax, y = longev, color = treat)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Relationship between Thorax Length and Longevity",
       x = "Thorax Length (mm)",
       y = "Longevity (days)",
       color = "Treatment") +
  theme_minimal() +
  theme(legend.position = "bottom")

Step 1: Test homogeneity of slopes

homo_slopes_model <- lm(longev ~ thorax * treat, data = p_df)
Anova(homo_slopes_model, type = 3)

Step 2: Fit the ANCOVA model

p_ancova_model <- lm(longev ~ thorax + treat, data = p_df)
summary(p_ancova_model)

Anova(p_ancova_model, type = "II")

Step 3: Check assumptions

par(mfrow = c(2, 2))
plot(p_ancova_model)
par(mfrow = c(1, 1))

shapiro.test(p_ancova_model$residuals)

leveneTest( ~ thorax * treat, data = p_df)

# Breusch-Pagan (BP) test for homoscedasticity
lmtest::bptest(p_ancova_model)

Step 4: Calculate adjusted means

p_means <- emmeans(p_ancova_model, "treat")
p_means

Step 5: Pairwise comparisons

pairs(p_means, adjust = "tukey")

plot(p_means, comparisons = TRUE)

multcomp::cld(p_means)

✏️ Your turn: According to the compact letter display, which treatment(s) had the highest adjusted mean longevity? Does that match what you’d predict biologically? ________________________

Part 3 · Example with Heterogeneous Slopes

Let’s look at an example where slopes are NOT homogeneous, using simulated sea urchin data.

▶ Run this:

set.seed(345)
n <- 72  # 24 urchins per group

treatments <- rep(c("Initial", "Low Food", "High Food"), each = n/3)
volume <- c(
  runif(n/3, 10, 40),  # Initial
  runif(n/3, 10, 40),  # Low Food
  runif(n/3, 10, 40)   # High Food
)

suture_width <- ifelse(
  treatments == "Initial", 0.05 + 0.002 * volume,
  ifelse(
    treatments == "Low Food", 0.04 + 0.0005 * volume,
    0.02 + 0.003 * volume  # High Food
  )
) + rnorm(n, 0, 0.01)

u_df <- data.frame(treatment = treatments, volume = volume, suture_width = suture_width)

# Explicitly set "Initial" as the reference level for the factor
u_df$treatment <- factor(u_df$treatment, levels = c("Initial", "Low Food", "High Food"))

ggplot(u_df, aes(x = volume, y = suture_width, color = treatment)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(x = "Cube Root Body Volume",
       y = "Suture Width (mm)",
       color = "Treatment") +
  theme_minimal()

Test for homogeneity of slopes

urchin_model <- lm(suture_width ~ volume * treatment, data = u_df)
Anova(urchin_model, type = 3)

Result: With p < 0.05, we have heterogeneous slopes! Standard ANCOVA is inappropriate here.

What to do with heterogeneous slopes

When slopes are not homogeneous, you have several options.

Option 1 — Analyze groups separately:

initial_model <- lm(suture_width ~ volume, data = filter(u_df, treatment == "Initial"))
low_food_model <- lm(suture_width ~ volume, data = filter(u_df, treatment == "Low Food"))
high_food_model <- lm(suture_width ~ volume, data = filter(u_df, treatment == "High Food"))

initial_model
low_food_model
high_food_model

Option 2 — Johnson-Neyman procedure: identifies the specific regions of the covariate where groups do (and do not) significantly differ.

library(interactions)
jn_model <- lm(suture_width ~ volume + treatment + volume * treatment, data = u_df)

# pred = "volume" (the continuous moderator)
# modx = "treatment" (the categorical predictor)
sim_slopes(jn_model,
           pred = "volume",
           modx = "treatment",
           johnson_neyman = TRUE)

✏️ Your turn: Based on the Johnson-Neyman output, over what range of volume (if any) do the treatment groups differ significantly? ________________________

Tip

🚀 If you finish early: Try emtrends(urchin_model, "treatment", var = "volume") to get the estimated slope for each treatment group directly, instead of fitting three separate models.

# Write your code here:

Summary checklist for ANCOVA

When conducting ANCOVA, always follow these steps:

Tip

ANCOVA checklist

Visualize your data — plot response vs covariate, colored by group.
Test homogeneity of slopes — fit a model with an interaction term.
- If p > 0.05: proceed with ANCOVA.
- If p < 0.05: use alternative approaches (Johnson-Neyman, separate regressions).
Fit the ANCOVA model — response ~ covariate + factor.
Check assumptions — use diagnostic plots.
Interpret results — focus on adjusted means, not raw means.
Conduct post-hoc tests — pairwise comparisons if needed.
Visualize results — show adjusted means with confidence intervals.

Key points to remember:

ANCOVA increases power by accounting for covariate variation.
Adjusted means are what we compare, not raw group means.
Homogeneity of slopes is the most critical assumption.
Parallel lines in the plot suggest homogeneous slopes.
Non-parallel lines indicate heterogeneous slopes — use alternative methods.

Important

Key points from ANCOVA analysis

Test homogeneity of slopes first — the most critical assumption.
ANCOVA compares adjusted means at the mean value of the covariate.
Increases statistical power by removing variation due to the covariate.
Choose appropriate methods based on whether slopes are homogeneous.
Visualize your results showing the relationship between variables.
Check all assumptions using diagnostic plots.
Interpret in biological context — what do the adjusted means tell us?

Remember: the covariate should be measured independently of the treatment and should not be affected by the treatment itself!

Review and checkpoint

At this point you can:

Explain what ANCOVA is and when to use it
Test for homogeneity of slopes before running an ANCOVA
Fit an additive ANCOVA model and interpret its ANOVA table
Calculate and compare adjusted (estimated marginal) means
Recognize when slopes are heterogeneous and know the alternatives (Johnson-Neyman, separate regressions)

Note

📤 What to turn in before next class

Upload both of these to the course management system:

Your code — the scripts/ folder (or just 08_ancova.R)
This worksheet, with your written answers

Part 4 · Take-Home Extension — Real Cricket Chirp Data

Due Wednesday, November 25 — before the Multivariate Statistics class.

Part 1 used simulated cricket data to introduce ANCOVA. Below is real data on the same question — two cricket species, chirp rate vs. temperature — from Walker’s classic study. Same workflow, real noise, and this time you decide whether standard ANCOVA is even appropriate before you run it.

▶ Run this:

cricket_df <- read_csv("data/walker_cricket_chirp_rate.csv")
head(cricket_df)

ggplot(cricket_df, aes(x = temp, y = pulse, color = species)) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE)

Your Task

Question: Does chirp rate (pulse) differ between the two cricket species after accounting for temperature (temp)?

✏️ Your turn: State your hypotheses. Test the homogeneity-of-slopes assumption yourself and decide whether standard (additive) ANCOVA is appropriate, or whether you need to keep the interaction and report separate slopes instead. Justify your choice, fit the model, check the remaining assumptions, interpret it, and report the result properly.

H0 =
________________________________________________________
Ha =
________________________________________________________
Model I will use and why:
________________________________________________________

# Write your code here:

Interpretation:
________________________________________________________
________________________________________________________

Final Figure

Produce one publication-quality figure showing chirp rate vs. temperature by species (adjusted means with confidence intervals, or the raw regression lines — whichever fits your model from above) — proper axis labels, no default ggplot grey background — and export it with ggsave().

# Write your code here:

Note

📤 What to turn in — due Nov 25

Your write-up should include:

The question
Your hypotheses
Which model you used, and your justification for choosing it
The code and its output
Your interpretation of the result
A final figure showing the result

Note

Reference

Walker, T.J. (1962). Factors responsible for intraspecific variation in the calling songs of crickets. Evolution, 16, 407–428.

Getting unstuck

When code breaks — and it will, that is normal:

Read the error message out loud. R usually names the line and the problem.
Check the usual suspects: did you run all the library() calls at the top? Spelling? A missing ) or %>% at the start of a line?
?function_name opens the built-in help page.
Bring the exact error (copy-paste it) to class or office hours.

💡 Key idea: Every working scientist googles error messages daily. Getting stuck is not failing — it is the job.

Other Formats

Worksheet: ANCOVA — Analysis of Covariance

How to use this worksheet

Part 1 · Cricket Chirping Analysis

Step 1: Test homogeneity of slopes

Step 2: Fit the ANCOVA model

Step 3: Check model assumptions

Step 4: Calculate adjusted means

Step 5: Pairwise comparisons

Step 6: Visualize results

Part 2 · Partridge Longevity Analysis

Step 1: Test homogeneity of slopes

Step 2: Fit the ANCOVA model

Step 3: Check assumptions

Step 4: Calculate adjusted means

Step 5: Pairwise comparisons

Part 3 · Example with Heterogeneous Slopes

Test for homogeneity of slopes

What to do with heterogeneous slopes

Summary checklist for ANCOVA

Review and checkpoint

Part 4 · Take-Home Extension — Real Cricket Chirp Data

Your Task

Final Figure

Getting unstuck