Compare Several Means with ANOVA

One-way ANOVA with residual checks and Tukey’s HSD, then chi-square for counts in categories

Coding 8 compared two means. Lecture A010 added two tests for other kinds of data, and this activity runs both: one-way ANOVA for a continuous response across three or more groups, and chi-square for counts in categories. Same recipe as always: H₀ as math, α first, look, check, test, conclude.

Part A uses Palmer Penguins. Part B uses the exact counts from the slides.

Which test? Match the data type (click to expand)

Your data	Test	H₀ as math	R
Continuous response, 2 groups	t-test (Welch, or paired)	μ₁ = μ₂	`t.test()`
Continuous response, 3+ groups	One-way ANOVA, then Tukey HSD	μ₁ = μ₂ = μ₃ = …	`aov()`, `TukeyHSD()`
Counts in categories, one variable vs. a predicted ratio	χ² goodness of fit	counts follow the ratio	`chisq.test(x, p = ...)`
Counts in categories, two variables	χ² test of independence	the variables are unrelated	`chisq.test(table)`

Pick the test from the data type and the number of groups. The research question tells you which row you are in.

ANOVA and chi-square functions (click to expand)

Function	What it does	Example
`aov(y ~ group, data = )`	Fits a one-way ANOVA; H₀ is that every group mean is equal	`aov(bill_depth_mm ~ species, data = penguin_data)`
`summary(model)`	Prints the F statistic and its p-value	`summary(anova_model)`
`residuals(model)`	Each value minus its own group mean, one number per row	`residuals(anova_model)`
`shapiro.test(residuals(model))`	Normality check on the residuals, one check for the whole model
`fligner.test(y ~ group, data = )`	Equal-variance check across groups; H₀ is equal variances	`fligner.test(bill_depth_mm ~ species, data = penguin_data)`
`TukeyHSD(model)`	Every pairwise difference with the family-wise error rate held at α	`TukeyHSD(anova_model)`
`choose(k, 2)`	How many pairs k groups make	`choose(6, 2)` is 15
`chisq.test(x, p = )`	Goodness of fit: observed counts vs. a predicted ratio	`chisq.test(c(428, 152), p = c(3/4, 1/4))`
`chisq.test(table)`	Independence: are two categorical variables related?	`chisq.test(infection_table)`
`matrix(..., nrow = 2, byrow = TRUE)`	Builds a contingency table from counts	see step 15

shapiro.test(), fligner.test(), and chisq.test() all return the same kind of object as t.test(): pull pieces out with $statistic, $parameter, and $p.value.

Part A: One-way ANOVA

Research question: Does mean bill depth differ among Adelie, Chinstrap, and Gentoo penguins?

H₀: μ_Adelie = μ_Chinstrap = μ_Gentoo. Hₐ: at least one mean differs. One-way ANOVA, F statistic, α = 0.05. Then, if the ANOVA is significant, Tukey’s HSD to see which species differ.

1. Load the packages

2. Build the analysis table

Keep every penguin with a non-missing bill_depth_mm. Save the result as penguin_data.

Check yourself

How many factors are in this analysis, and how many levels?
What is the response variable?
Why not run three t-tests (Adelie vs. Chinstrap, Adelie vs. Gentoo, Chinstrap vs. Gentoo)?

Answers: one factor (species) with three levels; bill_depth_mm; each test has a 5% false-positive rate, and the chances add up. Step 10 puts a number on it.

3. Boxplot by species

Make a boxplot of bill_depth_mm by species. Save it as box_species and display it.

Solution

box_species <- ggplot(penguin_data, aes(x = species, y = bill_depth_mm)) +
  geom_boxplot()
box_species

Gentoo bills look shallower than the other two. Adelie and Chinstrap look alike. The ANOVA will say whether chance could produce this.

4. Histogram of every group

Make a histogram of bill_depth_mm with one panel per species using facet_wrap(~ species). Save it as hist_species and display it.

Solution

hist_species <- ggplot(penguin_data, aes(x = bill_depth_mm)) +
  geom_histogram(binwidth = 0.5) +
  facet_wrap(~ species)
hist_species

Each species on its own looks roughly bell-shaped. Stacked together they would look bimodal, because Gentoo sits 3 mm lower. That is why ANOVA checks normality on the residuals, not on the pooled response.

5. Fit the ANOVA

Fit a one-way ANOVA of bill_depth_mm by species with aov(). Save the model as anova_model and print summary(anova_model).

Solution

anova_model <- aov(bill_depth_mm ~ species, data = penguin_data)
summary(anova_model)

F = 360 on 2 and 339 df, p < 2 × 10⁻¹⁶. The spread of the three group means is enormous compared with the spread within groups.

Check yourself

F = 360, p < 2 × 10⁻¹⁶. What have you learned?

Every species has a different mean bill depth
At least one species has a different mean bill depth
Gentoo bills are the shallowest

Answer: B. The F test says “at least one differs”. It does not say which, and it does not rank them. That is Tukey’s job in step 11.

6. Add the residuals

