Week 7: Comparing Groups

ANOVA and the multiple comparisons problem

Introduction

Last week you compared two groups with a t-test. What if you have four?

You could run every pairwise t-test — six of them for four groups. The content session showed you what that costs: at α = 0.05 each, the chance of at least one false positive when nothing is going on is 1 − 0.95⁶, about 26%. Run enough comparisons and you will always find something.

ANOVA asks a single question instead: does any group differ? If the answer is yes, a post-hoc test such as TukeyHSD says which pairs, with the correction for multiple looks built in.

Today’s data: wind farm capacity factors at 100 sites across four UK regions. Does geography matter for wind output?

Goals:

Visualize four groups at once, and order them usefully
Run an ANOVA and read its output
Find which pairs differ, with the correction applied
Understand what the correction does and does not buy you

How these exercises work

The code boxes start empty, except for comments to guide you through the code that’s expected.

Work one comment at a time. If you get stuck:

Hint 1 restates what you are trying to do.
Hint 2 sketches the shape the code might take.
Hint 3 gives you the code with the key pieces blanked out.
Solution gives you the lot.

Take a real attempt before opening a hint: you will learn by doing, not by reading. The demonstrators are here to help.

Three things worth knowing:

The boxes share one R session, exactly like a script. Anything you create in Exercise 1 is still there in Exercise 2.
Variable names are up to you. The last thing your block prints is what gets checked, so end each block with the answer.
Exercises have extensions underneath them. If you finish early, give them a stab. These are designed to equip you to tackle real coding problems you might face, using the documentation. We can help – but check the manuals first!

Stuck?

Ask on the class board, in Code Q&A. Someone else is likely stuck on the same thing, so asking in public helps them too. Demonstrators check it most weekdays, but they give classmates a chance to answer first.

The board is private to this class, so sign in to GitHub first. If the link shows “404 – page not found”, you’re not signed in, or you haven’t accepted the Classroom 50 invitation yet.

Exercise 1: Four groups at a glance

wind has 100 rows: a region and a capacity_pct — the percentage of its theoretical maximum output that a site actually delivered.

Draw violins by region, ordered by median rather than alphabetically, and finish by printing the four group means.

Hint 2

The shape of it:

ggplot(df, aes(x = reorder(group, value, FUN = median), y = value)) +
  geom_violin() +
  labs(x = "...", y = "...")

tapply(df$value, df$group, mean)

Reordering inside aes() leaves the data alone and changes only the plot. That is usually what you want — the alternative, converting the column to a factor with explicit levels, affects every later use of it.

Hint 3

ggplot(wind, aes(x = reorder(______, capacity_pct, FUN = ______),
                 y = ______)) +
  geom_violin() +
  labs(x = "Region",
       y = "Capacity factor (% of maximum output)",
       title = "Wind farm capacity factor by UK region")

tapply(wind$capacity_pct, wind$______, ______)

Solution

ggplot(wind, aes(x = reorder(region, capacity_pct, FUN = median),
                 y = capacity_pct)) +
  geom_violin() +
  labs(x = "Region",
       y = "Capacity factor (% of maximum output)",
       title = "Wind farm capacity factor by UK region")

tapply(wind$capacity_pct, wind$region, mean)

Two clusters, and ordering the axis is what makes them visible:

Higher: Northern England 29.9%, Scotland 29.6%
Lower: Wales 24.7%, Southern England 23.2%

That is geographically plausible — northern and western sites are windier. But the gaps are very different in size. Six points separate the clusters; three-tenths of a point separates Northern England from Scotland. The question is not “are the four numbers identical”, because four sample means never are. The question is which of these gaps is bigger than the sampling noise.

Note also that the violins are much the same height. ANOVA assumes roughly equal variance across groups, and you have just checked it without being asked to.

Extension: show the sites

Optional. ?geom_jitter, ?geom_violin, ?stat_summary, ?geom_errorbar.

Overlay the 25 individual sites on each violin. Does any region look like it has sub-groups within it? A bimodal region would make its mean nearly meaningless, and a violin’s smoothing can blur sub-groups away (or invent them); the points can’t.
Now replace the violins themselves with means and an approximate 95% confidence interval, but keep the jittered points from step 1 showing behind them. That combination — raw points plus mean ± CI — is what ANOVA actually reasons about, and it is noticeably less reassuring than the violins: the points overlap heavily between regions even where the CIs do not. stat_summary(fun.data = mean_se, fun.args = list(mult = 1.96)) draws it in one line (mean_cl_normal() would be the exact t-based version, but it needs the Hmisc package, which isn’t available here — 1.96 standard errors is close enough at 25 sites per group).
Add a horizontal line at the overall mean of 26.8%. ANOVA compares the spread of group means around that line with the spread of sites around their own group mean, so this one line makes the F statistic visible.
Order the regions north to south instead of by output, and decide which ordering you would publish. There is a real argument for each.

Further reading: on why we don’t use box plots — they are harder to read than they look — and what to show instead — Nick Desbarats, “I’ve Stopped Using Box Plots. Should You?”, Nightingale, and Claus Wilke, “Visualizing uncertainty”, Fundamentals of Data Visualization, ch. 9.

