17  ANOVA and the General Linear Model

The way statistics is usually taught tends to hide something important: ANOVA is not a different technique from regression. It is regression, the very same least-squares machinery from the last two chapters, applied to categorical predictors. Analysis of variance, t-tests, regression: they are all special cases of one thing, the General Linear Model (GLM). Once you see the unification, a large part of your statistics education collapses into a single idea.

The only new skill you need is coding: how to turn a category (“cat,” “dog”) into numbers a regression can chew on. That is the whole chapter.

17.1 Learning Objectives

  1. See that a categorical predictor is just a set of coded numbers in a regression.
  2. Interpret dummy-coded coefficients as differences from a reference group.
  3. Interpret effects-coded coefficients as differences from the grand mean.
  4. Recognize the ANOVA \(F\) as the same omnibus test regression already gave you.

17.2 A Categorical Predictor Is a Sliced-Up Continuous One

Let us start with a continuous relationship. Suppose some trait \(x\) makes people happy:

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(1925)  # Fisher's "Statistical Methods for Research Workers"
d <- tibble(x = sort(rnorm(100))) |>
  mutate(happy = x + rnorm(100))

Now imagine \(x\) measures something like “dogness,” where higher values mean more dog-like. We rarely measure dogness on a continuum; instead we lump animals into categories. Let us chop the continuum into four pet types:

d <- d |> mutate(pet = cut(x, 4, labels = c("cat", "fish", "pig", "dog")))

d |> count(pet)
pet n
cat 17
fish 47
pig 29
dog 7

Does pet type relate to happiness? A regression will happily tell us, even though pet is a factor:

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 ch16-pets. For example, book_data("ch16-pets") in R, Julia, or Python, or !bookdata name = "ch16-pets". in SPSS. Every language reads the same shipped files, so your numbers will match the ones printed here exactly.

lm_pet <- lm(happy ~ pet, data = d)
tidy2(lm_pet)
term estimate std.error statistic p.value
(Intercept) -1.32 0.24 -5.51 < .001
petfish 1.18 0.28 4.23 < .001
petpig 2.05 0.30 6.79 < .001
petdog 2.82 0.44 6.35 < .001 
INSERT FILE='data/sim/ch16-pets.sps'.
ONEWAY happy BY pet /STATISTICS DESCRIPTIVES.
                                  Descriptives
+-----------+---+-----+-----------+-------+-------------------+-------+-------+
|           |   |     |           |       |   95% Confidence  |       |       |
|           |   |     |           |       | Interval for Mean |       |       |
|           |   |     |           |       +---------+---------+       |       |
|           |   |     |    Std.   |  Std. |  Lower  |  Upper  |       |       |
|      pet  | N | Mean| Deviation | Error |  Bound  |  Bound  |Minimum|Maximum|
+-----------+---+-----+-----------+-------+---------+---------+-------+-------+
|happy cat  | 17|-1.32|        .91|    .22|    -1.79|     -.86|  -2.91|    .04|
|      dog  |  7| 1.50|       1.55|    .59|      .06|     2.93|   -.41|   4.14|
|      fish | 47| -.14|       1.05|    .15|     -.45|      .17|  -2.54|   1.89|
|      pig  | 29|  .73|        .76|    .14|      .44|     1.02|   -.53|   2.31|
|      Total|100|  .03|       1.25|    .13|     -.22|      .27|  -2.91|   4.14|
+-----------+---+-----+-----------+-------+---------+---------+-------+-------+

                             ANOVA
+--------------------+--------------+--+-----------+-----+----+
|                    |Sum of Squares|df|Mean Square|  F  |Sig.|
+--------------------+--------------+--+-----------+-----+----+
|happy Between Groups|         61.67| 3|      20.56|20.99|.000|
|      Within Groups |         94.01|96|        .98|     |    |
|      Total         |        155.69|99|           |     |    |
+--------------------+--------------+--+-----------+-----+----+
using CSV, DataFrames
df = CSV.read("data/sim/ch16-pets.csv", DataFrame)

