Generalized Additive Models: splines, penalties and smoothing

How to let the data choose the shape of a relationship — without overfitting. Animations, equations, and Python you can run directly in this page.

📅 July 2026 ⏱️ ~15 min read 🏷️ Nonparametric Regression • Splines • GAM

← Back to blog Go to the penalty 🇫🇷 Lire en français

Contents

For recruiters
Why additive models?
From linear to additive
A smooth is a sum of bases
The roughness penalty and λ
Choosing λ automatically (GCV)
Additive structure
The “G”: link functions
In practice with mgcv
Limits and extensions
References

For recruiters (quick read)

A GAM replaces the coefficients of a linear regression with estimated smooth functions, while remaining interpretable. It is the middle ground between linear regression (too rigid) and black-box methods (unreadable). What this topic demonstrates:

Regularisation: the roughness penalty is ridge regression in a spline basis.
Bias-variance trade-off: controlled continuously by a single parameter λ, measurable through effective degrees of freedom.
Model selection: λ chosen by GCV or REML, without expensive cross-validation.
Numerical linear algebra: hat matrix, trace, well-conditioned penalised systems.

Penalised splines Regularisation GCV / REML mgcv Interpretable models

Why additive models?

Linear regression makes a strong assumption: the effect of every variable is a straight line. Most real phenomena are not that polite. Temperature raises crop yield up to an optimum, then destroys it. A drug dose helps, then harms. Traffic to a website follows a shape across the day, not a slope.

One could add \(x^2\) or \(x^3\) terms, or bin the variable into categories — but those choices are arbitrary and fragile. The generalized additive model asks the question differently: what if we let the data decide the shape of the curve, penalising only its irregularity?

01 / 06Opening titles: splines, penalisation, smoothing.

The idea in one line. A GAM is written \(\;g(\mu) = \beta_0 + \sum_j f_j(x_j)\;\): each \(f_j\) is a smooth function estimated from the data, with its flexibility controlled by a penalty rather than fixed in advance.

From linear to additive

Start from the linear model \(y = \beta_0 + \beta_1 x + \varepsilon\). Fitted by least squares to curved data, it leaves obvious structure in the residuals: they are not scattered at random about the line, they follow a pattern. That is the signature of systematic information the model failed to capture.

02 / 06 The line's residuals (orange) collapse once \(\beta_1 x\) is replaced by a smooth function \(f(x)\).

Linear model

Shape imposed

\[y_i = \beta_0 + \beta_1 x_i + \varepsilon_i\]

A single degree of freedom for the shape: the slope. High bias whenever the true relationship curves.

Additive model

Shape learned

\[y_i = \beta_0 + f(x_i) + \varepsilon_i\]

\(f\) lives in a space of smooth functions. It remains to define that space — and to stop \(f\) from passing through every point.

A smooth is a sum of bases

To estimate a function, write it as a linear combination of known elementary functions, called basis functions:

Basis representation

\[f(x) = \sum_{k=1}^{K} \beta_k\, B_k(x)\]

This step is decisive: estimating a function (an infinite-dimensional object) becomes estimating a vector of coefficients \(\beta\) — an ordinary linear regression problem on the design matrix \(B\). The usual \(B_k\) are B-splines: piecewise polynomial bumps, each non-zero only over a short interval. That locality is what makes the fit stable.

03 / 06 Basis functions appear, are scaled by their \(\beta_k\), then summed into a single curve.

Building the basis by hand

The code below really runs in your browser (Python compiled to WebAssembly). Edit the values and run it again: the first run takes a few seconds while Python loads.

PYTHON · NUMPY

# B-spline basis via the Cox-de Boor recursion (numpy only)
import numpy as np

def make_knots(xl, xr, n_inner, degree=3):
    inner = np.linspace(xl, xr, n_inner + 2)[1:-1]
    return np.r_[[xl]*(degree+1), inner, [xr]*(degree+1)]

def _cox(x, t, i, k):
    if k == 0:
        return ((x >= t[i]) & (x < t[i+1])).astype(float)
    c1 = np.zeros_like(x); c2 = np.zeros_like(x)
    d1 = t[i+k] - t[i]; d2 = t[i+k+1] - t[i+1]
    if d1 > 0: c1 = (x - t[i]) / d1 * _cox(x, t, i, k-1)
    if d2 > 0: c2 = (t[i+k+1] - x) / d2 * _cox(x, t, i+1, k-1)
    return c1 + c2

