Lecture: Analysis of Variance (ANOVA)

Single-factor ANOVA

Bill Perry

Where We Left Off

Covered last time — Multiple Regression:

• Partial slopes, and why they are conditional
• Collinearity, and saying “these data cannot separate them”
• Interactions
• Factors, dummy variables, and fct_relevel()
• Model comparison with adjusted R² and AIC

Note

✅ Key idea from last lecture

You already built an ANOVA table for a regression, and you already turned a categorical predictor into dummy variables.

Today every predictor is categorical. The machinery does not change at all — only the name on the front does.

Part 1 · Overview

Today’s Objectives

One-way ANOVA, end to end:

• Why ANOVA is regression with a categorical predictor
• Partitioning variance, and the F-ratio
• Type I, II and III sums of squares, and the contrast setting Type III needs
• Assumptions and diagnostics
• Setting the reference level on purpose
• Post-hoc tests with emmeans — which pairs actually differ, and the letters that report it
• Reporting results
• What to do when the assumptions genuinely fail

📥 Crab data 📥 Nematode data

Note

📖 Reading

Whitlock & Schluter, Ch. 15 — Comparing Means of More Than Two Groups.

Gotelli & Ellison, Ch. 10 — The Analysis of Variance.

	Continuous X	Categorical X
Continuous Y	Regression	ANOVA
Categorical Y	Logistic	Contingency

The Study · Does a Big Claw Act as a Radiator?

A male fiddler crab has one hugely enlarged claw. Darnell & Munguia asked if it acts as a heat radiator, so the crab can stay out foraging in the sun. If the big claw sheds heat, removing it should make males hotter.

• Four groups of 21 crabs each, body temperature above ambient measured in the field:

  • intact male — the control
  • male major removed — the big claw taken off
  • male minor removed — a small claw taken off, the procedural control
  • female — never had a major claw at all

crab_df <- read_csv("data/fiddler_crabs.csv") %>%
  mutate(crab_type = factor(crab_type))

count(crab_df, crab_type)

# A tibble: 4 × 2
  crab_type              n
  <fct>              <int>
1 female                21
2 intact male           21
3 male major removed    21
4 male minor removed    21

Note

A balanced design

Equal n in every group — which makes the degrees-of-freedom arithmetic below clean, and the test robust.

🗳️ How Today Works

Warning

🗳️ Six times today, you make the call

Six slides are marked 🗳️ DECIDE. Each asks a multiple-choice question, and I will not tell you the answer on that slide.

• How it works:

  • write your letter down — actually write it
  • then hold up fingers 1–4
  • If the room splits, turn to the person next to you, argue for ninety seconds, and we vote again

Part 2 · ANOVA Is Regression

Same Machinery, Different Predictor

Regression

\[Y_i = \beta_0 + \beta_1X_i + \varepsilon_i\]

• X is continuous
• H₀: β₁ = 0
• One slope

ANOVA

\[Y_{ij} = \mu + A_i + \varepsilon_{ij}\]

• \(A_i\) is the effect of group i
• H₀: μ₁ = μ₂ = … = μₖ
• One coefficient per group, minus one

Both:

• Partition total variation into explained and left over
• Divide each sum of squares by its degrees of freedom
• Test with an F-ratio of the two mean squares

Important

✅ The same lm() call

# regression
lm(absorbance ~ nitrate, data = df)

# ANOVA
lm(body_temp_c ~ crab_type, data = df)

The only difference is that crab_type is a factor, so R builds dummy variables for you — exactly as it did last week.

What R Actually Fits

lm() cannot multiply by the word “female.” Behind the formula it builds a dummy variable for every level except the first — a column of 0s and 1s saying “is this crab in that group?”

crab_model <- lm(body_temp_c ~ crab_type, data = crab_df)

model.matrix(crab_model) %>% head(8)

  (Intercept) crab_typeintact male crab_typemale major removed
1           1                    0                           0
2           1                    0                           0
3           1                    0                           0
4           1                    0                           0
5           1                    0                           0
6           1                    0                           0
7           1                    0                           0
8           1                    0                           0
  crab_typemale minor removed
1                           0
2                           0
3                           0
4                           0
5                           0
6                           0
7                           0
8                           0

Four groups, three dummy columns. The level with no column of its own is the reference, carried by the intercept.

Note

✅ This is the whole trick

A one-way ANOVA is a multiple regression whose predictors happen to be 0/1 columns — the dummy variables you met last week.

aov() and lm() fit the identical model; aov() just prints it as an ANOVA table.

The Coefficients Are the Group Means

coefs <- coef(crab_model)

crab_df %>%
  group_by(crab_type) %>%
  summarise(group_mean = mean(body_temp_c)) %>%
  mutate(
    from_coefficients = c(
      coefs[1],
      coefs[1] + coefs[2],
      coefs[1] + coefs[3],
      coefs[1] + coefs[4]
    )
  )

# A tibble: 4 × 3
  crab_type          group_mean from_coefficients
  <fct>                   <dbl>             <dbl>
1 female                   1.68              1.68
2 intact male              1.29              1.29
3 male major removed       1.41              1.41
4 male minor removed       1.21              1.21

