ANOVA & Group Comparisons

1. The Concept: Why Not Just Run More t-Tests?

Suppose you measure PM2.5 at three sites (urban, suburban, rural) and want to know whether air quality differs by location. You could run three t-tests: urban vs. suburban, urban vs. rural, suburban vs. rural.

The problem: with three tests at α = 0.05, the probability of at least one false positive is \(1 - (1 - 0.05)^3 = 0.143\) — nearly 14%. With 10 groups the problem is catastrophic. ANOVA tests all groups simultaneously with a single α level, controlling the family-wise error rate.

📌 ANOVA’s null hypothesis \(H_0: \mu_1 = \mu_2 = \cdots = \mu_k\) (all group means are equal)

\(H_A\): at least one pair of means differs

ANOVA tells you whether any group differs, not which ones. Post-hoc tests answer the follow-up question.

2. Partitioning Variance

ANOVA decomposes the total variation in y into two components:

\[\underbrace{SS_\text{total}}_{\text{all variation}} = \underbrace{SS_\text{between}}_{\text{group effect}} + \underbrace{SS_\text{within}}_{\text{random error}}\]

Source	SS	df	MS	F
Between groups	\(\sum_j n_j(\bar y_j - \bar y)^2\)	\(k-1\)	\(SS_B / (k-1)\)	\(MS_B / MS_W\)
Within groups	\(\sum_j \sum_i (y_{ij} - \bar y_j)^2\)	\(N-k\)	\(SS_W / (N-k)\)	—
Total	\(\sum_j \sum_i (y_{ij} - \bar y)^2\)	\(N-1\)	—	—

The F-statistic is the ratio of between-group variance to within-group variance. Under \(H_0\), both estimate the same population variance, so F ≈ 1. A large F indicates the group means are more spread apart than chance alone would produce.

3. One-Way ANOVA in R

The dataset

We use plant weight data from a classic fertilizer experiment: three groups (control, treatment 1, treatment 2) with 10 plants each.

data(PlantGrowth)
str(PlantGrowth)

## 'data.frame':    30 obs. of  2 variables:
##  $ weight: num  4.17 5.58 5.18 6.11 4.5 4.61 5.17 4.53 5.33 5.14 ...
##  $ group : Factor w/ 3 levels "ctrl","trt1",..: 1 1 1 1 1 1 1 1 1 1 ...

table(PlantGrowth$group)

## 
## ctrl trt1 trt2 
##   10   10   10

Visualize first

ggplot(PlantGrowth, aes(x = group, y = weight, color = group)) +
  geom_boxplot(outlier.shape = NA, width = 0.4, linewidth = 0.8) +
  geom_jitter(width = 0.12, size = 2, alpha = 0.7) +
  scale_color_manual(values = c("#4a7c5e", "#c8832a", "#8b3a2a")) +
  labs(x = "Treatment group", y = "Dry weight (g)",
       title = "Plant growth by treatment") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "none")

Dry weight of plants by treatment group. Points are individual plants; box = IQR; line = median. Showing raw data alongside summaries is always best practice.

Fit the ANOVA

The recommended workflow in R is: fit with lm(), then pass to anova(). This keeps the output consistent with the regression framework you learned in Module 07.

mod_anova <- lm(weight ~ group, data = PlantGrowth)
anova(mod_anova)

## Analysis of Variance Table
## 
## Response: weight
##           Df  Sum Sq Mean Sq F value  Pr(>F)  
## group      2  3.7663  1.8832  4.8461 0.01591 *
## Residuals 27 10.4921  0.3886                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Equivalently, you can use aov():

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

##             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                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Both give identical results. The ANOVA table shows:

Df: degrees of freedom — 2 for between (3 groups − 1), 27 for within (30 obs − 3).
Sum Sq: the between-group and within-group sums of squares.
Mean Sq: SS / df for each source.
F value: the F-statistic = MS_between / MS_within.
Pr(>F): p-value. If < α, reject \(H_0\).

📌 Effect size: η² \(\eta^2 = SS_\text{between} / SS_\text{total}\) is the proportion of total variance explained by the group factor — analogous to R² in regression. Values of 0.01, 0.06, and 0.14 are conventionally described as small, medium, and large.

ss    <- anova(mod_anova)$"Sum Sq"
eta2  <- ss[1] / sum(ss)

4. Post-hoc Tests

A significant ANOVA F-test tells you that some groups differ. Post-hoc tests identify which pairs differ while controlling for multiple comparisons.

Tukey’s Honest Significant Difference (HSD)

Tukey’s HSD is the most common choice when you want to compare all pairwise combinations. It controls the family-wise error rate exactly at α.

tukey_res <- TukeyHSD(mod_aov, conf.level = 0.95)
tukey_res

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = weight ~ group, data = PlantGrowth)
## 
## $group
##             diff        lwr       upr     p adj
## trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
## trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
## trt2-trt1  0.865  0.1737839 1.5562161 0.0120064