using GLM; ftest(lm(@formula(happy ~ pet), df).model)
F-test against the null model:
F-statistic: 20.99 on 100 observations and 3 degrees of freedom, p-value: <1e-09
import pandas as pd
df = pd.read_csv("data/sim/ch16-pets.csv")

import statsmodels.formula.api as smf
smf.ols("happy ~ pet", data=df).fit().params.round(3)
Intercept     -1.323
pet[T.dog]     2.820
pet[T.fish]    1.185
pet[T.pig]     2.052
dtype: float64

R fit it without complaint. But look at those coefficient names - petfish, petpig, petdog. Where did they come from, and what do they mean? That is the coding question, and there is more than one answer.

17.3 Dummy Coding: Differences From a Reference

By default, R uses dummy coding (contr.treatment): it picks the first level (“cat”) as the reference, and each coefficient is the difference between that group’s mean and the reference group’s mean. Let us prove it by pulling the group means directly:

means <- d |>
  summarise(mean = mean(happy), .by = pet) |>
  arrange(pet)

means |> round2()
pet mean
cat -1.32
fish -0.14
pig 0.73
dog 1.50

Now rebuild every coefficient from those means:

coefs <- tidy(lm_pet) |> pull(estimate, name = term)

tibble(
  coefficient = c("(Intercept)", "petfish", "petdog"),
  from_model  = c(coefs[["(Intercept)"]], coefs[["petfish"]], coefs[["petdog"]]),
  from_means  = c(means$mean[1],                    # the reference (cat) mean
                  means$mean[2] - means$mean[1],    # fish mean - cat mean
                  means$mean[4] - means$mean[1])    # dog mean  - cat mean
) |>
  round2()
coefficient from_model from_means
(Intercept) -1.32 -1.32
petfish 1.18 1.18
petdog 2.82 2.82

There is no magic in the box. A dummy-coded coefficient is a difference between two group means, and the intercept is the reference group’s mean. That is exactly what a t-test between two groups reports, which is your first hint that t-tests, too, are just this same model.

17.4 Effects Coding: Differences From the Grand Mean

Dummy coding answers “how does each group differ from the reference?” Often we would rather ask “how does each group differ from the overall average?” For that we switch to effects coding (contr.sum):

d_e <- d |> mutate(pet_e = pet)
contrasts(d_e$pet_e) <- contr.sum(4)   # base R still owns contrast assignment

lm_pet_e <- lm(happy ~ pet_e, data = d_e)
tidy2(lm_pet_e)
term estimate std.error statistic p.value
(Intercept) 0.19 0.13 1.52 0.132
pet_e1 -1.51 0.21 -7.17 < .001
pet_e2 -0.33 0.16 -2.04 0.045
pet_e3 0.54 0.18 2.98 0.004

The coefficients changed, because the question changed. Now the intercept is the grand mean - and here is a subtlety worth stopping for: the grand mean of effects coding is the mean of the group means (each group weighted equally), not the mean of the raw scores. With unequal group sizes those two differ, and confusing them is a classic error.

GM <- mean(means$mean)      # mean of the group means (unweighted)
eff <- tidy(lm_pet_e) |> pull(estimate)

tibble(intercept = eff[1], grand_mean = GM) |> round2()
intercept grand_mean
0.19 0.19

Each effects coefficient is then that group’s mean minus the grand mean:

# first coefficient = cat mean - grand mean:
tibble(coef1 = eff[2], cat_minus_GM = means$mean[1] - GM) |> round2()
coef1 cat_minus_GM
-1.51 -1.51
ImportantWhere Did the Last Group Go?

Count the coefficients: with four groups, effects coding gives you an intercept plus only three effects. The fourth group (“dog”) seems to have vanished. It has not - it is implied. Because the effects are deviations from the grand mean, they must sum to zero, so the last group’s effect is minus the sum of the others:

tibble(
  implied_dog_effect = -sum(eff[2:4]),   # the three effects must sum to zero
  dog_minus_GM       = means$mean[4] - GM
) |>
  round2()