The intercept is the mean of the reference group. Every other coefficient is that group’s difference from it, so adding them back recovers each group mean exactly.

Important

✅ Why this matters in a moment

If the coefficients are “differences from the reference,” then which level is the reference decides what every coefficient means.

R chose it alphabetically. Part 7 is about taking that choice back.

🗳️ DECIDE · Why Not Just Run All the t-tests?

With four crab groups there are six possible pairwise t-tests. Suppose every group really has the same mean, and you run all six at α = 0.05.

What is the chance you get at least one “significant” result?

(a) 5% — that is what α = 0.05 means

(b) About 12%

(c) About 26%

(d) 30% — 6 × 5%

Warning

✍️ Write your letter down, then vote.

Note

The question to ask

What is the chance of getting no false positives in six independent tries?

The Multiple-Comparisons Problem

Important

✅ The answer is (c) — about 26%.

\[1 - (1 - 0.05)^6 = 0.265\]

Each test has a 95% chance of behaving. Six in a row behave only 73.5% of the time, so you have roughly a one in four chance of at least one false positive — before you have learned anything.

(d) is the common guess. Adding probabilities overshoots, because the ways of being wrong overlap.

ANOVA’s job is to ask one question first: is there any difference at all among these groups?

One test, one α, one p-value. Only if it rejects do you go looking for which groups differ — and then you use a post-hoc method that controls the error rate across that whole family of comparisons.

With 10 groups there are 45 pairs, and the chance of at least one false positive is over 90%.

Part 3 · The Logic

ANOVA’s Question

Even if every group had the same true mean, the sample means would differ, just from sampling error.

ANOVA asks whether the sample means differ more than sampling error alone would produce:

• \(MS_{among}\) — how spread out the group means are
• \(MS_{within}\) — how spread out individuals are inside their own group

\[F = \frac{MS_{among}}{MS_{within}}\]

• H₀: μ₁ = μ₂ = μ₃ = μ₄
• H₁: at least one mean differs

Where the Variation Goes

Every point’s distance from the grand mean splits into two pieces, exactly as in regression:

\[\underbrace{(Y_{ij} - \bar{Y})}_{\text{total}} = \underbrace{(\bar{Y}_i - \bar{Y})}_{\text{among}} + \underbrace{(Y_{ij} - \bar{Y}_i)}_{\text{within}}\]

Square and sum over every observation:

\[SS_{total} = SS_{among} + SS_{within}\]

The dashed line is the grand mean; the coloured bars are the group means.

Part 4 · Partitioning the Sum of Squares

The Sums of Squares

\[SS_{total} = \sum_{i}\sum_{j}(Y_{ij} - \bar{Y})^2\]

\[SS_{among} = \sum_{i} n_i(\bar{Y}_i - \bar{Y})^2\]

\[SS_{within} = \sum_{i}\sum_{j}(Y_{ij} - \bar{Y}_i)^2\]

with a groups and \(n_i\) observations in group i:

• \(df_{among} = a - 1\) = 3
• \(df_{within} = N - a\) = 80
• \(df_{total} = N - 1\) = 83

Tip

🖐 Note the \(n_i\)

Each group mean’s deviation is weighted by that group’s own sample size. With a balanced design every \(n_i\) is the same and it simplifies to \(n\sum(\bar{Y}_i - \bar{Y})^2\) — but only then.

Our design is balanced at n = 21, so both forms agree here.

Computing Them by Hand

grand_mean <- mean(crab_df$body_temp_c)

crab_parts <- crab_df %>%
  group_by(crab_type) %>%
  mutate(group_mean = mean(body_temp_c)) %>%
  ungroup()

ss_total <- sum((crab_parts$body_temp_c - grand_mean)^2)
ss_among <- sum((crab_parts$group_mean - grand_mean)^2)
ss_within <- sum((crab_parts$body_temp_c - crab_parts$group_mean)^2)

c(total = ss_total, among = ss_among, within = ss_within)

   total    among   within 
6.108929 2.641310 3.467619 

ss_among + ss_within

[1] 6.108929

2.641 + 3.468 = 6.109 — exactly, as it must.

Important

✅ The same identity as regression

• \(SS_{among}\) here is \(SS_{regression}\) there
• \(SS_{within}\) here is \(SS_{residual}\) there

Different names, same partition, same arithmetic.

Look at summary() First

crab_model <- lm(body_temp_c ~ crab_type, data = crab_df)
summary(crab_model)