par(las = 1, mar = c(4.5, 8, 3, 1))
plot(tukey_res, col = "#4a7c5e")

Tukey 95% family-wise confidence intervals for all pairwise group differences. Intervals that do not cross zero indicate a statistically significant difference.

Interpreting the output

Each row shows one pairwise comparison. The columns are:

Column	Meaning
`diff`	Estimated difference in group means
`lwr`, `upr`	95% confidence interval for the difference
`p adj`	Tukey-adjusted p-value

An interval that does not include zero, or a p-value < α, indicates the two groups are statistically different.

Bonferroni correction (manual)

For a smaller set of planned comparisons, the Bonferroni correction is straightforward: divide α by the number of comparisons.

pairwise.t.test(PlantGrowth$weight, PlantGrowth$group,
                p.adjust.method = "bonferroni")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  PlantGrowth$weight and PlantGrowth$group 
## 
##      ctrl  trt1 
## trt1 0.583 -    
## trt2 0.263 0.013
## 
## P value adjustment method: bonferroni

⚠️ Planned vs. unplanned comparisons If you specified your hypotheses before seeing the data, a smaller set of planned contrasts is appropriate and more powerful than Tukey HSD. If you are doing exploratory pairwise comparisons after a significant F-test, use Tukey HSD or Bonferroni to control inflation of the false positive rate.

5. Checking ANOVA Assumptions

One-way ANOVA requires: (1) independence of observations, (2) normally distributed residuals within each group, and (3) equal variances across groups (homoscedasticity).

Residual normality

par(las = 1, mar = c(4.5, 4.5, 3, 1))
qqnorm(residuals(mod_anova), main = "Q-Q Plot of ANOVA Residuals",
       col = "#4a7c5e", pch = 20)
qqline(residuals(mod_anova), col = "#c8832a", lwd = 2)

Q-Q plot of ANOVA residuals. Points close to the diagonal support the normality assumption.

shapiro.test(residuals(mod_anova))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_anova)
## W = 0.96607, p-value = 0.4379

Equal variances: Levene’s test

# Base R: Bartlett's test (sensitive to non-normality)
bartlett.test(weight ~ group, data = PlantGrowth)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  weight by group
## Bartlett's K-squared = 2.8786, df = 2, p-value = 0.2371

💡 ANOVA is robust to mild violations ANOVA is reasonably robust to moderate departures from normality, especially with equal group sizes. It is more sensitive to unequal variances when group sizes differ. If Levene’s or Bartlett’s test is significant, consider Welch’s one-way ANOVA (oneway.test()) which does not assume equal variances.

Welch’s ANOVA (robust to unequal variances)

oneway.test(weight ~ group, data = PlantGrowth, var.equal = FALSE)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  weight and group
## F = 5.181, num df = 2.000, denom df = 17.128, p-value = 0.01739

6. Nonparametric Alternative: Kruskal–Wallis Test

When the normality assumption is seriously violated and transformations do not help, the Kruskal–Wallis test is the nonparametric counterpart of one-way ANOVA. It ranks all observations and tests whether the rank distributions differ across groups.

kruskal.test(weight ~ group, data = PlantGrowth)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  weight by group
## Kruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842

Post-hoc: pairwise Wilcoxon tests with Bonferroni or BH correction.

pairwise.wilcox.test(PlantGrowth$weight, PlantGrowth$group,
                     p.adjust.method = "BH")

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  PlantGrowth$weight and PlantGrowth$group 
## 
##      ctrl  trt1 
## trt1 0.199 -    
## trt2 0.095 0.027
## 
## P value adjustment method: BH

7. Two-Way ANOVA

Two-way ANOVA tests the effects of two categorical predictors simultaneously and their interaction — whether the effect of one factor depends on the level of the other.

Model structure

# Use ToothGrowth: tooth length by supplement type and dose
data(ToothGrowth)
ToothGrowth$dose <- factor(ToothGrowth$dose)
str(ToothGrowth)

## 'data.frame':    60 obs. of  3 variables:
##  $ len : num  4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
##  $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
##  $ dose: Factor w/ 3 levels "0.5","1","2": 1 1 1 1 1 1 1 1 1 1 ...

mod_2way <- lm(len ~ supp * dose, data = ToothGrowth)
anova(mod_2way)

## Analysis of Variance Table
## 
## Response: len
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## supp       1  205.35  205.35  15.572 0.0002312 ***
## dose       2 2426.43 1213.22  92.000 < 2.2e-16 ***
## supp:dose  2  108.32   54.16   4.107 0.0218603 *  
## Residuals 54  712.11   13.19                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The ANOVA table now has three rows: main effect of supp, main effect of dose, and the supp:dose interaction. A significant interaction means the effect of supplement type differs across dose levels — the two factors are not independent.

Interaction plot

