18  Basic ANOVA

We closed the last chapter with a slogan: ANOVA is regression with categorical predictors. That is true, and it is the deep view. But ANOVA also has a classical view, one built entirely out of the variance ideas from early in the book, and seeing that view makes the machinery click. Analysis of variance earns its name honestly: it works by partitioning variance into the part explained by group membership and the part left over. Let’s take it apart.

18.1 Learning Objectives

  1. Partition the total sum of squares into between-group and within-group pieces.
  2. Build the ANOVA \(F\) ratio from those pieces and understand what it compares.
  3. Confirm that the classical ANOVA and the regression GLM give the identical test.
  4. Follow a significant omnibus test with appropriate post-hoc comparisons.

18.2 The Question and the Data

One-way ANOVA asks: do three or more groups differ on a continuous outcome? Say we test a plant fertilizer at three doses:

library(tidyverse)    # dplyr, ggplot2, purrr, tibble, readr, stringr, forcats
library(broom)        # tidy(), glance(), augment() for model output
source("_common.R")   # book-wide helpers: round2(), fmt_p(), tidy2()

set.seed(17)
d <- tibble(
  dose   = factor(rep(c("none", "low", "high"), each = 30),
                  levels = c("none", "low", "high")),
  growth = rnorm(90, mean = rep(c(10, 12, 15), each = 30), sd = 3)
)

d |> summarise(mean = mean(growth), .by = dose) |> round2()
dose mean
none 11.04
low 11.85
high 14.24

The group means differ in our sample - but is that real, or just sampling noise? ANOVA answers by comparing two kinds of variation.

18.3 Partitioning the Variance

Here is the whole idea in one line: the total spread of every observation around the grand mean can be split cleanly into spread between groups and spread within groups.

\[SS_{total} = SS_{between} + SS_{within}\]

  • \(SS_{total}\): how far each observation is from the overall (grand) mean.
  • \(SS_{between}\): how far each group mean is from the grand mean - the variation our grouping explains (the signal).
  • \(SS_{within}\): how far each observation is from its own group’s mean - the leftover variation (the noise).

Let’s build all three by hand and confirm they add up:

grand <- mean(d$growth)

groups <- d |> summarise(n = n(), group_mean = mean(growth), .by = dose)

SS_total   <- d |> summarise(s = sum((growth - grand)^2)) |> pull(s)
SS_between <- groups |> summarise(s = sum(n * (group_mean - grand)^2)) |> pull(s)
SS_within  <- d |>
  left_join(groups, by = "dose") |>
  summarise(s = sum((growth - group_mean)^2)) |>
  pull(s)

tibble(SS_total            = SS_total,
       between_plus_within = SS_between + SS_within) |>
  round2()   # they match
SS_total between_plus_within
1060.5 1060.5

The total variance really does split exactly into “explained by group” plus “unexplained.” That partition is analysis of variance.

18.4 The F Ratio

We cannot compare sums of squares directly - they depend on how many things went into each. We convert each to a mean square by dividing by its degrees of freedom, then take their ratio:

\[F = \frac{MS_{between}}{MS_{within}} = \frac{SS_{between}/(k-1)}{SS_{within}/(N-k)}\]

with \(k\) groups and \(N\) observations. \(F\) asks a simple question: is the variation between groups bigger than the variation within them? If groups matter, the between-group spread is large relative to the within-group noise and \(F\) climbs well above 1.

k <- 3; N <- 90

tibble(
  MS_between = SS_between / (k - 1),
  MS_within  = SS_within  / (N - k)
) |>
  mutate(
    F = MS_between / MS_within,
    p = pf(F, df1 = k - 1, df2 = N - k, lower.tail = FALSE)
  ) |>
  mutate(p = fmt_p(p)) |>
  round2()
MS_between MS_within F p
82.86 10.28 8.06 < .001

And of course, R will do the bookkeeping for us with aov() - confirm it matches our by-hand work:

NoteWorking in SPSS, Julia, or Python?

The code tabs below assume this chapter’s data is already loaded. Grab the one-file setup for your language from Getting the Book’s Data, run it once, then load what you need by name - this chapter uses ch17-dose. For example, book_data("ch17-dose") in R, Julia, or Python, or !bookdata name = "ch17-dose". in SPSS. Every language reads the same shipped files, so your numbers will match the ones printed here exactly.