def bspline_basis(x, knots, degree=3):
    xl, xr = knots[0], knots[-1]
    x = np.clip(np.atleast_1d(np.asarray(x, float)), xl, xr - 1e-12)
    n = len(knots) - degree - 1
    return np.column_stack([_cox(x, knots, i, degree) for i in range(n)])

# --- try it: 8 interior knots, cubic ---
xe = np.linspace(0, 1, 300)
t  = make_knots(0, 1, n_inner=8, degree=3)
B  = bspline_basis(xe, t, 3)

print("design matrix:", B.shape, " -> K =", B.shape[1], "coefficients")
print("partition of unity (bases sum to 1):", np.allclose(B.sum(1), 1))

plot_basis(xe, B)          # draws the K basis functions

Why partition of unity matters. B-splines sum to 1 at every point. A constant is therefore exactly representable, and each coefficient \(\beta_k\) reads as the local height of the curve — which makes the penalty of the next section entirely natural.

The roughness penalty and λ

With \(K\) large, ordinary least squares threads the curve through nearly every point: overfitting. The remedy is not to shrink \(K\) — a coarse, discrete choice — but to add a penalty term that charges the fit for being wiggly:

Penalised criterion

\[\hat{\beta} = \arg\min_{\beta}\; \underbrace{\|y - B\beta\|^2}_{\text{fidelity to data}} \;+\; \lambda \underbrace{\beta^{\top} S \beta}_{\text{roughness}}\]

The matrix \(S\) measures irregularity. In the P-spline approach of Eilers and Marx one takes \(S = D^{\top}D\), where \(D\) is the second-difference operator: penalising \(\beta^\top S\beta\) amounts to penalising abrupt changes between neighbouring coefficients. The solution is explicit:

Penalised solution (ridge in the basis)

\[\hat{\beta} = (B^{\top}B + \lambda S)^{-1} B^{\top} y\]

04 / 06 λ swept from \(10^{-4}\) to \(10^{4}\): the curve stiffens and the effective degrees of freedom melt from ≈ 29 down to ≈ 2.

Effective degrees of freedom

How do we quantify the complexity actually used? Through the trace of the hat matrix \(H = B(B^{\top}B + \lambda S)^{-1}B^{\top}\), which projects \(y\) onto the fit:

Effective degrees of freedom

\[\mathrm{edf}(\lambda) = \operatorname{tr}(H)\]

This is the continuous generalisation of the “number of parameters”. As \(\lambda \to 0\), \(\mathrm{edf} \to K\); as \(\lambda \to \infty\), the second-difference penalty forces a straight line and \(\mathrm{edf} \to 2\). In between, complexity varies continuously — exactly what the animation shows as a number.

PYTHON · NUMPY

Choosing λ automatically (GCV)

Leaving λ to the user's judgement would be an admission of defeat. Leave-one-out cross-validation has a remarkable closed form here: for a linear smoother, the cross-validation error can be computed without ever refitting. Its stabilised version is generalized cross-validation:

GCV criterion

\[\mathrm{GCV}(\lambda) = \frac{n\,\|y - Hy\|^2}{\left(n - \operatorname{tr}(H)\right)^2}\]

The numerator rewards fit, the denominator punishes complexity: the larger \(\operatorname{tr}(H)\), the smaller the denominator and the worse the score. We minimise over a logarithmic grid of λ. In practice mgcv defaults to REML, which is more resistant to undersmoothing, but GCV remains the easiest to understand and to code.

PYTHON · NUMPY

On this data GCV selects \(\lambda \approx 0.93\), giving \(\mathrm{edf} \approx 11.5\) out of the 29 available coefficients. The model therefore uses only a third of its flexibility: the penalty performed model selection, continuously.

Additive structure

Everything above concerned a single variable. The model's strength comes from the assembly: several smooths are summed, one per covariate, each with its own penalty \(\lambda_j\).

Additive model

\[y_i = \beta_0 + f_1(x_{1i}) + f_2(x_{2i}) + \dots + f_p(x_{pi}) + \varepsilon_i\]

05 / 06 Three partial effects, each readable on its own — the decisive advantage of GAMs over black boxes.