Call:
lm(formula = body_temp_c ~ crab_type, data = crab_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.57619 -0.11071  0.01429  0.12381  0.61429 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  1.67619    0.04543  36.895  < 2e-16 ***
crab_typeintact male        -0.39048    0.06425  -6.077 3.93e-08 ***
crab_typemale major removed -0.26190    0.06425  -4.076 0.000107 ***
crab_typemale minor removed -0.46667    0.06425  -7.263 2.20e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2082 on 80 degrees of freedom
Multiple R-squared:  0.4324,    Adjusted R-squared:  0.4111 
F-statistic: 20.31 on 3 and 80 DF,  p-value: 6.997e-10

Same model, two different tables. summary() gives you the coefficients; anova() gives you the partition.

• (Intercept) = 1.676 — the mean of the reference level, which R picked alphabetically: female
• Each slope — that level’s mean minus the reference, not its own mean
• Bottom line — the same overall F = 20.31 on the same two df that anova() is about to print

Note

🖐 Where to look for what

The intercept only ever appears in summary() and in a Type III table. anova() and Anova(type = "II") do not print an (Intercept) row at all.

The ANOVA Table from R

anova(crab_model)

Analysis of Variance Table

Response: body_temp_c
          Df Sum Sq Mean Sq F value    Pr(>F)    
crab_type  3 2.6413 0.88044  20.312 6.997e-10 ***
Residuals 80 3.4676 0.04335                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Every number is one you just built:

• Sum Sq — 2.64 and 3.47
• Df — a − 1 = 3, N − a = 80
• Mean Sq — Sum Sq ÷ Df
• F value — the ratio of those two mean squares

F = 20.31, p = 7.00e-10.

🗳️ DECIDE · What Does That F Tell You?

The ANOVA returns F = 20.31, p = 7.00e-10.

What have you just learned?

1. All four crab types have different mean temperatures
2. At least one crab type differs from at least one other
3. Removing the major claw raises body temperature
4. The major claw explains 43.2% of the variation in temperature, so the radiator hypothesis is supported

Warning

✍️ One letter. Write it down.

Note

The question to ask

• What was the null hypothesis, exactly?
• A test can only tell you about the hypothesis it tested.

What F Does and Does Not Say

Important

✅ The answer is (b).

H₀ was “all four means are equal.” Rejecting it says at least one differs — not which, not how many, and not in which direction.

(a) claims far more than the test supports. (c) is a directional claim about one specific pair; that needs a post-hoc comparison, and Part 8 shows it is not even true. (d) misreads R²: 0.432 is the variance explained by crab type, which bundles sex and two surgeries together, not by the claw alone.

Warning

⚠️ The most common ANOVA mistake

Reporting a significant F and then describing whichever group means look farthest apart, as though the F test had said so.

It did not. The F test is an entry ticket, not an answer.

Which Sums of Squares Did You Just Use?

anova() gives Type I — sequential sums of squares. It adds terms to the model in the order you wrote them and credits each one with whatever variation is left over at that point.

With three types in circulation:

• Type I (sequential) — each term adjusted for the terms before it. Order matters
• Type II (hierarchical) — each term adjusted for all other terms at the same level, ignoring interactions
• Type III (marginal) — each term adjusted for every other term, interactions included

anova() is Type I only. Types II and III come from car::Anova() — capital A.

Anova(crab_model, type = "II")

Anova Table (Type II tests)

Response: body_temp_c
          Sum Sq Df F value    Pr(>F)    
crab_type 2.6413  3  20.312 6.997e-10 ***
Residuals 3.4676 80                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sum Sq 2.6413, F = 20.31 — identical to the anova() table a few slides back.

Important

✅ In a one-way ANOVA all three types agree

One predictor, so there is nothing to adjust it for. Balanced or unbalanced, Type I = Type II = Type III here — imbalance alone is never enough; it takes two predictors sharing variation.

So this is the safe place to learn the machinery, before the choice changes your answer.

Type III Needs a Different Contrast Setting

Type III tests each term at the average of the others, and that is only well defined if the factor’s coding is centred on the grand mean. R’s default contr.treatment is not — it is centred on the reference level.

So a true Type III needs the option switched, and the model refit:

options(contrasts = c("contr.sum", "contr.poly"))

crab_sum <- lm(body_temp_c ~ crab_type, data = crab_df)
Anova(crab_sum, type = "III")

Anova Table (Type III tests)

Response: body_temp_c
             Sum Sq Df  F value    Pr(>F)    
(Intercept) 163.801  1 3778.986 < 2.2e-16 ***
crab_type     2.641  3   20.312 6.997e-10 ***
Residuals     3.468 80                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

options(contrasts = c("contr.treatment", "contr.poly"))

The crab_type row is still 2.6413 and F = 20.31. Only the (Intercept) row moved — it now tests the grand mean instead of the reference group’s mean. The next slide shows that row in summary().

Warning

⚠️ Contrasts are baked in at fit time

options() changes nothing about a model you already fitted. Change the option then refit — otherwise Anova(..., type = "III") silently reports Type III sums of squares computed from treatment coding, which are not the Type III anyone means.

Note

📖 What contr.sum actually changes

Nothing about the data or the fit — only how the factor’s dummy columns are built. contr.treatment codes each level against the reference level; contr.sum codes it against the grand mean.

Same model, same residuals, same F. Different bookkeeping.

It Is the Coding That Moves the Intercept

anova() never showed you the intercept, so run summary() on both fits and watch the one row that changed:

# contr.treatment — R's default
round(summary(crab_model)$coefficients[, 1:2], 3)

                            Estimate Std. Error
(Intercept)                    1.676      0.045
crab_typeintact male          -0.390      0.064
crab_typemale major removed   -0.262      0.064
crab_typemale minor removed   -0.467      0.064

# contr.sum — the model refit for Type III
round(summary(crab_sum)$coefficients[, 1:2], 3)

            Estimate Std. Error
(Intercept)    1.396      0.023
crab_type1     0.280      0.039
crab_type2    -0.111      0.039
crab_type3     0.018      0.039

1.676 is the reference group’s mean (female). 1.396 is the grand mean — and the slopes are now deviations from it, crab_type1…crab_type3.

Important

✅ Type III does not move the intercept — contr.sum does

They just travel together, because Type III is only well defined under grand-mean coding. Ask a treatment-coded model for Type III and its (Intercept) row still tests the reference mean: SS = 59 against 163.8 under contr.sum.

The crab_type row is 2.641 either way.

options() Is Global · Always Set It Back

options(contrasts = c("contr.sum", "contr.poly"))

# ... your Type III model and Anova() call ...

options(contrasts = c("contr.treatment", "contr.poly"))

options() is global and sticky for the whole session. It is not attached to a model, a data frame or a script — it stays changed until you change it back or restart R.

Put the two calls in the same block, with the Type III work sandwiched between them. Not at the top of the script and “later.” Not in a different chunk.

Warning

⚠️ What going home with contr.sum left on looks like

Every summary() you run afterwards prints

crab_type1   crab_type2   crab_type3

— deviations from the grand mean — instead of named contrasts against your control. emmeans and TukeyHSD still work, but your coefficient table is answering a different question.

Nothing errors. The output just quietly stops meaning what you think it means, and that is far worse than a crash.

Two-way ANOVA is where this starts to bite — that is where Type I and Type III give genuinely different answers, and where you will set this option on purpose.

Type I in Detail · Order Is the Whole Story

The slide before said all three types agree in a one-way ANOVA. That still holds — it had one predictor. Add a second one and imbalance finally has something to act on, so here is next week’s penguin data: body mass ~ species + sex, 333 birds, and the species–sex counts are not quite equal.

Type I walks the formula left to right. The first term is tested against nothing; the second is tested against whatever variation the first one left behind.

anova(lm(body_mass_g ~ species + sex, peng_df))["species", "Sum Sq"]

[1] 145190219

anova(lm(body_mass_g ~ sex + species, peng_df))["species", "Sum Sq"]

[1] 143401584

Same data. Same two predictors. Same model. Species is worth 145.2 million when it goes first and 143.4 million when it goes second.

Type I sums of squares do add exactly to \(SS_{total}\) — the one thing they have going for them, and why anova() still defaults to them.

Note

📖 Where the difference comes from

Species and sex are slightly correlated here — the species are not sampled in identical sex ratios. That overlap is 1.8 million of variation that either predictor could claim.

Type I gives all of it to whichever one you typed first. The second term gets the leftovers.

Warning

⚠️ This is a real decision, made by accident

Nobody chooses “species should get credit for the overlap.” They just type species + sex because species felt like the main thing.

Type I turns the order of your typing into a statistical claim.

Why Type II or III · What Type II Does

Type II asks a cleaner question: what is each term worth after adjusting for all the others? No order, because every term is adjusted for every other one.

Anova(lm(body_mass_g ~ species + sex, peng_df), type = "II")

Anova Table (Type II tests)

Response: body_mass_g
             Sum Sq  Df F value    Pr(>F)    
species   143401584   2  715.29 < 2.2e-16 ***
sex        37090262   1  370.01 < 2.2e-16 ***
Residuals  32979185 329                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Species gets 143.4 million — exactly its “went last” Type I value. That is the trick to remembering it:

Type II gives every term the sequential SS it would get if it went last.

Nobody pockets the overlap, so Type II sums of squares do not add to \(SS_{total}\). The 1.8 million of shared variation is simply not credited to anyone, which is the honest answer.

Tip

🖐 Which one to use

• Type II — your default when there is no interaction, or the interaction is not significant. It is the most powerful of the three, because it does not hold back the variation an absent interaction would have claimed
• Type III — when you are keeping an interaction. It goes one step further than Type II and adjusts each main effect for the interaction too, which is why it needs contr.sum to mean anything
• Type I — only when the order is your hypothesis: fitting a covariate first to strip it out before testing your treatment

Balance Makes the Whole Problem Vanish

Take the same penguins, force 34 per species-by-sex cell, and run both orders again:

set.seed(12)
peng_bal <- peng_df %>%
  group_by(species, sex) %>%
  slice_sample(n = 34) %>%
  ungroup()

anova(lm(body_mass_g ~ species + sex, peng_bal))["species", "Sum Sq"]

[1] 95895631

anova(lm(body_mass_g ~ sex + species, peng_bal))["species", "Sum Sq"]

[1] 95895631

Identical, to the last digit. With equal n the two predictors are perfectly uncorrelated, there is no shared variation to fight over, and Type I = Type II = Type III.

Important

✅ Why this matters to us specifically

The Type I / II / III problem needs both ingredients: two or more predictors and imbalance. Take either one away — one predictor, or equal n — and it disappears.

And ecological data is almost never balanced: animals die, plots flood, a site gets sampled twice, a treatment fails, a field season ends early. The penguins above are only barely unbalanced — counts of 73/73, 34/34, 58/61 — and the sums of squares still move by 1.8 million.

Plan for equal n. Report Type II or III when you do not get it.

Part 5 · The F-Ratio

What F Actually Is

\[F = \frac{MS_{among}}{MS_{within}}\]

• \(MS_{within}\) estimates the background variance σ² — pooled from inside the groups, where the treatment cannot reach
• \(MS_{among}\) estimates that same σ², plus anything real that separates the group means
• If H₀ is true, both estimate σ² and F ≈ 1

Here F = 20.31 — the group means are 20 times further apart than within-group noise alone can explain.

Tip

🖐 Why F is never negative

It is a ratio of two variances, and a variance is a sum of squares. The direction of a difference disappears the moment you square it — which is why a significant F cannot tell you which way anything went.

F and t Are the Same Test

Take just two of the four groups — intact males and females — and test them both ways:

crab_two <- crab_df %>%
  filter(crab_type %in% c("intact male", "female")) %>%
  droplevels()

t.test(body_temp_c ~ crab_type, data = crab_two, var.equal = TRUE)$statistic

       t 
5.958321 

anova(lm(body_temp_c ~ crab_type, data = crab_two))[1, "F value"]

[1] 35.50158

t = 5.958, so t² = 35.5 — and F = 35.5.

Identical, with identical p-values.

Important

✅ ANOVA generalises the t-test

With two groups they are the same test. With three or more, the t-test has nowhere to go and ANOVA does.

This is the same t² = F identity you proved on a regression slope in Week 6.

Variance Explained

\[R^2 = \frac{SS_{among}}{SS_{total}} = \frac{2.64}{6.11} = 0.432\]

Crab type accounts for about 43.2% of the variation in body temperature. The rest is individual crabs being individual crabs.

summary(crab_model)$r.squared

[1] 0.4323687

Note

Significant and partial, at the same time

p = 7.00e-10 says the group differences are not noise. R² = 0.43 says they explain under half of what is going on.

Both are true, and a results section should carry both.

Part 6 · Assumptions and Diagnostics

The Three Assumptions

One-way ANOVA assumes:

1. Independence — observations do not depend on each other. This one comes from the design, and no plot can rescue it
2. Normality — of the residuals, not of the raw response
3. Equal variance (homoscedasticity) across groups

ANOVA is fairly robust to mild departures from 2 and 3 when the design is balanced, which is one more reason to plan for equal n.

Warning

⚠️ Normality of what?

Not body_temp_c. Not each group separately. The residuals of the fitted model — exactly as in regression.

A response that looks bimodal because the groups differ is perfectly fine; that is the effect you are testing.

The Four Diagnostic Plots

The same four panels you read for a regression:

1. Residuals vs fitted — with a factor, the points stack in vertical columns, one per group. You are looking for equal spread in each column
2. QQ plot — residuals against normal
3. Scale-location — same check, easier to see
4. Residuals vs leverage — influential points

Levene’s Test — Equal Variances

leveneTest(body_temp_c ~ crab_type, data = crab_df)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  3  0.1587 0.9238
      80               

H₀ is “the variances are equal”, so a large p-value is the good outcome — the same backwards logic as Shapiro-Wilk.

p = 0.924 — nothing to worry about.

Tip

🖐 Levene vs Bartlett

Bartlett’s test is more powerful but assumes normality itself. Levene’s works on absolute deviations from the group median, so it is far more robust. Use Levene.

Shapiro-Wilk — Normality of the Residuals

shapiro.test(residuals(crab_model))

    Shapiro-Wilk normality test

data:  residuals(crab_model)
W = 0.97566, p-value = 0.112

Note what goes inside: residuals(crab_model).

p = 0.112 — the residuals are consistent with normal.

Warning

⚠️ Do not let the test overrule the plot

At small n it cannot detect real problems; at large n it flags departures too small to matter. The QQ plot is the primary evidence.

One Call for All of It

performance::check_model() draws every assumption check in one labelled panel, with a plain-English caption on each.

check_model(crab_model)

Useful for a quick look. Still worth knowing what each panel is, which is why we did them one at a time first.

Part 7 · The Reference Level, Again

R Sorts Alphabetically — Still

Read the file without thinking about it and the four crab types become a factor in alphabetical order:

levels(crab_df$crab_type)

[1] "female"             "intact male"        "male major removed"
[4] "male minor removed"

So the reference is female — and every coefficient in the model becomes “compared with a female crab.”

Warning

⚠️ That is the wrong baseline for this experiment

The question is what happens when you remove a claw from a male. The comparison that answers it is against an intact male, not against a female.

The model fits identically either way. The output only answers your question if you tell it what your question was.

🗳️ DECIDE · Which Crab Is the Control?

You are about to set the reference level deliberately.

Which level should it be, and why?

1. female — keep R’s default; the model fit is identical anyway
2. intact male — the untouched version of the animals that were manipulated
3. male minor removed — it is the procedural control, so it accounts for the surgery
4. It makes no difference, because the F-test and R² do not change

Warning

✍️ Careful — two of these are true statements. Only one is the right decision.

Note

The question to ask

Which single comparison do I most want to read straight off the coefficient table?

fct_relevel() · Saying What You Meant

Important

✅ The answer is (b).

4. is true — F, R², and the ANOVA table are identical whichever level is the reference. That is not a reason to accept a default chosen by spelling.

5. is tempting and it is the right comparison to make eventually — but as the baseline it hides the very contrast the experiment was built around.

crab_df <- crab_df %>%
  mutate(crab_type = fct_relevel(crab_type, "intact male"))

crab_model <- lm(body_temp_c ~ crab_type, data = crab_df)

summary(crab_model)$coefficients %>% round(3)

                            Estimate Std. Error t value Pr(>|t|)
(Intercept)                    1.286      0.045  28.300    0.000
crab_typefemale                0.390      0.064   6.077    0.000
crab_typemale major removed    0.129      0.064   2.001    0.049
crab_typemale minor removed   -0.076      0.064  -1.186    0.239

Now every row reads as “degrees hotter than an intact male”:

• female +0.39 — clearly hotter
• male major removed +0.13
• male minor removed -0.08

Tip

🖐 The habit

Set the reference in the same line that creates the factor. fct_relevel(x, "control") moves one level to the front and leaves the rest alone.

Part 8 · Post-Hoc Testing

What Post-Hoc Tests Are For

The F-test said something differs. Post-hoc tests say what, while keeping the family-wide error rate at α.

• Planned (a priori) comparisons — a small number of contrasts you committed to before seeing the data. More powerful, because there are fewer of them
• Unplanned (post-hoc) comparisons — all pairs, after the fact. Safer, less powerful

Tukey’s HSD is the standard for all pairwise comparisons with equal n. We run it twice: once with TukeyHSD(), then with emmeans — the same answer, from the tool we use for every model after this one.

Warning

⚠️ Deciding after you look is not a plan

Choosing which comparisons to make because they looked big re-introduces exactly the multiple-comparisons problem the correction exists to solve.

Tukey’s HSD

TukeyHSD(aov(crab_model))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = crab_model)

