Additive Models

an additive model replaces each linear term \(\beta_j x_j\) of linear regression with an unknown smooth function \(f_j(x_j)\), and adds them up. a generalised additive model (gam) does the same inside the link of a generalised linear model. the single-predictor smoother is the building block and lives on splines and smoothing; this page is about what changes when there are several smooths in one model: identifiability, one penalty per term, choosing several smoothing parameters at once, and the non-gaussian fit. 𐃏

the recipe in one line: give each smooth a basis and a penalty, fit by penalised least squares (additive model) or penalised irls (gam), and pick the penalties by generalised cross-validation.

the additive model#

additivity is a strong assumption: \(f_1(x) + f_2(v)\) is a very special surface. the interaction \(xv\) is not of this form, since every additive surface has \(\partial^2 / \partial x \partial v = 0\) while \(xv\) has mixed derivative \(1\).

basis and penalty for each term#

the notes use the most transparent basis: piecewise linear “tent” functions on \(k_1\) evenly spaced knots \(x_1^* < \dots < x_{k_1}^*\) with spacing \(h\), \(b_j(x) = \max(0, 1 - |x - x_j^*|/h)\), so that \[ f_1(x) = \sum_{j=1}^{k_1} b_j(x)\, \delta_j, \qquad f_2(v) = \sum_{j=1}^{k_2} \mathcal{B}_j(v)\, \gamma_j . \] two properties carry everything below. interpolation: \(b_i(x_j^*) = 1\) if \(i = j\) and \(0\) otherwise, so \(\delta_j = f_1(x_j^*)\) and \(f_1\) is the broken line through the points \((x_j^*, \delta_j)\). partition of unity: \(\sum_j b_j(x) = 1\) on the knot range. at the data, \(\mathbf{f}_1 = \mathbf{X}_1 \boldsymbol\delta\) with \([\mathbf{X}_1]_{ij} = b_j(x_i)\). 𐃏

the roughness penalty is the sum of squared second differences of the knot values, \[ \boldsymbol\delta^\top \mathbf{S}_1 \boldsymbol\delta = \|\mathbf{D}_1 \boldsymbol\delta\|^2 = \sum_{j=2}^{k_1-1} (\delta_{j-1} - 2\delta_j + \delta_{j+1})^2\,, \] with \(\mathbf{D}_1\) the \((k_1 - 2) \times k_1\) second-difference matrix. it is zero exactly when \(f_1\) is a straight line. any other basis–penalty pair (cubic regression splines, thin plate splines; see splines and smoothing) slots into the same algebra.

the identifiability constraint#

of all the linear constraints that would work, sum-to-zero gives the narrowest confidence bands: it pins the curve where the data are rather than at one end. the visible consequence is the “bow-tie” band of a term estimated as a straight line (see the mgcv example below).

penalised least squares#

stack the constrained bases into one design, \(\mathbf{X} = (\mathbf{1}, \mathbf{X}_1, \mathbf{X}_2)\) with \(\boldsymbol\beta = (\beta_0, \boldsymbol\delta^\top, \boldsymbol\gamma^\top)^\top\), and pad each penalty with zeros (\(\bar{\mathbf{S}}_1, \bar{\mathbf{S}}_2\)) so it acts on \(\boldsymbol\beta\) but touches only its own block.

this is ridge regression with a structured penalty matrix in place of \(\lambda \mathbf{I}\) (compare regularised regression). the intercept is never penalised, and each \(\lambda_j\) acts only on its own function: \(\lambda_1 \to \infty\) forces \(\mathbf{D}_1\boldsymbol\delta = \mathbf{0}\), a straight-line \(f_1\), while leaving \(f_2\) free. \(\mathbf{H}\) is no longer a projection, but \(\operatorname{tr}(\mathbf{H})\) is still the effective degrees of freedom, ranging from \(3\) (intercept plus two lines) up to \(1 + (k_1 - 1) + (k_2 - 1)\) (no penalty). the augmented form is how it is computed: lm() does the QR decomposition and \(\mathbf{X}^\top\mathbf{X}\) is never formed.

choosing the smoothing parameters#