ToothGrowth |>
  group_by(supp, dose) |>
  summarise(mean_len = mean(len), se_len = sd(len) / sqrt(n()), .groups = "drop") |>
  ggplot(aes(x = dose, y = mean_len, color = supp, group = supp)) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = mean_len - se_len, ymax = mean_len + se_len),
                width = 0.12) +
  scale_color_manual(values = c("#4a7c5e", "#c8832a")) +
  labs(x = "Dose (mg/day)", y = "Tooth length (mm)",
       title = "Tooth length by supplement and dose",
       color = "Supplement") +
  theme_minimal(base_size = 13)

Interaction plot: tooth length by dose and supplement. Non-parallel lines indicate an interaction between the two factors — the supplement effect is larger at lower doses.

⚠️ If the interaction is significant, do not interpret main effects alone A significant interaction means the effect of each factor depends on the level of the other. Report and interpret the interaction first; the main effects become secondary descriptions. Ignoring a significant interaction and reporting only main effects is a common and serious misinterpretation.

8. Visualizing Group Differences

A well-designed group comparison plot always shows: (1) the raw data, (2) a measure of center, and (3) a measure of uncertainty.

PlantGrowth |>
  group_by(group) |>
  summarise(
    mean_w = mean(weight),
    ci_lo  = mean_w - qt(0.975, n() - 1) * sd(weight) / sqrt(n()),
    ci_hi  = mean_w + qt(0.975, n() - 1) * sd(weight) / sqrt(n()),
    .groups = "drop"
  ) -> pg_ci

ggplot(PlantGrowth, aes(x = group, y = weight)) +
  geom_jitter(width = 0.1, alpha = 0.45, size = 2,
              color = "#4a7c5e") +
  geom_point(data = pg_ci, aes(y = mean_w),
             size = 4, color = "#1e3a2f") +
  geom_errorbar(data = pg_ci, aes(y = mean_w, ymin = ci_lo, ymax = ci_hi),
                width = 0.18, linewidth = 1, color = "#1e3a2f") +
  labs(x = "Treatment group", y = "Dry weight (g)",
       title = "Plant growth: mean ± 95% CI with raw data") +
  theme_minimal(base_size = 13)

Recommended visualization for group comparisons: raw jittered points (low alpha), group mean (filled point), and 95% confidence interval (error bar).

🔑 Key takeaways

ANOVA tests whether any of k group means differ using a single F-test, avoiding the multiple-comparison inflation of repeated t-tests.
A significant F only tells you some groups differ. Use Tukey HSD to identify which pairs are different.
Check the LINE assumptions as you would for regression — ANOVA is the same model, just with a categorical predictor.
If variances are unequal, use oneway.test(var.equal = FALSE). If normality is seriously violated, use the Kruskal–Wallis test.
A significant interaction in two-way ANOVA must be reported and interpreted before main effects.

What’s next?

Simple Linear Regression — continuous predictor version of the same F-test framework
Logistic Regression — extend to binary outcomes
Mixed-Effects Models — handle repeated measures and nested designs

✏️ Practice Exercises

Exercise 1 — One-way ANOVA

The built-in PlantGrowth dataset has plant dry weights under three conditions: ctrl, trt1, trt2. (a) Plot the data (boxplot + jitter). (b) Run aov() and report the F-statistic and p-value. (c) Run Tukey HSD post-hoc tests. Which pairs differ significantly?

Show solution

data(PlantGrowth)
boxplot(weight ~ group, data = PlantGrowth, col = c("#8aad96","#4a7c5e","#1e3a2f"))

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

ANOVA: F ≈ 4.85, p ≈ 0.016. Tukey: trt1 vs. trt2 is the only significant pairwise difference (trt2 grows better than trt1; neither differs from ctrl at α = 0.05).

Exercise 2 — Assumption check

Test the homogeneity of variance assumption for the PlantGrowth ANOVA using bartlett.test(). Does the assumption hold? What would you do if it failed?

Show solution

bartlett.test(weight ~ group, data = PlantGrowth)

p ≈ 0.24 — no evidence that variances differ. If the assumption failed, you could use Welch’s ANOVA (oneway.test(…, var.equal=FALSE)) or a non-parametric Kruskal–Wallis test.

Simple Linear Regression ↑ Module Index Tests for Association

ANOVA & Group Comparisons

Marc Los Huertos

EA030 · Compiled 2026-04-14

ANOVA & Group Comparisons

1. The Concept: Why Not Just Run More t-Tests?

2. Partitioning Variance

3. One-Way ANOVA in R

The dataset

Visualize first

Fit the ANOVA

4. Post-hoc Tests

Tukey’s Honest Significant Difference (HSD)

Interpreting the output

Bonferroni correction (manual)

5. Checking ANOVA Assumptions

Residual normality

Equal variances: Levene’s test

Welch’s ANOVA (robust to unequal variances)

6. Nonparametric Alternative: Kruskal–Wallis Test

7. Two-Way ANOVA

Model structure

Interaction plot

8. Visualizing Group Differences

🔑 Key takeaways

What’s next?