Activity: ANCOVA
Covariates
In-class Activity 17: ANCOVA
Today’s Objectives
- Say what a covariate is, and the two different jobs it does
- Test homogeneity of slopes, and know why it comes first
- Fit an additive ANCOVA and read its ANOVA table
- Measure what the covariate bought you — in residual SS and in F
- Read adjusted means, and say how they differ from raw means
- Recognise heterogeneous slopes and stop before fitting the wrong model
- Use
emtrends()to compare slopes, and contrast groups at covariate values - Say where an ANCOVA adjustment becomes extrapolation
📥 Script skeleton (.R) 📥 Fruitfly data 📥 Grayling data
- Save the skeleton as
scripts/17_ancova.R, and both.csvfiles in yourdata/folder. - You already have
gray_I3_I8.csvfrom week 5 — same fish, new question.
- All of today’s code is right here in this activity. Type it into your script as we go and run it chunk by chunk. The skeleton has the plumbing — libraries and the two data reads — and nothing else. Every line of statistics you type yourself.
- Section 6 has no code here either. By then you will have run the whole workflow once, so you start it again on your own.
- The Q boxes are things to think about, not inline homework. Read them, talk them over, answer them in your head. You do not have to type answers into the script, and they are not graded.
- The write-up at the bottom — sections 8, 9, and 10 — is graded. Those answers go into the script as
#comments. - Hand in
scripts/17_ancova.R. Run it top to bottom (Ctrl/Cmd + Shift + Enter) first — it must run with no errors.
Today has one decision in it, and everything downstream depends on getting it right: are the slopes parallel? You will meet one dataset where they are and one where they are not. Before each test, look at the plot and commit to an answer out loud.
invisible(NULL) lines in the skeleton
Every section ends with invisible(NULL). It is a bumper. Press Ctrl/Cmd + Enter on a section’s last real line and your cursor jumps to the next line that has code on it — without the bumper, that is the next section, way down past the next heading. With it, the cursor stops right there and you keep your place. It does nothing and prints nothing.
Keep the plots plain today — no labels, colors, or themes. The goal is to read them. The one exception is the final figure in section 10, which is the graded deliverable.
Where this fits
Week 10 changed the distribution. Week 11 changed the error structure. Today changes neither — it goes back to the ordinary linear model and asks what happens when you have one categorical predictor and one continuous predictor at the same time.
That is an ANCOVA, and it is the model that joins the ANOVA half of this course to the regression half.
| What it is | Example today | |
|---|---|---|
| Response | continuous | fruitfly longevity, days |
| Factor | categorical, the thing you care about | mating treatment |
| Covariate | continuous, the thing you adjust for | thorax length, mm |
A covariate does two different jobs, and they are worth keeping separate in your head:
| Job | What goes wrong without it | Fixed by |
|---|---|---|
| Confounding | groups sit at different covariate values, so a raw comparison of means confuses the two | comparing at a common covariate value |
| Noise | the covariate’s variation sits in the error term, shrinking your F | removing it from the residual |
The one rule. The covariate must be measured independently of the treatment and must not be affected by it. Adjusting for something your treatment caused deletes the effect you are trying to measure. No diagnostic will catch this — only knowing how the data were collected.
🧩 Block 1 — The data and the critical test (after lecture Parts 1–3)
Sections 1–2: meet the fruitflies, then run the one test that decides everything else.
Section 1 · The fruitfly data
Written for you in the skeleton — run it and read it.
# ---- Libraries -------------------------------------------
library(janitor)
library(car) # Anova()
library(emmeans) # emmeans(), emtrends()
library(broom)
library(lmtest) # bptest()
library(interactions) # sim_slopes()
library(tidyverse)
options(scipen = 999)
options(contrasts = c("contr.sum", "contr.poly")) # for Anova(type = 3)Partridge & Farquhar (1981) asked whether sexual activity shortens a male fruitfly’s life. They housed 125 males, 25 in each of five treatments, and recorded how long each lived.
The design is cleverer than it first looks:
| Treatment | Females | Do they mate? |
|---|---|---|
| No females | 0 | — |
| One inseminated female daily | 1 | no — newly inseminated females will not re-mate |
| Eight inseminated females daily | 8 | no |
| One virgin female daily | 1 | yes |
| Eight virgin females daily | 8 | yes |
# ---- 1 The fruitfly data ---------------------------------
p_df <- read_csv("data/partridge.csv", show_col_types = FALSE) %>%
clean_names() %>%
mutate(
treat = factor(
case_when(
partners == 0 ~ "no females",
partners == 1 & type == 0 ~ "1 inseminated",
partners == 8 & type == 0 ~ "8 inseminated",
partners == 1 & type == 1 ~ "1 virgin",
partners == 8 & type == 1 ~ "8 virgin"
),
levels = c("no females", "1 inseminated", "8 inseminated",
"1 virgin", "8 virgin")
)
)
head(p_df)
count(p_df, treat)⚠️ Why the labels are built from partners and type, not from treatmen
The file has a treatmen column coded 1–5, but the codes are not in the order you would guess: 1 is eight inseminated, 2 is no females, 3 is one inseminated, 4 is one virgin, 5 is eight virgin.
Labelling them 1→“none”, 2→“one virgin”, 3→“eight virgin” — the obvious guess — reverses the biological conclusion. Building the labels from the design columns (partners and type) instead means they cannot drift out of step with the data.
✏️ Q1: What is the response, what is the factor, and what is the covariate?
✏️ Q2: Why are the inseminated-female groups in this experiment at all? What would you be unable to conclude without them?
✏️ Q3: Thorax length is the covariate. Could the mating treatment have changed a fly’s thorax length? Why does the answer matter?
Section 2 · Look first, then test the slopes
▶ Type this yourself:
# ---- 2 Look first, then test the slopes ------------------
ggplot(p_df, aes(x = thorax, y = longev, color = treat)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE)🔮 Predict first: look at those five lines. Are they parallel? Commit to an answer before you run the test.
▶ Type this yourself — the homogeneity-of-slopes test:
slopes_model <- lm(longev ~ thorax * treat, data = p_df)
Anova(slopes_model, type = 3)✏️ Q4: Write down the F and p for the thorax:treat interaction. Are the slopes homogeneous?
✏️ Q5: This is a test you want to fail to reject. Name one other test you have run this semester where that was also true.
✏️ Q6: Suppose it had come back p = 0.001. What would you do next, and what would you not be allowed to report?
✏️ Q7: Look at the treat row in this table. Is it significant? Hold onto that — you are about to fit a different model and get a very different answer.
💡 Key idea: the interaction term asks whether the groups share one slope. If they do not, there is no single “difference between the groups” to report — the difference depends on where you are along the covariate. That is why this test comes first, before the model you intend to publish.
🧩 Block 2 — The additive model (after lecture Part 4)
Sections 3–5: fit it, measure what the covariate bought, then read the adjusted means.
Section 3 · Fit the ANCOVA
Slopes are homogeneous, so drop the interaction.
▶ Type this yourself:
# ---- 3 Fit the ANCOVA ------------------------------------
ancova_model <- lm(longev ~ thorax + treat, data = p_df)
summary(ancova_model)
Anova(ancova_model, type = 2)✏️ Q8: Both thorax and treat are significant now. In Q7 the treat row was not. Same data, same question — what changed?
✏️ Q9: The thorax coefficient is about 136. What are its units, and what does it mean in fruitflies?
✏️ Q10: summary() shows treat1, treat2, treat3, treat4 — four coefficients for five groups. Why four, and what is the fifth group’s effect?
▶ Type this yourself — check the assumptions on the final model:
par(mfrow = c(2, 2))
plot(ancova_model)
par(mfrow = c(1, 1))
shapiro.test(residuals(ancova_model))
bptest(ancova_model)✏️ Q11: Do the diagnostics look acceptable? Which plot tells you about equal variances?
Section 4 · What did the covariate buy you?
This is the whole power argument, in two models.
🔮 Predict first: run the same question with no covariate at all. Will the treatment effect get stronger, weaker, or stay the same?
▶ Type this yourself:
# ---- 4 What did the covariate buy you? -------------------
no_cov_model <- lm(longev ~ treat, data = p_df)
Anova(no_cov_model, type = 2)
Anova(ancova_model, type = 2)✏️ Q12: Fill this in from the two tables:
residual SS residual df treat F
no covariate
ANCOVA✏️ Q13: What fraction of the residual sum of squares did one degree of freedom of thorax remove?
✏️ Q14: The treatment effect was significant both ways. So was the covariate worth including? Give a reason that is not about this particular p-value.
✏️ Q15: The five groups have almost identical mean thorax lengths (check with group_by(treat) %>% summarise(mean(thorax))). Given that, which of the covariate’s two jobs — confounding or noise — is it doing here?
💡 Key idea: a covariate earns its place even when your groups are perfectly balanced on it. Balance kills the confounding problem; it does nothing about the noise problem. Measure the covariate anyway.
Section 5 · Adjusted means
▶ Type this yourself:
# ---- 5 Adjusted means ------------------------------------
p_emm <- emmeans(ancova_model, "treat")
p_emm
p_df %>%
group_by(treat) %>%
summarise(raw_mean = mean(longev), mean_thorax = mean(thorax))▶ Type this yourself — the comparisons:
pairs(p_emm, adjust = "tukey")
multcomp::cld(p_emm, Letters = letters)
plot(p_emm, comparisons = TRUE)✏️ Q16: Put the adjusted means next to the raw means. How far apart are they, and why is the gap small here?
✏️ Q17: What covariate value are these means adjusted to? (Check mean(p_df$thorax).)
✏️ Q18: Read the letter display. Which groups form one block, and what do those groups have in common biologically?
✏️ Q19: Say the result in one sentence about fruitflies — not about F and p. What exactly is costing them their lifespan?
✏️ Q20: Four of the ten pairwise comparisons are not significant. Which four, and does that pattern support or undercut the conclusion in Q19?
💡 Key idea: adjusted means are what ANCOVA compares — every group’s predicted response at the same covariate value. When the groups are already balanced on the covariate, they land close to the raw means and the covariate is pure power. When the groups are not balanced, the two can differ a lot, and the raw means are the misleading ones.
🧩 Block 3 — When the slopes are not parallel (after lecture Part 5)
Sections 6–7: the same workflow on data that refuses to cooperate.
Section 6 · Grayling allometry (you start this one)
No code for part of this page. You have just run the first two steps of the workflow once. Do them again.
Back to the Arctic grayling from week 5 — lakes I3 and I8. Fish mass scales with length as a power law, so on log-log axes the relationship is a straight line and its slope is the allometric exponent: mass ∝ lengthb. A slope of 3 means the fish keeps its shape as it grows.
Question: do grayling in the two lakes have the same length–mass relationship?
Written for you in the skeleton:
# ---- 6 Grayling allometry --------------------------------
g_df <- read_csv("data/gray_I3_I8.csv", show_col_types = FALSE) %>%
clean_names() %>%
filter(mass_g > 0) %>%
mutate(log_len = log10(length_mm),
log_mass = log10(mass_g))
count(g_df, lake)Now you:
- Plot
log_massagainstlog_len, colored bylake, with a line per lake. - Run the homogeneity-of-slopes test — the interaction model, then
Anova(type = 3). - Decide whether you are allowed to fit
log_mass ~ log_len + lake.
💡 Hint: the only things that change from sections 2–3 are the data frame and the three variable names.
✏️ Q21: What is the F and p for the interaction? What does that forbid you from doing?
✏️ Q22: Fit lm(log_mass ~ log_len + lake) anyway and look at the lake term. It is wildly significant. Why is that number not reportable?
✏️ Q23: Compare the two lakes’ length ranges (group_by(lake) %>% summarise(min(length_mm), max(length_mm))). What does that tell you about where you are allowed to compare them?
Section 7 · Reporting heterogeneous slopes
When slopes differ, the slopes are the result.
▶ Type this yourself:
# ---- 7 Reporting heterogeneous slopes --------------------
g_int <- lm(log_mass ~ log_len * lake, data = g_df)
emtrends(g_int, "lake", var = "log_len")
pairs(emtrends(g_int, "lake", var = "log_len"))▶ Type this yourself — is each slope different from 3, the isometric value?
test(emtrends(g_int, "lake", var = "log_len"), null = 3)▶ Type this yourself — where along the covariate do the lakes actually differ?
lens <- log10(c(200, 250, 300))
pairs(emmeans(g_int, ~ lake | log_len, at = list(log_len = lens)))✏️ Q24: Write down each lake’s slope with its SE. Which one is indistinguishable from 3, and what does isometry mean for a fish?
✏️ Q25: The other slope is below 3. In words: what happens to those fish as they get longer?
✏️ Q26: Look at the contrasts at 200, 250 and 300 mm. At which length are the lakes indistinguishable, and what happens to the gap as fish get bigger?
✏️ Q27: You might have reached for sim_slopes(..., johnson_neyman = TRUE) here. Try it. What does it tell you, and why? (Read the warning carefully.)
✏️ Q28: An additive ANCOVA would have reported one “lake effect.” Name the biological finding it would have thrown away.
💡 Key idea: a significant interaction is not a failed analysis — it is usually the more interesting one. “These groups differ by a constant amount” is a weaker, less informative claim than “these groups differ in how they scale.” When the data offer you the second, do not flatten it into the first.
Write-up — graded
These go into the bottom of your script, as # comments after ANSWER:.
Section 8 · Report the fruitfly model
✏️ Q29: Fill in:
homogeneity of slopes F( , ) = p =
covariate (thorax) F( , ) = p =
treatment F( , ) = p =
adjusted mean, 8 virgin ±
adjusted mean, no females ±✏️ Q30: A 3–4 sentence results paragraph. Report the slopes test, then the ANCOVA, then the adjusted means with the letter display. Say what value the means are adjusted to.
✏️ Q31: One sentence on the biology: what does this experiment show costs a male fruitfly its lifespan, and which treatment groups let you rule out the alternative?
Section 9 · Report the grayling model
✏️ Q32: The interaction test, both slopes with their SEs, and the contrast at 250 mm.
✏️ Q33: Two or three sentences reporting it. You are reporting an interaction, so report slopes — not an adjusted mean.
✏️ Q34: Why would it have been wrong to report a single “lake effect” for these fish? Answer in a way a reviewer would accept.
Section 10 · The final figure
One publication-quality figure. Pick either the fruitfly adjusted means with their intervals, or the grayling log-log regression lines by lake with the raw points. Proper axis labels, a real theme, readable sizing — this one is graded on presentation. Export it with ggsave().
# ---- 10 The final figure ---------------------------------
# ... your plot here, then:
# ggsave("figures/17_ancova.pdf", width = 7, height = 5)✏️ Q35: If you drew the grayling figure, why must the two lines stop where the data stop? If you drew the fruitfly figure, why does the caption have to name the thorax length the means were adjusted to?
Review and checkpoint
At this point you should be able to:
Run your whole script top to bottom with Ctrl/Cmd + Shift + Enter. It must run clean.
💡 Key idea: every model this semester has been the same linear model with different pieces swapped in. ANCOVA swaps in a continuous predictor alongside a categorical one — and the only genuinely new question it asks is whether those two predictors interact.
Homework
Chirp rate against temperature in two tree cricket species — the whole workflow on your own: test the interaction first to see whether the slopes differ, then interpret the species effect at a common temperature with emmeans(), and check the assumptions before you trust any of it.
Data: cricket_chirp_rate.csv — Walker (1962), Evolution 16: 407–428.
Start a new script for it — do not add it to 17_ancova.R.