$crab_type
                                             diff         lwr         upr
female-intact male                     0.39047619  0.22189182  0.55906056
male major removed-intact male         0.12857143 -0.04001294  0.29715579
male minor removed-intact male        -0.07619048 -0.24477484  0.09239389
male major removed-female             -0.26190476 -0.43048913 -0.09332040
male minor removed-female             -0.46666667 -0.63525103 -0.29808230
male minor removed-male major removed -0.20476190 -0.37334627 -0.03617754
                                          p adj
female-intact male                    0.0000002
male major removed-intact male        0.1962212
male minor removed-intact male        0.6375171
male major removed-female             0.0006106
male minor removed-female             0.0000000
male minor removed-male major removed 0.0108312

Each row is one pair: the difference in means, a confidence interval, and a p-value already adjusted for the whole family of six comparisons.

A pair differs when its interval excludes zero.

🗳️ DECIDE · Reading the Tukey Output

The radiator hypothesis predicts that a male whose major claw was removed should run hotter than an intact male.

Tukey gives male major removed − intact male a difference of +0.129 with an adjusted p = 0.196.

What should you conclude?

(a) The hypothesis is supported — the difference is in the predicted direction

(b) The hypothesis fails — removing the major claw did not significantly change temperature

(c) The hypothesis is supported, because females are significantly hotter than every male group

