Analysis of Variance
analysis of variance compares the means of groups. it was born as its own subject, with its own tables and vocabulary, but once a categorical predictor is written as a block of indicator columns it is simply a linear gaussian model — the same \(\hat\beta = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}\), the same nested \(F\)-test from estimation & inference. what is genuinely new is the bookkeeping: how to code a factor so the design matrix has full rank, what the coefficients then mean, and how the total variation splits into a between-groups and a within-groups part. 𐃏
the same machinery then carries covariates alongside factors (ancova), several factors crossed or nested (general linear models), and product and power columns that bend the additive linear model (interactions, polynomials).
factors as columns
a factor is a categorical predictor with \(J\) levels; it enters the design matrix through dummy variables, indicators of the levels. in a completely randomised experiment units are randomly allocated to the \(J\) groups, and the responses \(Y_{j1}, \dots, Y_{jn_j}\) at level \(j\) are replicates.
with \(K\) replicates per group (\(N = JK\)) there are three natural ways to write the same model.
- (1) cell means, \(\mathbb{E}(Y_{jk}) = \mu_j\): \(\mathbf{X}\) has one indicator column per level, \(\mathbf{X}^\top\mathbf{X} = K\mathbf{I}_J\), and \(\hat\mu_j = \bar Y_{j\cdot}\).
- (2) overparameterised, \(\mathbb{E}(Y_{jk}) = \mu + \alpha_j\): the intercept column is the sum of the \(J\) indicators, so \(\mathbf{X}^\top\mathbf{X}\) is singular and the normal equations have the one-parameter family of solutions \(\hat\mu = \lambda\), \(\hat\alpha_j = \bar Y_{j\cdot} - \lambda\). the sum-to-zero constraint \(\sum_j\alpha_j = 0\) picks \(\hat\mu = \bar Y_{\cdot\cdot}\), \(\hat\alpha_j = \bar Y_{j\cdot} - \bar Y_{\cdot\cdot}\).
- (3) corner point, \(\alpha_1 = 0\): drop the first indicator; then \(\hat\mu = \bar Y_{1\cdot}\) and \(\hat\alpha_j = \bar Y_{j\cdot} - \bar Y_{1\cdot}\), differences from a reference level.
all three have the same column space, hence the same fitted values (the group means) and the same RSS; they differ only in what a coefficient means.
R’s default is (3), with the first factor level as reference — so every factor coefficient in a summary() is a contrast with that level, and its \(t\)-test asks only whether that level differs from the reference. cell means (1) is the cleanest for one factor but has no overall level to which a second factor’s effects could be added; that is why the constrained forms dominate.
genetically similar seeds were randomly raised under control conditions or one of two enriched treatments, \(K = 10\) plants each, \(N = 30\), and the plants dried and weighed:
| group | weights (plants 1–5) |
|---|---|
| control | 4.17 5.58 5.18 6.11 4.50 |
| treatment A | 4.81 4.17 4.41 3.59 5.87 |
| treatment B | 6.31 5.12 5.54 5.50 5.37 |
| group | weights (plants 6–10) | total | mean |
|---|---|---|---|
| control | 4.61 5.17 4.53 5.33 5.14 | 50.32 | 5.032 |
| treatment A | 3.83 6.03 4.89 4.32 4.69 | 46.61 | 4.661 |
| treatment B | 5.29 4.92 6.15 5.80 5.26 | 55.26 | 5.526 |
grand total 152.19, grand mean 5.073. the three parameterisations give
- cell means: \((5.032, 4.661, 5.526)\);
- sum-to-zero: \(\hat\mu = 5.073\), \(\hat\alpha = (-0.041, -0.412, 0.453)\), summing to zero;
- corner point: \(\hat\mu = 5.032\), \(\hat\alpha_A = -0.371\), \(\hat\alpha_B = 0.494\) — exactly what
lm(weight ~ group)prints.
proofs
(1) column \(j\) is \(K\) ones on group \(j\): orthogonal columns, \(\mathbf{X}^\top\mathbf{X} = K\mathbf{I}\), \(\mathbf{X}^\top\mathbf{y} = (Y_{1\cdot},\dots,Y_{J\cdot})^\top\), so \(\hat\mu_j = Y_{j\cdot}/K\).
(2) the normal equations are
\begin{align} N\mu + K\textstyle\sum_j\alpha_j &= Y_{\cdot\cdot} && \text{intercept row}\\ K\mu + K\alpha_j &= Y_{j\cdot},\quad j = 1,\dots,J && \text{indicator rows} \end{align}
and the first is the sum of the others: \(J\) equations, \(J+1\) unknowns, so \(\alpha_j = Y_{j\cdot}/K - \lambda\) with \(\mu = \lambda\) free. imposing \(\sum_j\alpha_j = 0\): \(\sum_j Y_{j\cdot}/K = J\lambda\), so \(\lambda = Y_{\cdot\cdot}/(JK) = \bar Y_{\cdot\cdot}\).
(3) drop \(\alpha_1\): the indicator rows give \(\alpha_j = Y_{j\cdot}/K - \mu\) for \(j \ge 2\); substituting into the intercept row, \((N - (J-1)K)\mu = Y_{\cdot\cdot} - \sum_{j\ge2}Y_{j\cdot} = Y_{1\cdot}\), and \(N - (J-1)K = K\).
the variance decomposition and the f-test
the hypothesis “all group means are equal” is the intercept-only model nested inside the factor model, so the nested f-test applies with \(p = J\) and \(q = J - 1\). what earns the name analysis of variance is that the two sums of squares are an orthogonal split of the total.
\[ \underbrace{\sum_{j,k}(Y_{jk} - \bar Y_{\cdot\cdot})^2}_{\mathrm{TSS}} = \underbrace{\sum_{j}n_j(\bar Y_{j\cdot} - \bar Y_{\cdot\cdot})^2}_{\text{between}} + \underbrace{\sum_{j,k}(Y_{jk} - \bar Y_{j\cdot})^2}_{\text{within} \,=\, \mathrm{RSS}}, \qquad F = \frac{\text{between}/(J-1)}{\text{within}/(N-J)} \sim F_{J-1,\,N-J} \] under \(H_0: \mu_1 = \dots = \mu_J\) with gaussian errors of constant variance.
the within-groups mean square \(\mathrm{RSS}/(N-J)\) estimates \(\sigma^2\) whatever the means; the between-groups mean square estimates \(\sigma^2\) only if the means are equal and is inflated otherwise. \(F\) is their ratio.
\(\sum Y_{jk}^2 = 786.318\) and \(Y_{\cdot\cdot}^2/N = 152.19^2/30 = 772.060\), so \(\mathrm{TSS} = 14.258\). within: \[ 786.318 - 10(5.032^2 + 4.661^2 + 5.526^2) = 786.318 - 775.826 = 10.492\,. \] between: \(10(0.041^2 + 0.412^2 + 0.453^2) = 3.766\).
| source | df | SS | MS | F |
|---|---|---|---|---|
| between | 2 | 3.766 | 1.883 | 4.85 |
| within | 27 | 10.492 | 0.3886 | |
| total | 29 | 14.258 |
\(F = 4.85 > F_{0.95;2,27} = 3.35\), \(p = 0.016\): the mean dried weight differs between at least two groups. the rest of lm’s summary is the same arithmetic: \(\hat\sigma = \sqrt{0.3886} = 0.6234\), \(R^2 = 3.766/14.258 = 0.264\); the standard error of each treatment-minus-control coefficient is \(\hat\sigma\sqrt{2/K} = 0.2788\), and of the control mean \(\hat\sigma/\sqrt K = 0.1971\).
the trap in the printout: neither treatment differs significantly from control (\(p = 0.19\) for A, \(0.088\) for B), yet the \(F\)-test rejects. the coefficient \(t\)-tests are contrasts with the reference level; the largest contrast, B against A, is \(0.865\) with standard error \(0.2788\), \(t = 3.1\). “do the groups differ?” is the \(F\) question, not the \(t\) questions.
proofs
write \(Y_{jk} - \bar Y_{\cdot\cdot} = (Y_{jk} - \bar Y_{j\cdot}) + (\bar Y_{j\cdot} - \bar Y_{\cdot\cdot})\) and square:
\begin{align} \mathrm{TSS} &= \sum_{j,k}(Y_{jk} - \bar Y_{j\cdot})^2 + \sum_j n_j(\bar Y_{j\cdot} - \bar Y_{\cdot\cdot})^2 + 2\sum_j(\bar Y_{j\cdot} - \bar Y_{\cdot\cdot})\underbrace{\sum_k(Y_{jk} - \bar Y_{j\cdot})}_{=\,0} \end{align}
the within term is the RSS of the factor model (fitted values are group means), and TSS is the RSS of the intercept-only model, so between \(= \mathrm{RSS}_0 - \mathrm{RSS}_1\) and \(F\) is the nested \(F\) statistic with \(q = J - 1\), \(N - p = N - J\).
ancova: comparing adjusted means
models mixing factor indicators with continuous covariates; the group means are compared after adjustment for the covariates. with one factor and one covariate and a common slope, \(\mathbb{E}(Y_{jk}) = \mu + \alpha_j + \gamma x_{jk}\) — parallel lines, one per group.
in the common-slope model with corner-point coding, each fitted line passes through its group’s mean point \((\bar x_j, \bar y_j)\), and so \[ \hat\alpha_j = (\bar y_j - \bar y_1) - \hat\gamma\,(\bar x_j - \bar x_1): \] the adjusted difference is the raw difference minus the part explained by the groups’ different covariate averages.
on the training-methods data (fit and coefficient table on estimation & inference), the raw mean scores are A 4.43, B 7.57, C 6.71 and the mean aptitudes A 2.14, B 3.43, C 2.71. with \(\hat\gamma = 0.743\), B’s raw advantage of 3.14 becomes \(3.14 - 0.743 \times 1.29 = 2.19\) and C’s 2.29 becomes \(2.29 - 0.743 \times 0.57 = 1.86\): part of B’s lead was its abler intake. this is the whole point of ancova — a randomised experiment balances covariates on average, but a given sample does not.
proofs
the intercept and the indicators together span every group’s indicator, so \(\mathbf{X}^\top\mathbf{e} = \mathbf{0}\) forces the residuals to sum to zero within each group: \(\bar y_j = \hat\mu + \hat\alpha_j + \hat\gamma\bar x_j\) (with \(\hat\alpha_1 = 0\)). subtract the \(j = 1\) equation.
general linear models: crossed and nested factors
a gaussian linear model with any mix of factors and covariates. two factors are crossed if every combination of levels is observed (the two-factor anova: a second block of indicators, plus an interaction block if needed), and nested if the levels of one occur only within a single level of the other.
drug \(A_1\) is trialled in hospitals \(B_1, B_2, B_3\) and drug \(A_2\) in \(B_4, B_5\). the saturated model \(\mathbb{E}(Y_{jkl}) = \mu + \alpha_j + (\alpha\beta)_{jk}\) has eight symbols \(\mu\), \(\alpha_1\), \(\alpha_2\), \((\alpha\beta)_{11}\), \((\alpha\beta)_{12}\), \((\alpha\beta)_{13}\), \((\alpha\beta)_{24}\), \((\alpha\beta)_{25}\) and three corner-point constraints \(\alpha_1 = (\alpha\beta)_{11} = (\alpha\beta)_{24} = 0\), leaving five free parameters — one per hospital, as a saturated model must. reading them: \(\mu\) is the \(B_1\) mean; \((\alpha\beta)_{12}, (\alpha\beta)_{13}\) compare \(B_2, B_3\) with \(B_1\); \((\alpha\beta)_{25}\) compares \(B_5\) with \(B_4\); and \(\alpha_2\) compares \(B_4\) with \(B_1\), confounding the drug with the hospital baseline. there is no “\(B_1\) on drug \(A_2\)” cell, so no crossed drug-by-hospital interaction is estimable.
every model on this page inherits the gaussian assumptions — normal, linear, constant variance, independent — and they are checked through the residuals exactly as for any regression (regression diagnostics). with unbalanced data, factor sums of squares become order-dependent; fit the competing orders, or compare the nested models you actually mean.
interactions and non-linear terms
the additive assumption is that the effect of changing \(x_j\) does not depend on the other predictors. an interaction relaxes it with a product column: \[ Y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_1x_2 + \varepsilon = \beta_0 + (\beta_1 + \beta_3x_2)\,x_1 + \beta_2x_2 + \varepsilon, \] so the slope in \(x_1\) is \(\beta_1 + \beta_3x_2\), changing by \(\beta_3\) per unit of \(x_2\) (and symmetrically). the model is still linear in \(\beta\).
the hierarchical principle: if an interaction is in the model, its main effects stay in, significant or not. without \(\beta_1x_1\) the model forces the \(x_1\)-slope to vanish at \(x_2 = 0\) — an arbitrary constraint that moves when \(x_2\) changes units. interactions combine factors with factors, covariates with covariates, or one of each; a factor-by-covariate interaction gives each group its own slope, the non-parallel generalisation of ancova.
in ISLR’s Credit data, Balance on Income (thousands of dollars) and Student:
- additive: \[ \widehat{\text{Balance}} = 211.1 + 5.98\,\text{Income} + 382.7\,\mathbb{I}(\text{student}) \] — parallel lines 383 dollars apart at every income;
- with
Income * Student: \[ 200.6 + 6.22\,\text{Income} + 476.7\,\mathbb{I}(\text{student}) - 2.00\,\text{Income}\times\mathbb{I}(\text{student})\,. \] non-students get \(200.6 + 6.22\,\text{Income}\), students \(677.3 + 4.22\,\text{Income}\): a higher start but a flatter slope, so the lines converge as income rises.
the interaction has \(t = -1.15\), \(p = 0.25\): the data do not demand it, and the additive model stands. the example shows what an interaction would mean, not that this one is real.
polynomial and fractional-polynomial terms
a curved residual plot suggests a curved mean, and the quadratic \(\mathbb{E}(Y_i) = \beta_0 + \beta_1x_i + \beta_2x_i^2\) is still linear regression on the columns \(x\) and \(x^2\). centre and scale first, \(\tilde x_i = (x_i - \bar x)/\operatorname{sd}(x)\): raw \(x\) and \(x^2\) are highly collinear and centring removes most of it; the intercept becomes the mean response at the average \(x\) rather than at a possibly impossible \(x = 0\); and slopes read per standard deviation.
the quadratic is symmetric about its turning point. fractional polynomials \(\mathbb{E}(Y_i) = \beta_0 + \beta_1x_i^{q}\), with \(\log x_i\) at \(q = 0\), buy asymmetric shapes for one extra column: try \(q \in \{-2, -1, -\tfrac12, 0, \tfrac12, 1, 2, 3\}\) and keep the best fit. negative powers fall steeply then flatten, \(\log x\) and \(\sqrt x\) rise and flatten, \(x^2\) and \(x^3\) accelerate. beyond fixed powers lie piecewise polynomials — splines and smoothing.
see also
- estimation & inference — the nested \(F\)-test this page specialises, and the worked ancova fit
- regression diagnostics — checking the gaussian assumptions every anova inherits; collinearity and vif
- generalised linear models — anova as one cell of the response-scale taxonomy; log-linear models are its count analogue
- model selection — choosing which factors, interactions and powers to keep
- splines and smoothing — flexible alternatives to polynomial terms
- additive models — the additive assumption, kept, with each term made non-linear
Comments