A residual is each value minus its own group mean. residuals(anova_model) returns one per penguin, in the same row order as penguin_data. Add them as a column named resid with mutate() and save the result back as penguin_data.

Solution

penguin_data <- penguin_data |>
  mutate(resid = residuals(anova_model))
penguin_data

A new column. An Adelie with a 19.3 mm bill has a residual of about +1.0, because the Adelie mean is 18.3.

7. Histogram of the residuals

Make a histogram of resid in penguin_data. Save it as hist_resid and display it.

Solution

hist_resid <- ggplot(penguin_data, aes(x = resid)) +
  geom_histogram(binwidth = 0.25)
hist_resid

One bell centered on zero, running from about −3 to +3 mm. No skew, no second hump.

8. Shapiro-Wilk on the residuals

Run shapiro.test() on residuals(anova_model). Save it as shapiro_resid and print it.

Solution

shapiro_resid <- shapiro.test(residuals(anova_model))
shapiro_resid

shapiro.test(penguin_data$resid) gives the same result. W = 0.992, p = 0.047, just under 0.05. Read the callout below before you react.

p = 0.047 is not an order

The histogram is a clean bell. W = 0.992 is about as close to 1 as the statistic gets. The p-value is small because n = 342: with that many values, Shapiro-Wilk can detect a departure from normality far too small to matter. ANOVA is robust to this kind of departure when the groups are this large and this similar in shape.

Decision: proceed with the ANOVA and report both facts. “Residuals were approximately normal (histogram symmetric; Shapiro-Wilk W = 0.99, p = 0.047 with n = 342).” A test result is one piece of evidence, not the verdict. Coding 7 made the same point in the other direction, when p = 0.029 on 49 values turned out to be a mild skew that changed nothing.

9. Fligner test for equal variances

Check whether the three species have similar spread with fligner.test(). Save it as fligner_species and print it.

10. Why not three t-tests?

Three species make three pairs. If each pair gets its own t-test at α = 0.05, what is the chance of at least one false positive when H₀ is true for every pair? Compute n_pairs with choose(3, 2) and p_false_positive as 1 - 0.95^n_pairs. Print both.

Solution

n_pairs <- choose(3, 2)
p_false_positive <- 1 - 0.95^n_pairs
n_pairs
p_false_positive

3 pairs, 14% chance of at least one false positive. For the six feeds in chickwts it is 15 pairs and 54%. Tukey’s HSD keeps the rate for the whole set of comparisons at 5%.

11. Tukey’s HSD

The ANOVA said at least one species differs. Run TukeyHSD() on anova_model to find out which. Save it as tukey_result and print it.

Solution

tukey_result <- TukeyHSD(anova_model)
tukey_result

Three rows, one per pair. diff is the difference in means, lwr and upr are a 95% family-wise confidence interval, and p adj is the p-value after correction. Chinstrap − Adelie: 0.07 mm, p adj = 0.89. Gentoo − Adelie: −3.36 mm, p adj < 0.0001. Gentoo − Chinstrap: −3.44 mm, p adj < 0.0001.

12. Make the call

Using tukey_result, assign two logicals: adelie_chinstrap_differ (do Adelie and Chinstrap differ at α = 0.05?) and gentoo_shallower (is the Gentoo mean below both others, and significantly so?).

Solution

adelie_chinstrap_differ <- FALSE
gentoo_shallower <- TRUE

The overall ANOVA was significant, but only Gentoo stands apart. Adelie and Chinstrap bill depths are statistically indistinguishable (p adj = 0.89). Gentoo bills are about 3.4 mm shallower than both.

Part B: Chi-square

When the response is a count in categories, there is no mean to compare. Chi-square compares the counts you observed with the counts H₀ predicts. Two versions: goodness of fit (one categorical variable; expected counts come from a hypothesis) and independence (two categorical variables; expected counts assume the variables are unrelated).

13. Goodness of fit: Mendel’s pea pods

All parents were heterozygous for pod colour and green is dominant, so H₀ predicts 3 green : 1 yellow in the F₂. Mendel observed 428 green and 152 yellow. Create observed_counts as a named vector and run chisq.test() with p = c(3/4, 1/4). Save it as chisq_gof and print it.

Solution

observed_counts <- c(green = 428, yellow = 152)
chisq_gof <- chisq.test(observed_counts, p = c(3/4, 1/4))
chisq_gof

χ² = 0.45, df = 1, p = 0.50. Expected counts were 435 and 145; the observed 428 and 152 are well within chance.

14. Make the call

Do these counts support Mendel’s 3 : 1 prediction at α = 0.05? Assign supports_3_to_1 to TRUE or FALSE.

Solution

supports_3_to_1 <- TRUE

p = 0.50, so we fail to reject H₀. The counts are consistent with 3 : 1. That is support, not proof: a ratio of 2.9 : 1 would have fit about as well.

15. Independence: lizard infection by habitat

Sixty lizards from an urban site and sixty from a forest site were scored infected or clean. The contingency table is built for you with matrix(). H₀: infection status is independent of habitat. Run chisq.test() on infection_table with correct = FALSE, as on the slide. Save it as chisq_indep and print it.