(d) Nothing can be concluded, because the overall ANOVA was significant

Warning

✍️ One letter, then argue it with your neighbour.

Note

The question to ask

Which comparison does the hypothesis actually make a prediction about?

What Tukey Actually Found

Important

✅ The answer is (b).

Removing the major claw moved temperature by 0.13 °C, adjusted p = 0.2 — not distinguishable from zero. The direction is right and the evidence is not there.

(c) is the trap: females are hotter than every male group, but females differ from males in many ways besides the claw. That comparison cannot test a claw-removal hypothesis.

What did differ (adjusted p):

Comparison	p
female − intact male	0
major removed − intact	0.196
minor removed − intact	0.638
major − minor removed	0.011

The big effect in this dataset is sex, not surgery.

Estimated Marginal Means · emmeans

TukeyHSD() is the end of the road. emmeans is the tool you will use for every model from here on. Step one is the estimated marginal means — the model’s estimate of each group’s mean:

crab_emm <- emmeans(crab_model, ~crab_type)
crab_emm

 crab_type          emmean     SE df lower.CL upper.CL
 intact male          1.29 0.0454 80     1.20     1.38
 female               1.68 0.0454 80     1.59     1.77
 male major removed   1.41 0.0454 80     1.32     1.50
 male minor removed   1.21 0.0454 80     1.12     1.30