This is where interpretability is won. Gradient boosting or a neural network would capture these relationships too, but without letting you plot “the effect of \(x_2\), all else held equal”. Here each \(f_j\) is a standalone graphical object, with its confidence band. The price: interactions are not captured by default — they must be requested explicitly (see Limits).

The “G”: link functions

So far the response was continuous and the noise Gaussian. The “generalized” lifts that restriction, exactly as GLMs extend linear regression: the response may come from the exponential family (Poisson for counts, binomial for proportions, Gamma for durations), connected to the predictor by a link function \(g\).

Generalized additive model

\[g(\mu_i) = \beta_0 + \sum_{j} f_j(x_{ji}), \qquad y_i \sim \text{exponential family},\;\; \mathbb{E}[y_i] = \mu_i\]

06 / 06 Identity, log, logit: the link adapts the model to the nature of the response.

Nature of the response	Family	Canonical link
Continuous, symmetric	Gaussian	\(g(\mu) = \mu\)
Counts (0, 1, 2, …)	Poisson	\(g(\mu) = \log \mu\)
Binary or proportion	Binomial	\(g(\mu) = \log\frac{\mu}{1-\mu}\)
Positive, skewed	Gamma	\(g(\mu) = 1/\mu\) (ou log)
Overdispersed counts	Negative binomial	\(g(\mu) = \log \mu\)

Fitting then proceeds by penalised iteratively reweighted least squares (P-IRLS): the linear scheme above is iterated on a working response, with the weights updated at each pass.

In practice with mgcv

Everything we coded by hand is industrialised in the R package mgcv, written by Simon Wood himself. The code below does not run in the browser (it requires R), but it is the natural destination of this article:

# --- R / mgcv --- library(mgcv) # one smooth per covariate; k = basis dimension, bs = spline type m <- gam(y ~ s(x1, k = 20, bs = "tp") + s(x2) + x3, family = gaussian, method = "REML", data = df) summary(m) # edf per term, approximate p-values, deviance explained plot(m, pages = 1, shade = TRUE) # the partial effects, one panel each gam.check(m) # diagnostics: is k large enough? residuals sane? # count response, with a 2-D interaction and a cyclic effect m2 <- gam(n ~ te(temp, day) + s(month, bs = "cc", k = 12), family = nb, method = "REML", data = df)

Argument	Role	Guidance
`k`	Basis dimension (our \(K\))	Be generous; the penalty sorts it out. Verify with `gam.check`.
`bs = "tp"`	Thin-plate spline	Robust default, rotation-invariant, works in more than one dimension.
`bs = "cc"`	Cyclic spline	For month, hour, angle — enforces \(f(\text{start}) = f(\text{end})\).
`method = "REML"`	Smoothness selection	Preferred over GCV: less prone to undersmoothing.
`select = TRUE`	Extra shrinkage	Allows a term to be penalised to exactly zero.

Limits and extensions

Assumption or limit	If it bites → remedy
Purely additive effects	Tensor products: `te(x1, x2)` for a smooth interaction across different scales
Independent observations	GAMM: random effects, temporal or spatial correlation
Only the mean is modelled	GAMLSS: smooths on variance, skewness and kurtosis too
Large datasets	`bam()`: block-wise fitting, parallelisable
p-values for smooth terms	Approximate (the edf is itself estimated) — read them with care
Extrapolation beyond the data	Dangerous: the penalty constrains nothing outside the observed range

The classic trap. Choosing \(k\) too small breaks everything: no penalty can rescue a basis incapable of representing the true function. Choosing \(k\) too large costs computation time and nothing else. When in doubt, raise \(k\) and let λ decide.

References

Wood, S. N. (2017). Generalized Additive Models: An Introduction with R, 2nd ed., Chapman & Hall/CRC — the reference, and the source of this article.
Eilers, P. H. C. & Marx, B. D. (1996). Flexible smoothing with B-splines and penalties, Statistical Science — the founding P-spline paper.
Hastie, T. & Tibshirani, R. (1990). Generalized Additive Models, Chapman & Hall — where the model originates.
Online resources: mgcv on CRAN • GAMs in R (Noam Ross) • Pyodide (the Python engine of this page)

The animations in this article were produced with Manim; the source for all six scenes is available on GitHub.