Hint

Pass the whole table. No p =, because the expected counts come from the row and column totals.

chisq.test(my_table, correct = FALSE)

correct = FALSE turns off Yates’ continuity correction, which R applies by default to 2 × 2 tables. The slide ran without it.

Solution

chisq_indep <- chisq.test(infection_table, correct = FALSE)
chisq_indep

χ² = 7.67, df = 1, p = 0.0056. Under independence each habitat would have 25.5 infected lizards; urban had 18 and forest had 33.

16. Make the call

Urban lizards: 30% infected. Forest lizards: 55%. Is infection associated with habitat at α = 0.05? Assign infection_linked_to_habitat to TRUE or FALSE.

Solution

infection_linked_to_habitat <- TRUE

p = 0.0056 < 0.05. Reject H₀: infection is associated with habitat, and forest lizards are infected more often. It does not show that living in the forest causes infection. Habitat was not randomly assigned; forest lizards may differ in age, diet, or density.

17. Recap: functions you practiced

What does each function do?

Function	Role
`filter(!is.na(...))`	Drop rows with a missing response
`ggplot()` + `geom_boxplot()`	Compare a continuous response across factor levels
`geom_histogram()` + `facet_wrap(~ species)`	One distribution panel per level
`aov(y ~ group, data = )`	Fit a one-way ANOVA
`summary(anova_model)`	F statistic and p-value: does any mean differ?
`mutate(resid = residuals(anova_model))`	Attach residuals to the table
`geom_histogram()` on `resid`	Visual normality check for the whole model
`shapiro.test(residuals(anova_model))`	Formal normality check on the residuals
`fligner.test(y ~ group, data = )`	Equal-variance check
`choose(k, 2)` and `1 - 0.95^n`	How many pairs, and the family-wise false-positive rate
`TukeyHSD(anova_model)`	Which pairs differ, with the rate held at 5%
`chisq.test(x, p = )`	Goodness of fit against a predicted ratio
`chisq.test(table, correct = FALSE)`	Independence of two categorical variables

If you can explain each row in your own words, you have covered Coding Assignment 4, Part C and the chi-square half of lecture A010.

Three things to say in your write-up

ANOVA says “at least one”; Tukey says which. Reporting “Gentoo had the shallowest bills (ANOVA p < 0.001)” is wrong. The ANOVA p-value never tested Gentoo against anything. Cite the Tukey contrast: Gentoo − Adelie = −3.36 mm, 95% family-wise CI −3.68 to −3.04, p adj < 0.0001.
Checks are evidence, not verdicts. The residual histogram was a clean bell; Shapiro-Wilk gave p = 0.047 on 342 values; Fligner gave p = 0.13. Report all three and your decision to proceed.
Chi-square shows association. Forest lizards were infected more often (55% vs. 30%, χ² = 7.67, p = 0.0056). Nothing in the test says habitat is the cause.

18. Transfer to Coding Assignment 4

On DataHub, Part C fits aov(biomass_g ~ nitrogen, data = nitrogen_data) with four nitrogen levels. The assignment asks for plot(anova_model, which = 1), a base-R residuals-versus-fitted plot, instead of the residual histogram you built in step 7. Try it on the penguin model now so the picture is familiar:

Three vertical stripes, one per species mean, each a cloud of residuals centered on zero. Even spread across the stripes is the equal-variance assumption in picture form. Use only which = 1; the assignment does not ask you to interpret the other panels.

Then fligner.test(), TukeyHSD(), and a write-up naming at least one significant and one non-significant nitrogen contrast. Same functions as steps 5 through 12.

Keep playing

Six feeds, fifteen pairs. Run the lecture’s ANOVA on chickwts and read the Tukey table from the slide.

F = 15.4, p = 6 × 10⁻¹⁰. Of the 15 pairs, 8 are significant. Sunflower and casein do not differ from each other (p adj = 0.9999) but beat horsebean, linseed, and soybean. That is the answer to “which feed would you give your chicks?”

When the variance check fails. Try body_mass_g instead of bill depth.

Fligner p = 0.01. Body mass varies more in Gentoo than in the other two species, so the equal-variance assumption fails. That is why this activity used bill depth. On the assignment, report the Fligner result either way; a failed check is a finding, not a failure.

Which test? For each design, name the test.

Check yourself

Seedlings grown under four light treatments, 15 pots each, height measured after 30 days.
For each of 200 nests, whether it was in a nest box or a natural cavity, and whether it fledged at least one chick.
A bag of 500 M&Ms sorted by colour, compared with the proportions the company publishes.
A one-way ANOVA on five diets gives p = 0.03, and a student reports “Diet 3 produced the heaviest mice (p = 0.03).” What is wrong?

Answers: (1) one-way ANOVA, then Tukey HSD if significant. (2) χ² test of independence: two categorical variables, one count per combination. (3) χ² goodness of fit: one categorical variable against a predicted ratio. (4) The ANOVA p-value does not test diet 3 against the others; a post-hoc test is needed before naming a winner.