implied_dog_effect dog_minus_GM
1.31 1.31

They match. The “missing” group is never missing; it is the one the constraint solves for. (Work the algebra once and reference levels stop being confusing.)

NoteFrom the Mailbag (a real question from PSYC 754)

“I don’t see the difference between dummy coding and contrast coding. Why would labeling a 0 vs. -1 be a big change? … How is effects coding different from regular contrast coding?”

A perennial knot, so let’s cut it. Dummy coding uses 0s and 1s and compares each group to a single reference group (the one coded 0). Effects coding swaps that reference 0 for a \(-1\), and that small change quietly rewrites the question: each coefficient now compares a group to the grand mean rather than to a reference group. Contrast coding is the general family both belong to - any set of codes chosen to test a specific comparison - and effects coding is simply the member that yields deviations from the grand mean. So the “big change” from \(0\) to \(-1\) is really a change in what you compare against: one baseline group, or the average of them all.

17.5 The Punchline: This Is ANOVA

Whether you dummy-code or effects-code, the overall test - “do the group means differ at all?” - is identical, and it is identical to what analysis of variance calls the omnibus \(F\). Watch the same number fall out three ways:

tibble(
  source = c("1) the regression's own F test",
             "2) anova() on the regression object",
             "3) the classic aov() function"),
  F = c(
    glance(lm_pet)$statistic,
    anova(lm_pet) |> tidy() |> filter(term == "pet") |> pull(statistic),
    aov(happy ~ pet, data = d) |> tidy() |> filter(term == "pet") |> pull(statistic)
  )
) |>
  round2()
source F
1) the regression’s own F test 20.99
2) anova() on the regression object 20.99
3) the classic aov() function 20.99

One model, one \(F\), three names. ANOVA is regression with categorical predictors: the coding scheme only changes how we slice up and label the group differences, never the overall test. This is the General Linear Model: pick your predictors (continuous, categorical, or both), pick your coding, and it is all least squares underneath.

17.6 Two Categorical Predictors and Their Interaction

Add a second factor and you get a factorial design - still just regression. Suppose messiness also matters:

d <- d |>
  mutate(
    mess = factor(if_else(rnorm(100) > 0, "messy", "tidy")),
    sad  = -x + as.numeric(mess) + rnorm(100)
  )

lm_fac <- lm(sad ~ pet * mess, data = d)
tidy2(lm_fac)
term estimate std.error statistic p.value
(Intercept) 2.08 0.42 4.96 < .001
petfish -1.14 0.50 -2.27 0.026
petpig -1.61 0.50 -3.20 0.002
petdog -2.79 0.69 -4.02 < .001
messtidy 0.85 0.55 1.56 0.123
petfish:messtidy 0.61 0.64 0.95 0.343
petpig:messtidy 0.12 0.69 0.17 0.864
petdog:messtidy -0.40 1.01 -0.40 0.690 
INSERT FILE='data/sim/ch16-pets.sps'.

* GLM and UNIANOVA are the same procedure in SPSS; GLM is the name that also
* works in PSPP, so that is what the book runs. A factor has to be numeric
* here, which is what AUTORECODE does - it hands back the same groups with
* integer codes and the labels still attached.
AUTORECODE VARIABLES=pet mess /INTO petn messn.

GLM sad BY petn messn
  /DESIGN=petn messn petn*messn.
                  Tests of Between-Subjects Effects
+---------------+-----------------------+---+-----------+-----+----+
|               |Type III Sum Of Squares| df|Mean Square|  F  |Sig.|
+---------------+-----------------------+---+-----------+-----+----+
|Corrected Model|                  98.58|  7|      14.08|11.47|.000|
|petn           |                  53.15|  3|      17.72|14.42|.000|
|messn          |                  13.09|  1|      13.09|10.66|.002|
|petn × messn   |                   2.49|  3|        .83|  .67|.570|
|Error          |                 113.00| 92|       1.23|     |    |
|Total          |                 454.87|100|           |     |    |
|Corrected Total|                 211.58| 99|           |     |    |
+---------------+-----------------------+---+-----------+-----+----+
using CSV, DataFrames
df = CSV.read("data/sim/ch16-pets.csv", DataFrame)