the target is \(M = \frac1n \sum_i (\hat f_i - f_i)^2\), unknowable since \(f\) is. leave-one-out cross-validation, \(\mathcal{V}_{\mathrm{CV}} = \frac1n \sum_i (\hat f_i^{[-i]} - y_i)^2\), estimates \(\mathrm{E}[M] + \sigma^2\) up to the leave-one-out swap: the cross term \(\mathrm{E}[(\hat f_i^{[-i]} - f_i)\varepsilon_i]\) vanishes precisely because the left-out fit never saw \(\varepsilon_i\). (on the general idea see cross-validation.)

small \(\lambda\) shrinks the rss but inflates \(\operatorname{tr}(\mathbf{H})\); large \(\lambda\) the reverse; \(\mathcal{V}_g\) balances them, needs one trace rather than \(n\) leverages, and (unlike ocv) is invariant to rotations of the data. with two or more \(\lambda\)’s a grid is wasteful, so minimise \(\mathcal{V}_g\) numerically over \(\log \lambda_j\) — the log keeps every \(\lambda_j > 0\) and the search unconstrained.

in R, and the trees data#

am.fit <- function(y, x, v, sp, k = 10) {
  # tf.XD (lecture helper): tent basis at the data, last column dropped, columns centred; plus D
  xdx <- tf.XD(x, seq(min(x), max(x), length = k))
  xdv <- tf.XD(v, seq(min(v), max(v), length = k))
  nD <- nrow(xdx$D) * 2; sp <- sqrt(sp)
  X <- cbind(c(rep(1, nrow(xdx$X)), rep(0, nD)),           # unpenalised intercept
             rbind(xdx$X, sp[1] * xdx$D, xdv$D * 0),
             rbind(xdv$X, xdx$D * 0, sp[2] * xdv$D))       # augmented design
  b <- lm(c(y, rep(0, nD)) ~ X - 1)                       # OLS on augmented data
  n <- length(y); trA <- sum(influence(b)$hat[1:n])       # tr(H)
  rss <- sum((y - fitted(b)[1:n])^2)
  list(b = coef(b), edf = trA, fitted = fitted(b)[1:n], gcv = rss / (n - trA) * n / (n - trA))
}
am.gcv <- function(lsp, y, x, v, k) am.fit(y, x, v, exp(lsp), k)$gcv      # GCV as a function of log sp
fit.am <- optim(c(0, 0), am.gcv, y = trees$Volume, x = trees$Girth, v = trees$Height, k = 15)

proofs#

generalised additive models#

penalised irls#

a gam is fitted by penalised likelihood: maximise the log-likelihood minus \(\tfrac12\sum_j \lambda_j \boldsymbol\beta^\top\bar{\mathbf{S}}_j\boldsymbol\beta\). in practice that is the glm’s irls loop, unchanged except that each weighted least-squares step carries the penalty. from the current \(\hat{\boldsymbol\eta}\), \(\hat\mu_i = g^{-1}(\hat\eta_i)\):

1. weights and working response, exactly as for a glm: \(w_i = 1/[V(\hat\mu_i)\, g’(\hat\mu_i)^2]\), \(z_i = g’(\hat\mu_i)(y_i - \hat\mu_i) + \hat\eta_i\);
2. solve the penalised weighted least squares problem, \[ \min \|\sqrt{\mathbf{W}}(\mathbf{z} - \mathbf{X}\boldsymbol\beta)\|^2 + \sum_j \lambda_j \boldsymbol\beta^\top\bar{\mathbf{S}}_j\boldsymbol\beta\,, \] whose solution is \[ \hat{\boldsymbol\beta} = (\mathbf{X}^\top\mathbf{W}\mathbf{X} + \lambda_1\bar{\mathbf{S}}_1 + \lambda_2\bar{\mathbf{S}}_2)^{-1}\mathbf{X}^\top\mathbf{W}\mathbf{z}; \]
3. update \(\hat{\boldsymbol\eta} = \mathbf{X}\hat{\boldsymbol\beta}\) and repeat to convergence.