aov(growth ~ dose, data = d) |>
  tidy() |>
  mutate(p.value = fmt_p(p.value)) |>
  round2()
term df sumsq meansq statistic p.value
dose 2 165.73 82.86 8.06 < .001
Residuals 87 894.78 10.28 NA NA
INSERT FILE='data/sim/ch17-dose.sps'.
ONEWAY growth BY dose /STATISTICS DESCRIPTIVES /POSTHOC=TUKEY.
                                  Descriptives
+------------+--+-----+-----------+-------+-------------------+-------+-------+
|            |  |     |           |       |   95% Confidence  |       |       |
|            |  |     |           |       | Interval for Mean |       |       |
|            |  |     |           |       +---------+---------+       |       |
|            |  |     |    Std.   |  Std. |  Lower  |  Upper  |       |       |
|       dose | N| Mean| Deviation | Error |  Bound  |  Bound  |Minimum|Maximum|
+------------+--+-----+-----------+-------+---------+---------+-------+-------+
|growth high |30|14.24|       3.38|    .62|    12.98|    15.50|   5.04|  20.46|
|       low  |30|11.85|       3.76|    .69|    10.45|    13.26|   1.33|  19.33|
|       none |30|11.04|       2.30|    .42|    10.18|    11.90|   6.95|  15.15|
|       Total|90|12.38|       3.45|    .36|    11.66|    13.10|   1.33|  20.46|
+------------+--+-----+-----------+-------+---------+---------+-------+-------+

                             ANOVA
+---------------------+--------------+--+-----------+----+----+
|                     |Sum of Squares|df|Mean Square|  F |Sig.|
+---------------------+--------------+--+-----------+----+----+
|growth Between Groups|        165.73| 2|      82.86|8.06|.001|
|       Within Groups |        894.78|87|      10.28|    |    |
|       Total         |       1060.50|89|           |    |    |
+---------------------+--------------+--+-----------+----+----+

                         Multiple Comparisons (growth)
+--------------------------+-----------------+--------+----+------------------+
|                          |                 |        |    |  95% Confidence  |
|                          |                 |        |    |     Interval     |
|                          |                 |        |    +---------+--------+
|        (I)       (J)     | Mean Difference |  Std.  |    |  Lower  |  Upper |
|        Family    Family  |     (I - J)     |  Error |Sig.|  Bound  |  Bound |
+--------------------------+-----------------+--------+----+---------+--------+
|Tukey   high      low     |             2.39|     .83|.014|      .41|    4.36|
|HSD               none    |             3.20|     .83|.001|     1.22|    5.17|
|       -------------------+-----------------+--------+----+---------+--------+
|        low       high    |            -2.39|     .83|.014|    -4.36|    -.41|
|                  none    |              .81|     .83|.594|    -1.17|    2.78|
|       -------------------+-----------------+--------+----+---------+--------+
|        none      high    |            -3.20|     .83|.001|    -5.17|   -1.22|
|                  low     |             -.81|     .83|.594|    -2.78|    1.17|
+--------------------------+-----------------+--------+----+---------+--------+
using CSV, DataFrames
df = CSV.read("data/sim/ch17-dose.csv", DataFrame)

using GLM
ftest(lm(@formula(growth ~ dose), df).model)
F-test against the null model:
F-statistic: 8.06 on 90 observations and 2 degrees of freedom, p-value: 0.0006
import pandas as pd
df = pd.read_csv("data/sim/ch17-dose.csv")

import statsmodels.formula.api as smf
from statsmodels.stats.anova import anova_lm
anova_lm(smf.ols("growth ~ dose", data=df).fit())
            df      sum_sq    mean_sq         F    PR(>F)
dose       2.0  165.727050  82.863525  8.056913  0.000616
Residual  87.0  894.775284  10.284773       NaN       NaN

Same \(F\), same p, exactly as we computed by hand.

18.5 It Really Is the Same as Regression

To honor the slogan from the last chapter, let’s prove the classical ANOVA and the regression GLM are one and the same by pulling the identical \(F\) out of lm():

fit <- lm(growth ~ dose, data = d)

# the regression's F test == the ANOVA
anova(fit) |> tidy() |> mutate(p.value = fmt_p(p.value)) |> round2()
term df sumsq meansq statistic p.value
dose 2 165.73 82.86 8.06 < .001
Residuals 87 894.78 10.28 NA NA
tibble(
  regression_F = glance(fit)$statistic,
  aov_F = aov(growth ~ dose, data = d) |> tidy() |> slice(1) |> pull(statistic)
) |>
  round2()
