ANOVA in R: aov(), the ANOVA Table and Tukey HSD

Why ANOVA and Not a Pile of T-Tests

A t-test compares two means. With three or more groups, the tempting move is to t-test every pair - and that's exactly the trap. Each test run at the 0.05 level carries a 5% false-positive risk, and the risks accumulate: with 4 groups that's 6 pairwise tests and roughly a 26% chance of at least one spurious "significant" result; with 5 groups (10 tests), roughly 40%. You'd manufacture discoveries out of pure noise.

ANOVA (ANalysis Of VAriance) fixes this by asking one question with one p-value: are all the group means the same, or does at least one differ? Despite the name, it compares means - it just does so by analyzing variance: if the group averages are more spread out than the noise inside groups can explain, something real is going on.

One-Way ANOVA with aov()

PlantGrowth is built for this: plant weights under a control and two treatment conditions. Fit with aov(), then - important - print the table with summary():

fit <- aov(weight ~ group, data = PlantGrowth)
summary(fit)

Read the formula as "does weight depend on group?" The grouping variable must be a factor - PlantGrowth$group already is one, but if your groups are coded as numbers (doses, batch IDs), wrap them: aov(y ~ factor(dose), ...). Otherwise aov() silently fits a regression line through the group codes instead of comparing group means - wrong model, no error message.

Reading the ANOVA Table

The table is two rows - the factor and the residuals - and five columns. Cell by cell:

            Df Sum Sq Mean Sq F value Pr(>F)
group        2  3.766  1.8832   4.846 0.0159
Residuals   27 10.492  0.3886

• Df - degrees of freedom. The factor row gets groups − 1 (3 groups → 2); the residual row gets observations − groups (30 − 3 = 27). A quick sanity check that R saw the design you intended.
• Sum Sq - variation, split into two piles. The group row is the between-group pile: how far the group means sit from the overall mean. Residuals is the within-group pile: how much plants vary around their own group's mean. Together they add up to the total variation in the data.
• Mean Sq - each Sum Sq divided by its Df, turning piles of variation into comparable per-degree-of-freedom rates. The residual Mean Sq (0.389) is the noise level.
• F value - the ratio: Mean Sq of group over Mean Sq of Residuals (1.8832 / 0.3886 ≈ 4.85). If all group means were truly equal, this ratio hovers around 1. The bigger F gets, the more the between-group differences outrun what noise can explain.
• Pr(>F) - the p-value: the probability of an F this large if all three true means were identical. Here 0.016 - small enough at the conventional 0.05 level to reject "all equal." As always, it's not the probability the null is true, and it says nothing about which groups differ or by how much.

That last point is the crucial limitation: a significant F says "some difference, somewhere." Nothing more.

Which Groups Differ? TukeyHSD()

The follow-up question needs a post-hoc test. Tukey's Honest Significant Difference tests every pair while keeping the family-wise error rate at 5% across all comparisons at once:

fit <- aov(weight ~ group, data = PlantGrowth)
TukeyHSD(fit)

One row per pair, four numbers per row:

• diff - the estimated difference in means (second-named group minus first).
• lwr, upr - the 95% family-wise confidence interval for that difference.
• p adj - the p-value, already adjusted for making three comparisons.

For PlantGrowth: trt2-trt1 shows a difference of about 0.87 with p adj ≈ 0.012 and an interval clear of zero - treatment 2 outgrows treatment 1. Both treatment-vs-control rows (trt1-ctrl, trt2-ctrl) have intervals straddling zero and p adj well above 0.05 - neither treatment is distinguishable from control on this sample. The rule of thumb mirrors confidence intervals everywhere: interval excludes zero ⇔ p adj below 0.05.

Notice how the overall F said "a difference exists" while Tukey locates it in exactly one pair - the two-step structure is the whole design: one honest overall test, then properly-corrected pairwise detective work.

Two-Way ANOVA: Two Factors and Their Interaction

With two grouping variables, one call tests both - plus whether they interact. ToothGrowth crosses supplement type (supp) with dose (0.5, 1, 2 mg - numeric, so it needs factor()):

fit2 <- aov(len ~ supp * factor(dose), data = ToothGrowth)
summary(fit2)

supp * factor(dose) expands to three effects, one table row each:

• supp - averaged over doses, does supplement type matter? (Yes: p ≈ 0.0002.)
• factor(dose) - averaged over supplements, does dose matter? (Emphatically: p is tiny.)
• supp:factor(dose) - the interaction: does the effect of supplement change depending on dose? Here p ≈ 0.022 - it does. Orange juice beats ascorbic acid at low doses, but the gap closes at the 2 mg dose.

A significant interaction is a warning label on the main effects: "supplement type matters" is only true on average, and the average hides a dose-dependent story. When the interaction is significant, describe the combinations (a TukeyHSD(fit2, "supp:factor(dose)") or a table of group means via group-wise summaries) rather than reporting the main effects alone. Use + instead of * only when you deliberately want a no-interaction model.

Assumptions - and the Rank-Based Fallback

ANOVA's math leans on three assumptions, in decreasing order of negotiability:

• Independence - observations don't influence each other. Nothing fixes a violation after the fact; it's a property of the study design.
• Equal variances across groups - the residual Mean Sq is one pooled noise estimate, so groups should be about equally noisy. One-line check: bartlett.test(weight ~ group, data = PlantGrowth) (a large p-value means no evidence of unequal variances).
• Roughly normal residuals - matters least with balanced designs and decent group sizes; eyeball a histogram or boxplot of residuals via residuals(fit).

When the data is heavily skewed, ordinal, or riddled with outliers, the rank-based alternative to one-way ANOVA is kruskal.test(weight ~ group, data = PlantGrowth) - the many-group sibling of wilcox.test(). It's one line, and it trades some power for robustness.

What You Take Away

• ANOVA asks one question about 3+ group means with one p-value - the fix for t-test multiplicity.
• fit <- aov(y ~ group, data = df); summary(fit) - and the grouping variable must be a factor.
• The F value is between-group signal over within-group noise; Pr(>F) small means "some difference somewhere," and nothing more.
• TukeyHSD(fit) finds which pairs differ, with the multiple-comparison correction built in.
• a * b fits two factors plus their interaction; a significant interaction means the main effects don't tell the story alone.
• Check equal variances (bartlett.test) and residual normality; kruskal.test() is the rank-based fallback.

Next up: from comparing group means to modeling a relationship - linear regression with lm().