ANOVA & Experimental-Design Basics

1Objectives¶

By the end of this lesson you will be able to:

Fit a one-way ANOVA model with aov(y ~ g, data=) and read its table with anova(model).
Interpret an $F$ -statistic and its p-value to decide whether group means differ.
Define the core experimental-design vocabulary: factor, levels, randomization, replication.
Connect a small designed experiment to the ANOVA table it produces.
Explain the multiple-comparisons problem — why a significant omnibus ANOVA doesn’t say which groups differ — and run and interpret a post-hoc pairwise comparison with TukeyHSD().

2From two groups to several¶

L11 compared exactly two group means with a t-test. ANOVA (“Analysis of Variance”) answers the same kind of question — do group means differ? — for three or more groups at once, with a single test instead of running t-tests on every pair. survey’s five declared majors are a ready-made example: does mean Exam 1 score differ across Business, Kinesiology, Nursing, Other, and STEM?

suppressMessages({library(mosaic); library(BSDA)})
set.seed(2200)
survey <- read.csv("data/survey_sim.csv")
survey$major_area <- factor(survey$major_area)

31. EDA first, as always¶

favstats(exam_score ~ major_area, data = survey)

   major_area  min     Q1 median    Q3  max     mean        sd  n missing
1    Business 42.2 56.150  60.30 66.40 76.7 60.88065  8.811145 31       0
2 Kinesiology 40.0 60.250  63.40 69.25 86.1 64.42500 10.353166 24       0
3     Nursing 49.3 59.200  62.20 68.55 78.5 63.37895  7.844430 19       0
4       Other 43.6 51.475  55.85 65.95 82.6 58.84091 10.820160 22       0
5        STEM 39.7 54.375  58.00 67.85 80.7 60.57917 10.646819 24       0

Side-by-side boxplots of Exam 1 score for five declared majors -- Business, Kinesiology, Nursing, Other, and STEM. Medians range narrowly from about 56 to 63 points with heavily overlapping interquartile ranges across all five groups; Kinesiology has two high outlier points near 84 and 86. — Figure 1:Exam 1 score by declared major area — EDA before the ANOVA test.

Five means, all sitting within about six points of each other (58.8 to 64.4), with heavily overlapping boxes — EDA alone doesn’t show an obvious winner or loser. That’s exactly the situation a formal test is for: is that six-point spread more than you’d expect from sampling variability alone, or is it a real difference?

42. Fitting and reading the ANOVA table¶

fit <- aov(exam_score ~ major_area, data = survey)
anova(fit)

Analysis of Variance Table

Response: exam_score
            Df  Sum Sq Mean Sq F value Pr(>F)
major_area   4   459.9 114.976  1.2055 0.3123
Residuals  115 10967.8  95.372

aov(y ~ g, data=) — same y ~ g, data= shape as t.test() in L11, just with a grouping variable that now has more than two levels — fits the model; anova() (note: no x in the name, unlike xchisq.test()) prints its table. Read each column:

Df (degrees of freedom): major_area’s row gets $k - 1 = 4$ (five groups); Residuals’ row gets $n - k = 115$ (120 students minus 5 groups).
Sum Sq (sum of squares): major_area’s 459.9 is variation between the five group means; Residuals’ 10967.8 is variation within groups (students scoring differently even inside the same major). ANOVA is literally comparing these two sources of variation.
Mean Sq: each Sum Sq divided by its own Df — an average squared variation per degree of freedom.
F value: Mean Sq(major_area) / Mean Sq(Residuals) = $114.976 / 95.372 = 1.2055$ . An $F$ near 1 means the between-group variation is about the same size as the natural within-group noise — no real evidence of a group effect. A large $F$ (well above 1) means the groups differ by more than noise alone would explain.
Pr(>F): the p-value, 0.3123 here.

$H_0$ : all five major means are equal ( $\mu_{Business} = \mu_{Kinesiology} = \cdots = \mu_{STEM}$ ); $H_a$ : at least one differs. With $p = 0.312 > \alpha = 0.05$ : fail to reject $H_0$ . This sample gives no evidence that mean Exam 1 score differs across declared major — the six-point spread in the means is small enough to be ordinary sampling variability, exactly what the near-1 $F$ -statistic already signaled.

53. Experimental-design vocabulary¶

ANOVA’s own name comes from designed experiments, not just observational survey data like survey. Four terms describe how such an experiment is built:

Term	Meaning
Factor	The variable a researcher deliberately manipulates (e.g., which fertilizer a plant receives).
Levels	The specific values/categories the factor takes (e.g., Control, Treatment A, Treatment B).
Randomization	Assigning experimental units to levels by chance, so no hidden variable systematically favors one level (e.g., shuffling which tree gets which fertilizer, instead of choosing by hand).
Replication	Measuring more than one unit per level, so you can tell a real effect from one unusually high or low measurement.