regression_F aov_F
8.06 8.06
INSERT FILE='data/sim/ch17-dose.sps'.
REGRESSION /STATISTICS COEFF R ANOVA /DEPENDENT growth /METHOD=ENTER dose.
/tmp/RtmpKQhPo7/file34f06e064a32.sps:2.70-2.73: warning: REGRESSION: dose is
not a numeric variable.  It will not be included in the variable list.
    2 | REGRESSION /STATISTICS COEFF R ANOVA /DEPENDENT growth /METHOD=ENTER
dose.
      |
^~~~
using CSV, DataFrames
df = CSV.read("data/sim/ch17-dose.csv", DataFrame)

using GLM
coeftable(lm(@formula(growth ~ dose), df))
────────────────────────────────────────────────────────────────────────
                Coef.  Std. Error      t  Pr(>|t|)  Lower 95%  Upper 95%
────────────────────────────────────────────────────────────────────────
(Intercept)  14.2399     0.585513  24.32    <1e-39   13.0761   15.4036
dose: low    -2.3878     0.828041  -2.88    0.0050   -4.03362  -0.741977
dose: none   -3.19644    0.828041  -3.86    0.0002   -4.84226  -1.55062
────────────────────────────────────────────────────────────────────────
import pandas as pd
df = pd.read_csv("data/sim/ch17-dose.csv")

import statsmodels.formula.api as smf
from statsmodels.stats.anova import anova_lm
fit = smf.ols("growth ~ dose", data=df).fit()
anova_lm(fit)                  # the regression's F test == the ANOVA
fit.fvalue
            df      sum_sq    mean_sq         F    PR(>F)
dose       2.0  165.727050  82.863525  8.056913  0.000616
Residual  87.0  894.775284  10.284773       NaN       NaN
np.float64(8.056913085287713)

Identical. The “two techniques” were always one. ANOVA is just the variance-partitioning story told about a regression on categorical predictors.

18.6 After a Significant F: Post-Hoc Comparisons

The omnibus \(F\) only says “somewhere among these groups there is a difference.” It does not say which groups differ. For that we run post-hoc comparisons - and because we are now making several comparisons at once, we must control the inflated false-alarm rate (recall Type I error from the hypothesis-testing chapter). Tukey’s Honest Significant Difference does exactly that:

TukeyHSD(aov(growth ~ dose, data = d)) |>
  tidy() |>
  mutate(adj.p.value = fmt_p(adj.p.value)) |>
  round2()
term contrast null.value estimate conf.low conf.high adj.p.value
dose low-none 0 0.81 -1.17 2.78 0.594
dose high-none 0 3.20 1.22 5.17 < .001
dose high-low 0 2.39 0.41 4.36 0.014

Each row is a pairwise comparison with a confidence interval already adjusted for the fact that we peeked at all of them. Any interval that excludes zero is a real difference; any that includes zero is not.

ImportantDo Not Skip Straight to the Pairwise Tests

The two-step procedure, omnibus \(F\) first and post-hoc comparisons only if it is significant, is not just a formality. Running every pairwise t-test without the omnibus gate (and without correction) is how you manufacture false positives, one comparison at a time. The \(F\) is the gatekeeper, and Tukey (or a similar correction) keeps the gate honest.

18.7 Challenge

TipDo One Yourself
  1. Simulate four groups on a continuous outcome. Compute \(SS_{total}\), \(SS_{between}\), and \(SS_{within}\) by hand, confirm they sum, and build \(F\) from the mean squares. Check against summary(aov(...)).
  2. Confirm the same \(F\) falls out of anova(lm(y ~ group)).
  3. Make the group means nearly equal, then very different, keeping the within-group SD fixed. Watch what happens to \(F\) - and explain it in terms of the between-versus-within comparison.

18.8 Where We Go Next

You now hold the thread that runs through all of it: quantify your uncertainty, and remember that variance is the currency you quantify it in. From the first histogram to this \(F\) ratio, it was one idea the whole way down. But there is a loose end worth pulling. Every ANOVA you just ran handed R a factor and trusted it to build the right columns of numbers, and how those columns are built quietly decides what each coefficient means. That is the subject of the next chapter, where we stop trusting the defaults and take the coding into our own hands.