using GLM; lm(@formula(sad ~ pet * mess), df)
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}

sad ~ 1 + pet + mess + pet & mess

Coefficients:
────────────────────────────────────────────────────────────────────────────────────
                            Coef.  Std. Error      t  Pr(>|t|)  Lower 95%  Upper 95%
────────────────────────────────────────────────────────────────────────────────────
(Intercept)              2.07932     0.418881   4.96    <1e-05   1.24739    2.91125
pet: dog                -2.78944     0.694635  -4.02    0.0001  -4.16905   -1.40984
pet: fish               -1.13909     0.50222   -2.27    0.0257  -2.13654   -0.141637
pet: pig                -1.60861     0.50222   -3.20    0.0019  -2.60606   -0.611159
mess: tidy               0.849402    0.546153   1.56    0.1233  -0.235306   1.93411
pet: dog & mess: tidy   -0.402716    1.00735   -0.40    0.6902  -2.4034     1.59796
pet: fish & mess: tidy   0.613691    0.643947   0.95    0.3431  -0.665244   1.89263
pet: pig & mess: tidy    0.11734     0.68522    0.17    0.8644  -1.24357    1.47825
────────────────────────────────────────────────────────────────────────────────────
import pandas as pd
df = pd.read_csv("data/sim/ch16-pets.csv")

import statsmodels.formula.api as smf
smf.ols("sad ~ pet * mess", data=df).fit().params.round(3)
Intercept                   2.079
pet[T.dog]                 -2.789
pet[T.fish]                -1.139
pet[T.pig]                 -1.609
mess[T.tidy]                0.849
pet[T.dog]:mess[T.tidy]    -0.403
pet[T.fish]:mess[T.tidy]    0.614
pet[T.pig]:mess[T.tidy]     0.117
dtype: float64

The pet:messmessy terms are the interaction, and the same logic applies here. An interaction coefficient is a difference of differences: how the effect of pet changes across levels of mess. It is the cross-product of the two predictors’ codes. A “significant interaction” is nothing to worry about; it just means the group-mean pattern is not purely additive.

17.7 Beyond the Straight Line: The Generalized Linear Model

Everything so far assumes a continuous outcome with roughly normal residuals. But plenty of outcomes are not like that - a yes/no, a count, a rare event. The generalized linear model extends the same framework by letting the outcome follow a different distribution through a link function. The most common case is a binary outcome via logistic regression:

# 0/1 outcome driven by x
d <- d |> mutate(passed = rbinom(100, size = 1, prob = plogis(x)))

glm_fit <- glm(passed ~ x, family = binomial, data = d)
tidy2(glm_fit)
term estimate std.error statistic p.value
(Intercept) -0.18 0.21 -0.86 0.392
x 0.73 0.26 2.79 0.005 
INSERT FILE='data/sim/ch16-pets.sps'.
LOGISTIC REGRESSION VARIABLES passed WITH x.
  Dependent Variable Encoding
+--------------+--------------+
|Original Value|Internal Value|
+--------------+--------------+
|.000000       |           .00|
|1.000000      |          1.00|
+--------------+--------------+

      Case Processing Summary
+--------------------+---+-------+
|Unweighted Cases    | N |Percent|
+--------------------+---+-------+
|Included in Analysis|100| 100.0%|
|Missing Cases       |  0|    .0%|
|Total               |100| 100.0%|
+--------------------+---+-------+

note: Estimation terminated at iteration number 4 because parameter estimates
changed by less than 0.001

                          Model Summary
+----+-----------------+--------------------+-------------------+
|Step|-2 Log likelihood|Cox & Snell R Square|Nagelkerke R Square|
+----+-----------------+--------------------+-------------------+
|1   |           129.26|                 .09|                .11|
+----+-----------------+--------------------+-------------------+

                      Classification Table