survey$major_area is not a designed-experiment factor in this sense — students chose their major, nobody randomly assigned it, so “major” is an observational grouping variable. That distinction matters for what you can claim: Section 2’s result (no evidence of a difference) describes an association question, and even a significant result there could never, by itself, support a claim that changing your major causes a score change. A true experiment, with random assignment to levels, is what lets you make that stronger causal claim.

64. A small worked experiment¶

Here’s what a genuine designed experiment looks like end to end, worked through with simulated data — 24 young almond trees, randomly assigned 8 each to three fertilizer levels (Control, Treatment A, Treatment B), yield in pounds measured once per tree at harvest (one replication per tree, 8 replications per level). This is clearly labeled classroom-simulation data, not a real Kern County yield record — the generating code is right here, reproducible with set.seed(2200):

set.seed(2200)
fert_yield <- data.frame(
  fertilizer = factor(rep(c("Control", "Treatment A", "Treatment B"), each = 8),
                       levels = c("Control", "Treatment A", "Treatment B")),
  yield_lb   = round(c(rnorm(8, mean = 28, sd = 3),
                        rnorm(8, mean = 33, sd = 3),
                        rnorm(8, mean = 31, sd = 3)), 1)
)
favstats(yield_lb ~ fertilizer, data = fert_yield)

   fertilizer  min     Q1 median     Q3  max    mean       sd n missing
1     Control 25.3 26.875  28.45 30.350 32.0 28.6375 2.485350 8       0
2 Treatment A 25.5 29.700  32.65 34.250 37.8 31.8000 4.314428 8       0
3 Treatment B 30.8 30.875  32.35 33.425 34.9 32.4250 1.603345 8       0

Side-by-side boxplots of simulated almond yield in pounds per tree for three fertilizer treatments, 8 simulated trees each. Control has the lowest median around 28.5 pounds; Treatment A and Treatment B have higher medians near 32 to 33 pounds, with Treatment A showing the widest spread of the three groups. — Figure 2:Simulated almond yield by fertilizer treatment — a small designed experiment.

Control’s mean (28.6 lb) looks visibly lower than either treatment’s (31.8 and 32.4 lb). Fit and read the ANOVA table exactly as in Section 2:

fert_fit <- aov(yield_lb ~ fertilizer, data = fert_yield)
anova(fert_fit)

Analysis of Variance Table

Response: yield_lb
           Df  Sum Sq Mean Sq F value  Pr(>F)  
fertilizer  2  65.966  32.983  3.6163 0.04471 *
Residuals  21 191.534   9.121                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$Df = 2$ (3 fertilizer levels $- 1$ ) and 21 (24 trees $- 3$ levels), $F = 3.6163$ , $p = 0.0447$ . This time $p < \alpha = 0.05$ : reject $H_0$ — this simulated experiment shows evidence that mean yield differs across the three fertilizer levels. R even flags it with a * under Signif. codes. Unlike Section 2’s major_area result, this conclusion — because fertilizer was randomly assigned — can support a genuinely causal claim: in this simulated trial, changing the fertilizer treatment plausibly changed the yield, not just correlated with it.

75. Which pair differs? Post-hoc comparisons with `TukeyHSD()`¶

Section 4’s ANOVA rejected $H_0$ : at $p = 0.0447$ , this sample gives evidence that mean yield differs across some of the three fertilizer levels. But the omnibus $F$ -test only ever asks one broad, pooled question — “is there a group effect anywhere?” It never says which group (or pair of groups) is driving that result. Answering that is a genuinely different, harder question, and it’s tempting to just run a separate t-test on every pair — Control vs. A, Control vs. B, A vs. B — but that shortcut has a real cost.

7.1Why not just run three t-tests? The multiple-comparisons problem¶

Every hypothesis test run at $\alpha = 0.05$ carries a 5% chance of a false positive (rejecting a true $H_0$ ) on that one test. Run several such tests on the same data and those 5%-chances stack: for $k$ groups there are $\binom{k}{2}$ pairs to compare — $\binom{3}{2} = 3$ pairs here, but $\binom{5}{2} = 10$ pairs for Section 1’s five majors — and the chance that at least one of those pairwise tests turns up a false positive, purely by chance, climbs well past 5% as the number of pairs grows. That inflated risk, across the whole family of comparisons, is the family-wise error rate, and it’s exactly what running uncorrected pairwise t-tests ignores.

TukeyHSD() (“Tukey’s Honest Significant Difference”) runs every pairwise comparison at once and adjusts both the p-values and the confidence intervals so the family-wise error rate — the chance of any false positive across all the pairs together — stays at the stated level (95% by default), no matter how many pairs there are. It takes the fitted aov object directly — no new model to fit:

TukeyHSD(fert_fit)


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

Fit: aov(formula = yield_lb ~ fertilizer, data = fert_yield)