Confidence level used: 0.95 

Note

🖐 Not the same as group_by() means

The SEs come from the pooled \(MS_{within}\), not each group’s own SD. These are the means you plot.

Tip

🖐 A lot of machinery for nothing?

With one factor each EMM equals its group mean exactly — only the SEs differ. That is what “averaging over the other factors” does when there are none.

Next week there will be, and that is where emmeans stops agreeing with group_by() %>% mean().

Pairwise Contrasts with emmeans

Step two: feed those means to pairs().

pairs(crab_emm, adjust = "tukey")

 contrast                                estimate     SE df t.ratio p.value
 intact male - female                     -0.3905 0.0643 80  -6.077 <0.0001
 intact male - male major removed         -0.1286 0.0643 80  -2.001  0.1962
 intact male - male minor removed          0.0762 0.0643 80   1.186  0.6375
 female - male major removed               0.2619 0.0643 80   4.076  0.0006
 female - male minor removed               0.4667 0.0643 80   7.263 <0.0001
 male major removed - male minor removed   0.2048 0.0643 80   3.187  0.0108

P value adjustment: tukey method for comparing a family of 4 estimates 

Same six comparisons, same adjusted p-values as TukeyHSD() — but watch the subtraction order, which runs the other way, so the signs flip. Wrap it in confint() if you want the intervals back.

