How REML chooses smoothness#
A spline can follow the data as closely as you let it. The smoothing penalty decides how closely, and REML decides the penalty. This page says what that means, first in plain words, then in the maths, then in pictures on data where the true curve is known.
In plain words
A spline with forty pieces can copy the noise in the data as easily as the signal. The penalty charges the fit for wiggliness, and lambda is the price per unit of wiggle. REML sets the price by asking one question of the data: which amount of smoothness makes what we observed most probable, once the curve’s own uncertainty has been averaged out? No holdout set, no grid search; one criterion, maximised.
The maths#
The model is an additive predictor with one smooth per covariate:
where the \(b_{jk}\) are basis functions (here B-splines) and \(K_j\) is the
basis size, the k you pass to Spline. Fitting maximises a penalised
log-likelihood,
in which \(\ell\) is the ordinary log-likelihood of the family and each \(\mathbf S_j\) is a penalty matrix measuring the wiggliness of \(f_j\). For a P-spline the penalty is the sum of squared second differences of neighbouring coefficients,
a discrete stand-in for \(\int f''(x)^2\,\mathrm{d}x\). A straight line has zero second differences, so the penalty cannot charge for one: lines are the penalty’s null space.
How much of the basis the fit actually uses is the effective degrees of freedom,
where \(\mathbf X\) holds the basis functions evaluated at the data and
\(\mathbf W\) the working weights of the fit. With \(\lambda = 0\) the trace is
the rank of the identifiable design, one less than k here: the smooth is
centred, so its constant direction belongs to the intercept and is constrained
away, leaving 39 free degrees of freedom for a basis of 40. As
\(\lambda \to \infty\) it falls to the size of the null space, one line’s worth.
The EDF in every figure title below is this number.
REML chooses \(\boldsymbol\lambda\) by maximising the criterion Wood (2011) writes, for a fitted \(\hat{\boldsymbol\beta}\) at the given \(\boldsymbol\lambda\), as
with \(\mathbf H\) the negative Hessian of \(\ell\) at \(\hat{\boldsymbol\beta}\)
(for a GLM, \(\mathbf X^\top \mathbf W \mathbf X\)), \(|\cdot|_+\) the product
of the non-zero eigenvalues, and \(M_p\) the dimension of the null space. For
a Gaussian response this is exactly the restricted likelihood; for every
other family it is the Laplace approximation to it, which is what
fit_reml maximises.
Term |
What it does |
Why it matters |
|---|---|---|
\(\ell(\hat{\boldsymbol\beta})\) |
Rewards a fit that follows the data. |
On its own it would always choose \(\lambda = 0\). |
\(-\tfrac12 \hat{\boldsymbol\beta}^\top \mathbf S_{\boldsymbol\lambda} \hat{\boldsymbol\beta}\) |
Charges the fitted curve for its wiggle at the current price. |
The price is what is being chosen. |
\(+\tfrac12 \log\lvert\mathbf S_{\boldsymbol\lambda}\rvert_+\) |
Grows with \(\lambda\): the volume of curves the penalty considers plausible shrinks as the price rises. |
This is the term that rewards simplicity. |
\(-\tfrac12 \log\lvert\mathbf H + \mathbf S_{\boldsymbol\lambda}\rvert\) |
Falls with \(\lambda\): the volume of curves the data leave plausible. |
Together with the previous term it is the Occam factor: complexity is paid for automatically. |
The Bayesian reading makes the balance intuitive. The penalty is a prior \(\boldsymbol\beta \sim N(\mathbf 0, \mathbf S_{\boldsymbol\lambda}^{-})\) that prefers smooth curves, and \(\mathcal V\) is the log probability of the data with the curve integrated out. Maximising it is asking which smoothness makes the observed data most probable.
Why REML and not a holdout#
Cross-validation needs many refits and a holdout that is not always available. Generalised cross-validation is a single criterion but Reiss and Ogden (2009) showed it has more local optima than REML and tends to under-smooth. Wood (2011) gave a stable Newton method for \(\mathcal V\) that also yields the smoothing-parameter uncertainty, and Wood, Pya and Säfken (2016) extended it to any regular likelihood. The holdout curve in the figures below is the check an actuary still trusts; on this data REML lands where it bottoms out.
Removing a term altogether#
The penalty cannot shrink its own null space, so an ordinary smooth can
never disappear: at most it becomes a straight line. Marra and Wood (2011)
add a second penalty on that null space, \(\lambda_j^{*} \mathbf S_j^{*}\)
with \(\mathbf S_j^{*} = \mathbf U_j \mathbf U_j^\top\) built from the
null-space eigenvectors of \(\mathbf S_j\), and let REML estimate both
prices. That is select=True: a term with no signal can then be shrunk to
zero, which the last figure shows.
What this means for a tariff
A curve that follows noise is a price that follows noise, and a price that follows noise is one a competitor can pick off. REML gives a reproducible, defensible choice of smoothness that a reviewer can read off the summary: the penalty, the effective degrees of freedom, and the criterion value at the optimum.
Words used above#
Word |
Meaning here |
|---|---|
Basis size, |
How many pieces the spline is built from; the most flexible the curve can be. |
Penalty |
A number that grows with the wiggliness of the curve. |
Lambda |
The price per unit of penalty; large means smooth. |
EDF |
Effective degrees of freedom: how many of the |
Null space |
The shapes the penalty cannot charge for: straight lines. |
REML |
The criterion that chooses lambda from the data; for non-Gaussian families its Laplace approximation, sometimes written LAML. |
Holdout deviance |
The model’s error on rows it never saw; lower is better. |
See it happen#
Four hundred points on a sine wave with a gentle slope, plus Gaussian noise with standard deviation 0.45. The basis is a P-spline with 40 functions, deliberately generous, so that an unpenalised fit has room to misbehave.
Move the price yourself#
Drag the slider and the left panel refits at that lambda; the right panel says how many effective degrees of freedom that price leaves. Run REML walks the optimiser’s own iterates, numbered, ending on the yellow tick, which is the lambda REML settled on.
No penalty, REML, far too much penalty#
Three fits of the same model. The first fixes lambda at zero, the second lets
fit_reml choose it, the third fixes it at ten thousand.
No penalty, lambda = 0 EDF 39.0
REML EDF 8.1
Lambda = 10,000 EDF 2.4
The same 400 points three times. With no penalty the spline chases every point. REML picks a penalty that follows the truth. A penalty far too large leaves the curve a slope and one gentle bend, and it misses the peaks.#
With no penalty the fit spends all 39 of its free degrees of freedom on the noise. With the REML lambda it spends about eight, and the curve sits on the truth. At lambda ten thousand it has two left, enough for a slope and a little bend, and misses the peaks entirely.
What the penalty buys#
Now sweep lambda over a log grid of twelve fixed values and ask two things of each fit: how many effective degrees of freedom it keeps, and how well it predicts rows it did not see, measured by mean deviance over five folds.
REML chose lambda = 61.8
What the penalty buys. Left: how many effective parameters the spline keeps. Right: held-out deviance. The red line is the lambda REML chose without ever seeing a holdout.#
The left panel is the dial: each factor of ten in lambda takes away a few degrees of freedom. The right panel is the reason the dial matters. Held-out deviance is flat and high on the left, where every fit reproduces the noise, falls to a minimum, and rises steeply on the right, where the fits are too stiff to reach the peaks. REML never touched a fold and still landed at the bottom.
How REML gets there#
REML is an optimisation the solver drives directly, not a grid search.
superglm minimises the negative REML criterion, so lower is better and the
right-hand curve below falls. reml_diagnostics keeps the path the optimiser
took: one lambda per outer step plus the starting value, and one criterion
value per step.
REML is an optimisation the solver drives directly, not a grid search: a handful of steps from the starting value to the optimum on this example. Left, lambda after each step, where step 0 is the starting value, so the left panel carries one point more than the right. Right, the criterion superglm minimises, which is why the curve falls. In both panels the red marker is the value REML settled on.#
The path is not a steady climb. The optimiser probes downwards once, then climbs three orders of magnitude in two steps, overshoots, and settles back to a lambda near the bottom of the held-out curve in the previous figure. The criterion is within a fraction of a point of its final value after four steps; everything past that is refinement. A grid over twelve values, as in the previous figure, costs twelve fits and five folds each; the optimiser costs a handful of fits and no folds.
Removing a term that carries no signal#
Add a second column z that has nothing to do with y, and fit both columns
as splines. The ordinary penalty charges for bending, so a term it cannot
justify is shrunk to a straight line, and a straight line costs nothing under
that penalty, so it stays. select=True adds a second penalty on the straight
part as well, and REML can then take the term out altogether.
select=False x EDF 6.745
select=False z EDF 1.000
select=True x EDF 6.199
select=True z EDF 0.000
A term with no signal. The black curve is the fitted effect and the yellow band its 95% confidence interval. All four panels share one y axis, so the z term’s effect can be compared with the x term’s. Without the double penalty the z term keeps a slope and one degree of freedom, because the ordinary penalty cannot charge for a straight line; with select=True REML shrinks it to flat and its EDF to zero. The x term is untouched either way.#
References#
Wood, S. N. (2011). Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society: Series B, 73(1), 3–36. doi:10.1111/j.1467-9868.2010.00749.x
Wood, S. N., Pya, N., and Säfken, B. (2016). Smoothing parameter and model selection for general smooth models. Journal of the American Statistical Association, 111(516), 1548–1563. doi:10.1080/01621459.2016.1180986
Marra, G., and Wood, S. N. (2011). Practical variable selection for generalized additive models. Computational Statistics & Data Analysis, 55(7), 2372–2387. doi:10.1016/j.csda.2011.02.004
Reiss, P. T., and Ogden, R. T. (2009). Smoothing parameter selection for a class of semiparametric linear models. Journal of the Royal Statistical Society: Series B, 71(2), 505–523. doi:10.1111/j.1467-9868.2008.00695.x
Wood, S. N. (2017). Generalized Additive Models: An Introduction with R, 2nd edition. Chapman and Hall/CRC. doi:10.1201/9781315370279
Wahba, G. (1985). A comparison of GCV and GML for choosing the smoothing parameter in the generalized spline smoothing problem. Annals of Statistics, 13(4), 1378–1402. doi:10.1214/aos/1176349743
Main takeaways#
No penalty means the spline reproduces the noise; the effective degrees of freedom climb towards the basis size.
REML picks the penalty from the data alone, and on this example it lands where held-out deviance is lowest.
The optimiser reaches that value in a handful of steps; it is an optimisation the solver drives directly, not a grid search, and the criterion it drives down is the negative REML criterion, so lower is better.
select=Truelets REML remove a term that carries no signal instead of leaving it a straight line.
Next steps#
Solvers and internals for the REML criterion itself