gam.fit <- function(y, x, v, sp, k = 10) {
  eta <- log(y); old.gcv <- -100
  repeat {
    mu <- exp(eta); z <- (y - mu) / mu + eta     # pseudo-data, w = 1
    fit <- am.fit(z, x, v, sp, k)                # one penalised LS step
    if (abs(fit$gcv - old.gcv) < 1e-5 * fit$gcv) break
    old.gcv <- fit$gcv; eta <- fit$fitted
  }
  fit$fitted <- exp(fit$fitted); fit
}

proofs#

gams in practice: mgcv#

mgcv::gam() has the glm() interface plus smooth terms s(x); by default each term’s smoothness is chosen by GCV (for the families used here), with a thin plate regression spline basis of dimension k = 10.

library(mgcv)
gam.mod <- gam(Volume ~ s(Height) + s(Girth), family = Gamma(link = log), data = trees)

Estimated degrees of freedom:
1.00 2.42  total = 4.42
GCV score: 0.008082356

height gets edf \(1.00\) (a centred line has one free parameter, its slope; the penalty has removed every wiggle), girth \(2.42\) (a little more than a quadratic), and the total adds one for the intercept. this matches the hand-built gam despite a different basis: the data choose the shape, not the basis. plot(gam.mod) draws each \(\hat f_j\) with \(\pm 2\) s.e. bands and partial residuals, \(\hat f_j(x_{ij}) + r_i^P\) — the data with every other term removed, on the working scale.

the defaults, one at a time#

• basis: thin plate is the default but computationally heavy; bs = "cr" gives an indistinguishable fit. the basis is a computational device.
• basis dimension k is the ceiling on a term’s df (\(k - 1\) after the constraint); the penalty decides how many are used. raising k once it is big enough changes nothing; too small a k cannot represent a wiggly effect at all.
• gamma multiplies the df in the GCV denominator, \(n\sum_i (r_i^P)^2 / [n - \gamma\operatorname{tr}(\mathbf{H})]^2\), charging each df \(\gamma\) times to counter GCV’s tendency to overfit. scores at different \(\gamma\) are not comparable (n.c.).
• several variables: s(x, v) is an isotropic smooth — rotation-invariant but not scale-invariant, so it struggles when covariates are in different units (inches against feet); te(x, v) is a tensor product with one penalty per direction — more parameters, but invariant to rescaling. its basis dimension is the product of marginal ones (\(5 \times 5 = 25\), comparable with k = 25).

reading the table: the plain additive fits have the lowest GCV, and the 2-d smooths buy nothing — the empirical form of the multiplicative argument (on the log scale girth and height already act additively), with the isotropic smooth worst, as its unit-sensitivity predicts. an edf of 1.00 for height also suggests the semi-parametric Volume ~ Height + s(Girth); gam() mixes parametric and smooth terms freely.

backfitting, and the gam package#

the tutorial uses Hastie’s gam package on ISLR’s Wage data. with natural-spline terms a gam is just one large least-squares regression on basis columns (lm(wage ~ ns(year, 4) + ns(age, 5) + education)). with smoothing-spline terms s(age, 5) it is fitted by backfitting: cycle through the terms, re-smoothing each partial residual \(y_i - \beta_0 - \sum_{k \ne j}\hat f_k(x_{ik})\) against \(x_{ij}\) until nothing changes. its summary() tests linear against non-linear per term (year \(p = 0.35\): linear suffices; age \(p < 2 \times 10^{-16}\)), and nested \(F\)-tests agree (no year to linear year \(F = 14.48\), \(p = 0.00014\); linear to spline \(F = 1.10\), \(p = 0.35\)).

see also#

• splines and smoothing — the one-predictor smoother each term is built from, and effective degrees of freedom
• generalised linear models — the irls loop that penalised irls extends by one penalty term
• regularised regression — ridge is penalised least squares with \(\mathbf{S} = \mathbf{I}\); the augmented-data trick works there too
• cross-validation — the general method behind ocv and GCV
• model selection — choosing between parametric and smooth terms
• curse of dimensionality — why additivity, not a full surface \(f(x_1, \dots, x_p)\), is the default
• locally weighted regression — the other route to a smooth term (lo() in the gam package)