Tip

🖐 Why emmeans and not TukeyHSD()

TukeyHSD() only does all-pairs on a balanced aov. emmeans handles unbalanced designs, covariates, several factors and planned contrasts — and it is what you will use from Week 9 onward.

adjust = also takes "bonferroni", "sidak", "BH" or "none".

Compact Letter Display

Six pairwise p-values is more than a reader can hold at once. cld() collapses them into letters:

multcomp::cld(crab_emm, Letters = letters,
              adjust = "tukey")

 crab_type          emmean     SE df lower.CL upper.CL .group
 male minor removed   1.21 0.0454 80     1.09     1.33  a    
 intact male          1.29 0.0454 80     1.17     1.40  ab   
 male major removed   1.41 0.0454 80     1.30     1.53   b   
 female               1.68 0.0454 80     1.56     1.79    c  

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 4 estimates 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Groups that share a letter were not shown to differ. These are the letters that go above the bars in your final figure.

Warning

⚠️ Do not library(multcomp)

It pulls in MASS, which masks dplyr::select() and will break your pipes with a baffling error three chunks later. Use the multcomp:: prefix and nothing gets attached.

Important

✅ A shared letter is not “the same”

intact male and male major removed share b: we failed to show a difference. That is not evidence they are equal.

cld() prints that warning itself, at the bottom of its own output. It is the most misread table in biology.

Part 9 · Reporting

Reporting an ANOVA

A results sentence needs the test, the statistic with both df, the p-value, and the direction — then the post-hoc detail.

Body temperature differed among crab types (one-way ANOVA, F~3,80~ = 20.31, p < 0.001, R² = 0.43). Females were warmer than intact males by 0.39 °C (Tukey HSD, adjusted p < 0.001), but removing the major claw did not detectably change male temperature (0.13 °C, adjusted p = 0.2).

Note

✅ What that sentence does

• Names the test
• Gives both degrees of freedom
• Reports effect size, not just p
• States the post-hoc method and that its p-values are adjusted
• Reports the non-significant comparison the hypothesis was about, rather than quietly dropping it

Part 10 · When the Assumptions Fail

A Second Dataset · Nematode Lifespan

The crab data behaved. Most data do not.

Here are 200 nematodes, 50 in each of four treatments, exposed to a parasite as larvae, as adults, as both, or to water as a control. The response is lifespan in days.

worm_df <- read_csv("data/nematode_lifespan.csv") %>%
  mutate(treatment = fct_relevel(treatment, "water"))

count(worm_df, treatment)

# A tibble: 4 × 2
  treatment       n
  <fct>       <int>
1 water          50
2 adult          50
3 larva          50
4 larva+adult    50

Note that the water control sorts last alphabetically — another reason to relevel on purpose.

These Assumptions Genuinely Fail

worm_model <- lm(lifespan_days ~ treatment, data = worm_df)

shapiro.test(residuals(worm_model))

    Shapiro-Wilk normality test

data:  residuals(worm_model)
W = 0.97626, p-value = 0.001808

leveneTest(lifespan_days ~ treatment, data = worm_df)

Levene's Test for Homogeneity of Variance (center = median)
       Df F value    Pr(>F)    
group   3  22.481 1.493e-12 ***
      196                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Shapiro-Wilk p = 0.00181 — residuals are not normal
• Levene p = 1.49e-12 — variances are nowhere near equal

Look back at the boxplots: group SDs run from 3.3 to 8.7 days.

Warning

⚠️ Nothing here is simulated

This is the published dataset, unmodified. Real data fails these tests routinely.

🗳️ DECIDE · What Do You Do Now?

An ANOVA on these data returns F = 15.25, p = 5.87e-09 — but the residuals are non-normal and the variances are unequal.

What is the right next step?

(a) Report the ANOVA — it is significant, and ANOVA is robust

(b) Drop the group with the biggest variance and re-run

(c) Try a transformation, and if that fails use a rank-based test such as Kruskal-Wallis

(d) Report nothing — the data cannot be analysed

Warning

✍️ Write your letter down, then vote.