+--------------------------+-----------------------------------+
|                          |             Predicted             |
|                          +----------------+------------------+
|                          |     passed     |                  |
|                          +-------+--------+                  |
|       Observed           |.000000|1.000000|Percentage Correct|
+--------------------------+-------+--------+------------------+
|Step 1 passed   .000000   |     37|      16|             69.8%|
|                1.000000  |     24|      23|             48.9%|
|      --------------------+-------+--------+------------------+
|       Overall Percentage |       |        |             60.0%|
+--------------------------+-------+--------+------------------+

           Variables in the Equation
+---------------+----+----+----+--+----+------+
|               |  B |S.E.|Wald|df|Sig.|Exp(B)|
+---------------+----+----+----+--+----+------+
|Step 1 x       | .73| .26|7.77| 1|.005|  2.07|
|       Constant|-.18| .21| .73| 1|.392|   .83|
+---------------+----+----+----+--+----+------+
using CSV, DataFrames
df = CSV.read("data/sim/ch16-pets.csv", DataFrame)

using GLM; glm(@formula(passed ~ x), df, Binomial(), LogitLink())
StatsModels.TableRegressionModel{GeneralizedLinearModel{GLM.GlmResp{Vector{Float64}, Binomial{Float64}, LogitLink}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}

passed ~ 1 + x

Coefficients:
─────────────────────────────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)  Lower 95%  Upper 95%
─────────────────────────────────────────────────────────────────────────
(Intercept)  -0.18051     0.210769  -0.86    0.3918  -0.593609   0.232589
x             0.728377    0.261186   2.79    0.0053   0.216461   1.24029
─────────────────────────────────────────────────────────────────────────
import pandas as pd
df = pd.read_csv("data/sim/ch16-pets.csv")

import statsmodels.formula.api as smf
smf.logit("passed ~ x", data=df).fit().params.round(3)
Optimization terminated successfully.
         Current function value: 0.646319
         Iterations 5
Intercept   -0.181
x            0.728
dtype: float64

Same formula syntax, same partialling logic, same goal of quantifying the uncertainty. We have just swapped lm() for glm() and told it the outcome is binomial. Linear model, general linear model, generalized linear model: it is the same core idea extended step by step.

17.8 Challenge

TipDo One Yourself

Using the pet factor and happy from this chapter:

  1. Fit the model with dummy coding and, without reading the coefficients off the summary, reproduce all four coefficients (intercept + three) from the group means.
  2. Refit with effects coding and reproduce the intercept (grand mean) and all four group effects - including the implied “dog” effect - from the group means.
  3. Confirm the omnibus \(F\) is identical under both coding schemes, and identical to summary(aov(happy ~ pet)).
  4. In one sentence each, state what the intercept means under dummy coding versus effects coding. If you can do that cleanly, you understand why “ANOVA vs. regression” was always a false choice.

17.9 Where We Go Next

You now hold the unifying idea of the whole enterprise: the General Linear Model. Correlation, t-tests, regression, ANOVA, ANCOVA - all one model wearing different outfits. We build these coding schemes by hand in the dedicated coding chapter; for the definitive treatment, see Cohen, Cohen, West, and Aiken (Aiken 2002) and the accessible comparison of dummy and effects coding in Alkharusi (Alkharusi 2012); for analysis of variance specifically, Iversen and Norpoth (Iversen and Norpoth 1987) remains a compact classic. From here, the road leads into the specific ANOVA designs (one-way, factorial, repeated measures), each of which is just another coded regression.

Aiken, Stephen G. West, Patricia Cohen. 2002. Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences. 3rd ed. Routledge. https://doi.org/10.4324/9780203774441.
Alkharusi, Hussain. 2012. “Categorical Variables in Regression Analysis: A Comparison of Dummy and Effect Coding.” International Journal of Education 4 (April): 202–10. https://doi.org/10.5296/ije.v4i2.1962.
Iversen, Gudmund R., and Helmut Norpoth. 1987. Analysis of Variance. 2nd ed. Sage University Papers. Quantitative Applications in the Social Sciences, no. 07-001. Sage Publications.