Activity: Two-Way ANOVA
Factorial designs
In-class Activity 13: Two-Way (Factorial) ANOVA
Today’s Objectives
- Recognize a factorial design and state its three null hypotheses
- Say what balanced means, and what it costs to force it
- Fit a two-way factorial ANOVA with
lm()and read the interaction first - Check the assumptions on the residuals, not the raw data
- Compute estimated marginal means, and see why they are not just group means
- Compare Type I, II, and III sums of squares on a balanced design
- Run the same model on the unbalanced design and watch the three types come apart
- Report effect sizes with
effectsize, not just p-values - Report the result, and show the interaction as a figure
- Save the skeleton as
scripts/13_two_way_anova.Rin your project. - No data file. The penguin data comes from the
palmerpenguinspackage.
- 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, the data read, the summarise block, the balancing, and the section 7 EMM comparison — and nothing else. Every line of statistics you type yourself.
- 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 9 and 10 — is graded. Those answers go into the script as
#comments. - Hand in
scripts/13_two_way_anova.R. Run it top to bottom (Ctrl/Cmd + Shift + Enter) first — it must run with no errors.
Factorial ANOVA is easy to fake your way through by reading p-values off a table after the fact. Before each model, say out loud what you expect — which terms will be significant, whether the lines will be parallel — then run it. A wrong prediction points straight at the thing you misunderstood.
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 9, which is the graded deliverable.
🧩 Block 1 — The design and the data (after lecture Parts 1–2)
Sections 1–2: see the factorial structure, then build a balanced copy of it.
Section 1 · The data
Written for you in the skeleton — run it and read it.
# ---- Libraries -------------------------------------------
library(palmerpenguins)
library(car) # Anova(), leveneTest()
library(emmeans) # emmeans(), emmip()
library(effectsize) # eta_squared(), omega_squared()
library(broom) # augment()
library(tidyverse)# ---- 1 The data ------------------------------------------
p_df <- penguins %>%
select(spp = species, sex, mass_g = body_mass_g) %>%
drop_na(spp, sex, mass_g)
head(p_df)
p_df %>%
group_by(spp, sex) %>%
summarise(
n = n(),
mean_mass = mean(mass_g),
sd_mass = sd(mass_g),
se_mass = sd_mass / sqrt(n),
.groups = "drop"
)Two factors — spp with 3 levels, sex with 2 — crossed, so every species appears with both sexes. That is a 3 × 2 factorial design, and it has six cells.
✏️ Q1: Six groups, one response. A factorial ANOVA tests three separate null hypotheses, not one. State all three in words.
✏️ Q2: Look at the n column. Are the six cells the same size?
Section 2 · Count the cells, then force balance
Written for you in the skeleton — this is bookkeeping, not statistics.
# ---- 2 Count the cells, then force balance ---------------
p_df %>%
count(spp, sex) %>%
pivot_wider(names_from = sex, values_from = n)
min_n <- min(count(p_df, spp, sex)$n)
min_n
set.seed(123)
pb_df <- p_df %>%
group_by(spp, sex) %>%
slice_sample(n = min_n) %>%
ungroup()
pb_df %>%
count(spp, sex) %>%
pivot_wider(names_from = sex, values_from = n)slice_sample(n = min_n) inside group_by() takes a random min_n rows from each cell, so every cell ends up the same size. set.seed(123) makes that random draw repeatable — everyone in the room gets the same pb_df.
✏️ Q3: What was min_n? p_df has 333 birds; how many are in pb_df, and how many did we throw away to get balance?
✏️ Q4: Look carefully at the original counts. How many of the six cells were actually unbalanced? (Count them — the answer is smaller than you think.)
✏️ Q5: Throwing away a third of your data to get balance is a real cost. Why would anyone pay it on purpose? What does balance buy you?
💡 Key idea: we are about to run the same analysis twice — once on
pb_df(balanced) and once onp_df(unbalanced) — so that every difference you see between them is caused by the imbalance and nothing else. That is the whole design of today’s activity.
🧩 Block 2 — Fit it and check it (after lecture Parts 3–5)
Sections 3–4: fit the balanced factorial model, then check the assumptions.
Section 3 · Fit the balanced factorial model
One line fits it. We fit under R’s default contrasts for now — section 6 is where the contrast coding starts to matter.
🔮 Predict first: three terms — species, sex, and the interaction. Which will be significant? Will the interaction be?
▶ Type this yourself:
# ---- 3 Fit the balanced factorial model ------------------
pb_model <- lm(mass_g ~ spp * sex, data = pb_df)
summary(pb_model)
anova(pb_model)The formula mass_g ~ spp * sex is shorthand. R expands the * into three terms: spp + sex + spp:sex — the two main effects and the interaction.
✏️ Q6: Read the spp:sex row of the anova() table first, before the two main effects. Why does the interaction get read first? (What does a significant interaction do to the meaning of a main effect?)
✏️ Q7: Report the interaction: F, df, and p. Significant at α = 0.05?
✏️ Q8: summary() and anova() are both describing the same model. What does anova() give you that summary()’s coefficient table does not?
Section 4 · Assumptions
The four-panel plot is the one you have run every week. It is in your skeleton.
# ---- 4 Assumptions ---------------------------------------
par(mfrow = c(2, 2))
plot(pb_model)
par(mfrow = c(1, 1))▶ Type this yourself:
pb_resid_df <- augment(pb_model)
ggplot(pb_resid_df, aes(sample = .resid)) +
stat_qq() +
stat_qq_line()
ggplot(pb_resid_df, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed")▶ Type this yourself:
shapiro.test(residuals(pb_model))
leveneTest(mass_g ~ spp * sex, data = pb_df)✏️ Q9: Normality is an assumption about the residuals, not about the raw body masses. Why does that distinction matter here? (The three species have very different mean masses — what would a histogram of all 204 raw masses look like?)
✏️ Q10: leveneTest() is testing equal variance across how many groups here? Look at the df.
✏️ Q11: Both tests passed, or both failed? Either way — a factorial ANOVA on balanced data is fairly robust to unequal variance. Why would that robustness disappear once the design is unbalanced?
🧩 Block 3 — Marginal means and sums of squares (after lecture Parts 6–7)
Sections 5–6: get the estimated marginal means, then compare the three SS types on balanced data.
Section 5 · Estimated marginal means and post-hoc
An estimated marginal mean is the model’s estimate of a group mean, averaging over the other factor — and averaging the cell means equally, not the raw rows. The lecture’s fertilizer-and-light table is the whole idea in four numbers: when one cell is much bigger than the others, a raw group_by() mean lets that cell shout, and the EMM does not.
Last week, with one factor, the two were always identical. From here on they are not.
▶ Type this yourself:
# ---- 5 Estimated marginal means and post-hoc -------------
spp_emm <- emmeans(pb_model, ~ spp)
spp_emm
pairs(spp_emm)
sex_emm <- emmeans(pb_model, ~ sex)
sex_emm
pairs(sex_emm)
cell_emm <- emmeans(pb_model, ~ spp * sex)
cell_emm⚠️ Watch out!
emmeans()printsNOTE: Results may be misleading due to involvement in interactions. That is not an error — it is reminding you of Q6. With a significant interaction, the marginal mean for a species is an average across two sexes that behave differently.
🔮 Predict first: before you run the next chunk — on a balanced design, will the EMM for each cell equal the raw mean of that cell? Exactly, or just close?
▶ Type this yourself:
as.data.frame(cell_emm) %>%
select(spp, sex, emmean) %>%
pivot_wider(names_from = sex, values_from = emmean)
pb_df %>%
group_by(spp, sex) %>%
summarise(raw_mean = mean(mass_g), .groups = "drop") %>%
pivot_wider(names_from = sex, values_from = raw_mean)✏️ Q12: Compare the two tables. Do the EMMs match the raw cell means? Was your prediction right? Why does it come out that way on a balanced design?
▶ Type this yourself:
cld_df <- multcomp::cld(cell_emm, Letters = letters, adjust = "sidak") %>%
as.data.frame() %>%
arrange(spp, sex)
cld_df
emmip(pb_model, sex ~ spp, CIs = TRUE)emmip() is the one-line version of the interaction plot the lecture built by hand with ggplot(). Use it here to look; you will build the polished one in section 9.
✏️ Q13: In a compact letter display, cells sharing a letter are not distinguishable. Which cells share a letter?
✏️ Q14: emmip() draws one line per sex across the three species. If there were no interaction, what would those two lines look like? Do yours look like that?
✏️ Q15: Does the interaction plot tell you the same thing as the spp:sex p-value from Q7? Say how.
Section 6 · The three sums-of-squares types, balanced
Three different ways to decide what “the effect of species” means when there is more than one term in the model:
- Type I — sequential. Each term gets credit for whatever it explains after the terms listed before it. Order in the formula matters.
- Type II — each main effect conditions on the other main effect, but not on the interaction.
- Type III — each term conditions on everything else, interaction included.
🔮 Predict first: on this balanced dataset, will the three types give the same F-values, or different ones?
▶ Type this yourself:
# ---- 6 The three sums-of-squares types, balanced ---------
anova(pb_model) # Type I
Anova(pb_model, type = 2) # Type IIType I and Type II are contrast-free — they give the same answer however R coded the factors. Type III is not. It needs orthogonal contrasts (contr.sum), and the contrasts are baked into the model when it is fitted, so changing the option is only half the job: you have to refit.
▶ Type this yourself — set it, refit, then put R back the way you found it:
options(contrasts = c("contr.sum", "contr.poly"))
pb_model_3 <- lm(mass_g ~ spp * sex, data = pb_df)
Anova(pb_model_3, type = 3) # Type III
options(contrasts = c("contr.treatment", "contr.poly"))⚠️ Watch out! This is the single most common silent error in a factorial ANOVA. Run
Anova(pb_model, type = 3)on the default-contrast model and you get F = 196.1 forsppinstead of 477.2 — a real number, no error, no warning, and wrong. The reset line at the end matters too: leavecontr.sumswitched on and everysummary()coefficient table for the rest of your session reads against the grand mean instead of a reference level.
✏️ Q16: Write down the F-value for spp under all three types. What do you notice? Was your prediction right?
✏️ Q16b: Type III on pb_model_3 matches Types I and II. Type III on pb_model does not. Both are “Type III on a balanced design” — so what exactly is the options() line changing?
✏️ Q17: Explain why they come out that way. (What is it about equal cell sizes that makes “order of entry” stop mattering?)
💡 Key idea: this is the slide where the lecture computed \(SS_A\), \(SS_B\), \(SS_{AB}\) and \(SS_{res}\) by hand and they added up to \(SS_{total}\) exactly. That addition is what “the three types agree” means. On a perfectly balanced design the three SS types are the same number, so the choice is irrelevant and nobody has to defend it. Everything written about Type I vs. II vs. III is about what happens when that stops being true — which is section 7.
🧩 Block 4 — The payoff (after lecture Parts 8–10)
Sections 7–8: the same analysis, unbalanced, then the effect sizes that go with it. Section 7 is the one you do on your own.
The lecture warned you the penguin imbalance is mild — five of six cells are exactly even. So the three SS types will separate, but only a little, and no conclusion will change. Your job is to find where they separate and say why, not to find a dramatic result.
Section 7 · The same model, unbalanced (you write this one)
This is where we are headed all day. You have run four moves on pb_df: fit, check, EMMs, compare SS types. Now run the same four on the full p_df — all 333 birds, imbalance and all — and see what changes.
Only one of the six cells is unbalanced (Q4), so the honest expectation is that very little changes.
🔮 Predict first: with just one slightly-uneven cell out of six, will the three SS types still agree? Will the EMMs still match the raw means?
▶ Type this yourself — fit it, then the three SS types with the same contrast dance as section 6:
# ---- 7 The same model, unbalanced ------------------------
pu_model <- lm(mass_g ~ spp * sex, data = p_df)
summary(pu_model)
anova(pu_model) # Type I
Anova(pu_model, type = 2) # Type II
options(contrasts = c("contr.sum", "contr.poly"))
pu_model_3 <- lm(mass_g ~ spp * sex, data = p_df)
Anova(pu_model_3, type = 3) # Type III
options(contrasts = c("contr.treatment", "contr.poly"))▶ Type this yourself — the assumption checks:
shapiro.test(residuals(pu_model))
leveneTest(mass_g ~ spp * sex, data = p_df)The EMM comparison below is written for you in the skeleton — it is the same code as section 5 with p_df swapped in. Run it and read it.
pu_cell_emm <- emmeans(pu_model, ~ spp * sex)
as.data.frame(pu_cell_emm) %>%
select(spp, sex, emmean) %>%
pivot_wider(names_from = sex, values_from = emmean)
p_df %>%
group_by(spp, sex) %>%
summarise(raw_mean = mean(mass_g), .groups = "drop") %>%
pivot_wider(names_from = sex, values_from = raw_mean)
emmeans(pu_model, ~ spp)✏️ Q18: Fill in the F-values and compare the two runs (Type III from the _3 models):
balanced unbalanced
I II III I II III
spp
sex
spp:sex✏️ Q19: Which terms’ F-values moved apart across the three types, and which one did not move at all? Why is the interaction term immune?
✏️ Q20: Compare the marginal EMM for Gentoo to the raw Gentoo mean. They differ by about 10 g. Where does that 10 g come from? (Gentoo has 58 females and 61 males — and males are heavier.)
✏️ Q21: Which SS type would you report for this model, and why? Which would you report if the interaction had come out non-significant?
✏️ Q22: You had two versions of this analysis: one on 204 balanced birds, one on all 333. Which would you actually publish, and what did forcing balance buy you that was worth 129 birds?
💡 Key idea: imbalance is not a formatting problem, it is an interpretation problem. Once cells are uneven, “the effect of species” stops being one number and becomes a question about what you are holding constant. The three SS types are three different answers to that question. On balanced data they collapse into one, which is why balance is worth something — and why Type III is the safe default when you do not have it.
Section 8 · Effect size — how big, not just whether
A p-value answers “is there an effect?” It does not answer “how big?” With 333 birds, a trivial effect can still clear p < 0.05 — so every ANOVA you report needs an effect size next to it.
The lecture built partial η² by hand: \(\eta^2_p = SS_{effect} / (SS_{effect} + SS_{residual})\). You will not do that again. The effectsize package takes the Anova() object directly and returns all three terms with confidence intervals.
▶ Type this yourself — note it is pu_model_3, the orthogonal-contrast refit from section 7:
# ---- 8 Effect size ---------------------------------------
eta_squared(Anova(pu_model_3, type = 3))
omega_squared(Anova(pu_model_3, type = 3))effectsize prints a reminder about Type III needing orthogonal contrasts every time you call it — the same point as section 6. Here it is satisfied, because pu_model_3 was refit under contr.sum. Hand it pu_model instead and the warning is the only thing standing between you and a wrong effect size.
▶ Type this yourself — turn the numbers into words:
pu_eta <- eta_squared(Anova(pu_model_3, type = 3))
interpret_eta_squared(pu_eta$Eta2_partial, rules = "field2013")You should get roughly:
partial eta^2 95% CI label
spp 0.82 [0.80, 1.00] large
sex 0.49 [0.43, 1.00] large
spp:sex 0.05 [0.02, 1.00] smallSpecies accounts for about 82% of the variation it could have accounted for; sex about 49%; the interaction about 5%.
✏️ Q23: The spp:sex interaction was significant in section 7, and its partial η² is 0.05 — the smallest of the three. Say in one sentence what the effect size tells you that the p-value did not.
✏️ Q24: Compare eta_squared() to omega_squared(). They agree for spp but the interaction drops from 0.05 to 0.04. η² is biased upward in small samples and ω² corrects for it. Which term got corrected most, and why would the smallest effect be the one that moves?
✏️ Q25: The field2013 labels (“small / medium / large”) are cutoffs borrowed from psychology — 0.01, 0.06, 0.14. Is a partial η² of 0.05 for a sex-by-species interaction in penguin body mass “small” in any biological sense? Defend your answer without using the word “significant.”
💡 Key idea: report an effect size with every F-test, and report the same one throughout — partial η² pairs with Type III SS, which is why the lecture used it all day. The canned “small / medium / large” labels are a crutch for reading the scale, not a verdict. A 5% interaction in a body-mass model may be biologically real and worth a paragraph; a 5% effect somewhere else may be noise. The number is the finding, the label is not.
Homework
A two-factor design on rusty crayfish growth: native (Ohio) versus invasive (Wisconsin) populations, grown in enclosures in three northern Wisconsin lakes. Population origin × lake, with the interaction as the interesting term, plus the assumption checks and a publication-quality figure — and this time you decide which SS type to trust, and defend it.
Data: sargent_lodge_crayfish.csv — Sargent & Lodge (2014), Evolutionary Applications 7: 949–961. You meet these data again in Homework 9, where lake becomes a random effect.
Start a new script for it — do not add it to 13_two_way_anova.R.
Write-up — graded
These go into the bottom of your script, as # comments after ANSWER:.
Section 9 · Report the model
Use your pu_model_3 from section 7.
✏️ Q26: Fill in the Type III table:
effect F df p
spp
sex
spp:sex✏️ Q27: One sentence per term — species, sex, interaction — saying what each result means for penguins, not for the model.
✏️ Q28: Which SS type did you report, and why was it the right choice for an unbalanced design?
Section 10 · The final figure
One publication-quality interaction plot: body mass by species, with sex as separate lines or colors, error bars from the EMMs. Proper axis labels, a real theme, readable sizing — this one is graded on presentation. Build it from pu_cell_emm, not from the raw means, and export it with ggsave().
# ---- 10 The final figure ---------------------------------
pub_df <- as.data.frame(pu_cell_emm)
# ... your plot here, then:
# ggsave("figures/13_penguin_interaction.pdf", width = 7, height = 5)✏️ Q29: Your figure shows the interaction term as a picture. Are the lines parallel? Say how the figure and the spp:sex p-value agree.
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: real data is almost never balanced — something dies, a sample is lost, a site is unreachable. The point of today was not that imbalance is fatal; it is that imbalance makes you choose, and you should know which choice you made and why.