Note

The question to ask

Is the problem with the data, or with the test I chose for it?

The Kruskal-Wallis Test

Important

✅ The answer is (c).

Transform first — a log or square root often fixes both problems at once, and keeps you in a framework that extends to more factors. If it does not, move to a rank-based test.

(b) is data deletion to make a test pass, which is not a repair.

kruskal.test(lifespan_days ~ treatment, data = worm_df)

    Kruskal-Wallis rank sum test

data:  lifespan_days by treatment
Kruskal-Wallis chi-squared = 29.463, df = 3, p-value = 1.79e-06

Kruskal-Wallis replaces every value by its rank and asks whether the rank totals differ more than chance — so it needs no normality and is not thrown by unequal spread.

χ² = 29.46, df = 3, p = 1.79e-06.

The conclusion survives. What changes is that you can now defend it.

Dunn’s Test — the Post-Hoc for Ranks

Kruskal-Wallis has the same limit as F: it says something differs, not what. Dunn’s test is its post-hoc.

dunnTest(lifespan_days ~ treatment, data = worm_df, method = "bh")$res %>%
       mutate(across(where(is.numeric), ~ signif(.x, 3)))

           Comparison      Z  P.unadj    P.adj
1       water - adult -3.120 1.78e-03 2.67e-03
2       water - larva  0.917 3.59e-01 4.31e-01
3       adult - larva  4.040 5.32e-05 1.60e-04
4 water - larva+adult -3.510 4.56e-04 9.11e-04
5 adult - larva+adult -0.381 7.03e-01 7.03e-01
6 larva - larva+adult -4.420 9.76e-06 5.85e-05

Note

Which adjustment?

"bonferroni" is the strictest and costs the most power. "bh" (Benjamini-Hochberg) controls the false discovery rate instead, and is the usual choice when you have more than a handful of comparisons.

pairwise.wilcox.test() does a similar job; Dunn’s is preferred because it uses the same ranking as the Kruskal-Wallis it follows.

Which Test Should You Use?

Situation	Parametric route	Non-parametric route
Compare 2 groups	t-test	Mann-Whitney U / Wilcoxon rank-sum
Compare 3+ groups	one-way ANOVA	Kruskal-Wallis
Post-hoc, all pairs	emmeans (Tukey-adjusted)	Dunn’s test, pairwise Wilcoxon
What it compares	means	mean ranks (≈ medians)
Needs	normal residuals, equal variance	neither
Costs you	nothing, when assumptions hold	power, and the original units

Try a transformation before you give up the parametric route — it keeps your effect sizes interpretable and extends to the multi-factor designs coming in Week 9.

Reporting a Non-Parametric Result

Nematode lifespan differed among parasite-exposure treatments (Kruskal-Wallis χ²~3~ = 29.46, p < 0.001). A Kruskal-Wallis test was used because the residuals of the equivalent ANOVA were non-normal (Shapiro-Wilk p = 0.00181) and variances were unequal (Levene p = 1.49e-12). Dunn’s tests with Benjamini-Hochberg adjustment showed that larva-exposed worms outlived worms exposed as adults or as both.

Important

✅ Say why you went non-parametric

A reader cannot tell whether you chose a rank test because your data required it or because it gave you the answer you wanted — unless you show the assumption checks that made the decision.

Report the numbers that drove the choice, not just the choice.

Wrap-up

Today you:

• Saw that ANOVA is regression with a categorical predictor — the same lm(), the same partition, the same F
• Worked out why six t-tests carry a 26% chance of a false positive, and what ANOVA does about it
• Built \(SS_{among}\) and \(SS_{within}\) by hand and matched them to anova()
• Confirmed t² = F on a two-group subset, for the second time this term
• Checked independence, normality and equal variance — the last two on the residuals
• Set the reference level with fct_relevel() so the coefficients answered the experiment’s question
• Used Tukey HSD and then emmeans to find which pairs differ, and learned that the claw-removal comparison was not one of them
• Moved to Kruskal-Wallis and Dunn’s test on a dataset whose assumptions genuinely failed

The Three Things to Carry Out of Here

Important

⚠️ A significant F is an entry ticket

It says at least one group differs. It never says which, how many, or in which direction. Everything past that comes from post-hoc comparisons you planned or corrected for.

Note

✅ Assumptions live in the residuals

Except independence, which lives in the design — and no plot will ever show it to you.

Important

✅ Set the reference level

fct_relevel(x, "control"), in the line that makes the factor. The fit never changes; the readability always does.

Note

Up next — Two-Way ANOVA

Two factors at once, and the thing that only exists when you have two: the interaction.

Readings

• Whitlock & Schluter, Ch. 15 — Comparing Means of More Than Two Groups (§15.4 assumptions and transformations; §15.5 planned and unplanned comparisons; §15.6 the Kruskal-Wallis test).
• Gotelli & Ellison, Ch. 10 — The Analysis of Variance.
• Darnell, M. Z., & Munguia, P. (2011). Thermoregulation as an alternate function of the sexually dimorphic fiddler crab claw. The American Naturalist, 178(3), 419–428.