$fertilizer
                          diff         lwr      upr     p adj
Treatment A-Control     3.1625 -0.64361423 6.968614 0.1152184
Treatment B-Control     3.7875 -0.01861423 7.593614 0.0512801
Treatment B-Treatment A 0.6250 -3.18111423 4.431114 0.9102904

mosaic loads the pipe operator %>%, so the same call also reads left-to-right, the same style as this book’s gf_* chains:

fert_fit %>% TukeyHSD()

produces the identical table above — TukeyHSD() is a base R function (in the stats package, always available, no library() needed), and piping into it just hands fert_fit in as its first argument, nothing new to learn.

7.2Reading the table¶

Each row is one pair, in the order (later level) - (earlier level):

diff: the difference in that pair’s sample means — e.g. Treatment A-Control = $31.8 - 28.6375 = 3.1625$ lb, matching Section 4’s favstats() means exactly.
lwr, upr: the 95% family-wise confidence interval for that pairwise difference — wider than an uncorrected pairwise CI would be, because the correction is spending some of its confidence on guarding all three intervals at once, not just one.
p adj: the p-value after the multiple-comparisons correction — compare this to $\alpha = 0.05$ , not an uncorrected pairwise p-value.

Tukey HSD 95 percent family-wise confidence intervals for the three pairwise differences in mean almond yield among fertilizer treatments. Treatment A minus Control spans about negative 0.6 to 7.0 pounds. Treatment B minus Control spans about negative 0.02 to 7.6 pounds, barely touching zero on the low end. Treatment B minus Treatment A spans about negative 3.2 to 4.4 pounds. All three horizontal interval lines cross the dashed vertical reference line at zero. — Figure 3:Tukey HSD pairwise confidence intervals for the fertilizer trial — every interval that crosses the dashed zero line is a pair not distinguishable at the 95% family-wise level.

A pair is declared to differ (at the family-wise 95% level) exactly when its interval excludes 0 — equivalently, when p adj $< 0.05$ . Here, none of the three intervals excludes 0, and no p adj drops below 0.05 — Treatment B-Control comes closest, at p adj = 0.051, with a lower bound of -0.019 lb, a hair’s width from excluding zero.

7.3An honest, real result: significant omnibus, no significant pair¶

This is not a contradiction, and it’s the single most important thing to take from this section: the omnibus $F$ -test in Section 4 pools all the data’s evidence into one broad question, while each post-hoc comparison answers a narrower question with only part of that evidence and pays a correction for being one of several such questions. With only $n = 8$ trees per fertilizer level, this experiment has enough combined evidence to reject “all three means are equal,” but not quite enough left over, once divided three ways and corrected, to pin down which specific pair crosses the family-wise line — even though Treatment B-Control (the largest gap) comes within a hair of it. The honest conclusion: this trial shows real evidence of some fertilizer effect, most likely driven by Treatment B outperforming Control, but this sample size doesn’t quite let Tukey’s correction confirm that specific pair at 95% family-wise confidence. A larger trial — more replications per level — is exactly what would sharpen that answer, since bigger samples shrink every confidence interval in the table above.

8Summary¶

aov(y ~ g, data=) fits a one-way ANOVA; anova(fit) prints the table — Df, Sum Sq, Mean Sq, F value, Pr(>F) — comparing variation between groups to variation within groups.
An $F$ -statistic near 1 signals no real group effect; a large $F$ with a small Pr(>F) signals the group means likely differ. $H_0$ : all group means equal; $H_a$ : at least one differs.
Factor = the manipulated variable; levels = its specific values; randomization = chance assignment of units to levels; replication = more than one measurement per level. All four appear together only in a true designed experiment, not in an observational grouping variable like major_area.
A significant ANOVA from a randomized experiment can support a causal claim; the same significant-or-not result from observational data (students choosing their own major) only ever describes an association.
A significant omnibus ANOVA only says at least one group differs — it never says which. Uncorrected pairwise t-tests inflate the family-wise error rate as the number of pairs grows; TukeyHSD(fit) compares every pair at once while holding that family-wise rate at the stated level (95% by default), reading a pair as “different” only when its adjusted confidence interval excludes 0 (equivalently, p adj $< 0.05$ ).
A significant ANOVA and a clean post-hoc winner don’t always arrive together: with small samples, TukeyHSD() can fail to confirm any single pair even after a significant omnibus test, exactly what happened with this lesson’s 8-tree-per-group fertilizer trial.
The full script that generated every figure and every number through Section 4 of this lesson, including the simulated fertilizer-trial generator, is committed at data/make_L13_figures.R — run it yourself to reproduce all of it exactly. Section 5’s TukeyHSD(fert_fit) reuses that same script’s fert_fit object with no new seed or data — run it right after Section 4’s code to reproduce the post-hoc numbers too.