Exercise 2: One test instead of six

ANOVA compares the variation between group means with the variation within groups. If the group means are further apart than the noise inside each group would explain, F is large and p is small.

Run it, and print the summary.

Hint 2

The shape of it:

fit <- aov(outcome ~ group, data = df)
summary(fit)

Save the fit to a name. TukeyHSD in the next exercise needs the aov object itself, not its summary, and re-running the fit there would work but would also invite the two to drift apart.

Solution

wind_aov <- aov(capacity_pct ~ region, data = wind)
summary(wind_aov)

F = 14.21 on 3 and 96 degrees of freedom, p = 9.8 × 10⁻⁸.

The degrees of freedom are 3 (four groups, minus one) and 96 (100 sites, minus four group means). If those numbers are not what you expected, something is wrong with the data before anything is wrong with the test — a common cause is R having read a grouping column as something other than four distinct values.

F = 14.2 means the group means are spread about fourteen times more widely than the within-region scatter would produce by chance. The p-value is tiny.

And it tells you almost nothing you can use. “At least one of four regions differs from at least one other” is not a planning recommendation. You already suspected it from the plot. What ANOVA buys is permission to go looking for the specific differences without the 26% false-positive problem — which is the next exercise.

Extension: what ANOVA assumes

Optional. ?plot.lm, ?kruskal.test, ?oneway.test, ?bartlett.test.

plot(wind_aov) gives four diagnostic plots. The first two are the ones to read: residuals against fitted values (look for a fan shape) and a QQ plot of residuals (look for curvature). Does this ANOVA deserve to be trusted?
Run kruskal.test() on the same data — the rank-based version that assumes no normality. It gives p = 1.5 × 10⁻⁶. Both are tiny, so the conclusion survives, which is the most useful possible outcome of a robustness check and the least reportable.
oneway.test() is the ANOVA equivalent of Welch’s t-test, dropping the equal-variance assumption. Compare its answer. Then work out from the group standard deviations why it barely moves here.
ANOVA with exactly two groups is a t-test. Confirm it: run both on Scotland and Wales alone, and check that F equals t².

Exercise 3: Which pairs actually differ?

ANOVA said something differs. TukeyHSD says what, comparing all six pairs while holding the overall false-positive rate at 5% rather than letting it climb to 26%.

Run it, then finish by printing how many of the six pairs are significant at the 0.05 level.

Hint 2

The shape of it:

tukey <- TukeyHSD(fit)
tukey

sum(tukey$grouping_variable[, "p adj"] < 0.05)

The column is called p adj, with a space, so it has to be quoted inside the brackets. The comma before it means “all rows, this column”.

Solution

tukey_result <- TukeyHSD(wind_aov)
tukey_result

sum(tukey_result$region[, "p adj"] < 0.05)

Four of the six:

Comparison	diff	p adj
Scotland – Northern England	−0.33	0.994	no difference
Wales – Southern England	+1.43	0.673	no difference
Wales – Northern England	−5.21	0.0005	differs
Wales – Scotland	−4.88	0.0012	differs
Southern England – Scotland	−6.32	<0.001	differs
Southern England – Northern England	−6.64	<0.001	differs

Every north–south pair differs; neither within-cluster pair does. The finding is not “region matters” but something sharper and more useful: latitude matters, and the distinction between Scotland and Northern England is not supported by these data. A developer choosing between those two on capacity factor alone would be choosing on noise.

One honest footnote. Run the six comparisons uncorrected and you get the same four. The correction changed nothing here, because these p-values are not close to the threshold. That is the usual case, and it is exactly why the correction is easy to skip: it looks like it never matters, right up until the study where your best pair sits at 0.03 and the corrected value is 0.18.

Extension: the correction, and what it costs

Optional. ?p.adjust, ?pairwise.t.test, ?plot.TukeyHSD.

Plot the Tukey result. Six horizontal intervals; the two that cross the zero line are the two non-significant pairs. Then work out why this figure is more use in a report than the table is.
Run pairwise.t.test() with p.adjust.method = "none", then "bonferroni", then "holm". Four comparisons stay significant throughout. Construct — or simulate — a dataset where they do not agree, which is the only way to see what the correction is for.
Corrections are not free: holding the false-positive rate down raises the false-negative rate. With 28 pairs instead of 6, how large would a real difference have to be to survive Bonferroni? ?power.t.test with a smaller sig.level will show you.
The deeper problem the correction does not solve: it assumes you declared your six comparisons in advance. Describe, in one sentence, how someone could run this exact analysis honestly and still end up reporting a false positive.

Quick check

With four groups there are six pairwise comparisons. How many are there with eight groups? Print the number.

Solution

choose(8, 2)

Doubling the groups nearly quintuples the comparisons, so the false-positive problem grows much faster than the study does:

Groups	Pairs	Chance of at least one false positive at α = 0.05
2	1	5%
4	6	26%
8	28	76%
12	66	97%

At 12 groups you are essentially guaranteed a “finding”. This is the arithmetic behind the HolmesCo mineral survey, and behind a good deal of published science.

Save your work

Copy the code you are most proud of into your week7.R file. Commit and push via GitHub Desktop. Write a commit message that describes what you found — not just “week 7”.