Splines, Basis Expansions, and Generalized Additive Models

A generalized additive model replaces a single linear coefficient with a smooth function for selected features, then adds those functions together through a link. That captures nonlinear main effects without surrendering the model to a black box.

The additive assumption is powerful and limiting. Interactions must be added deliberately rather than appearing automatically.

The class has a founding paper, and its definition is precise enough to argue with. Hastie and Tibshirani wrote it in Statistical Science in 1986. They introduce “the class of generalized additive models which replaces the linear form ∑ βjXj by a sum of smooth functions ∑ sj(Xj)”, where “the sj(·)’s are unspecified functions that are estimated using a scatterplot smoother, in an iterative procedure we call the local scoring algorithm”. Two words in that sentence carry the whole trade. The functions are unspecified, which is where the flexibility comes from. And they are summed, which is where the interactions do not.

Spline and basis expansions turn nonlinear shape into a linear combination of constructed columns.

The curve is never drawn freehand. It is a set of columns and a set of weights. That is why the same regression or likelihood machinery that fits a straight line also fits a bend, and why a penalty on the weights is enough to control how much bending is allowed.

Each representation allocates flexibility differently, and the choice is not academic. Two public-health agencies faced the same task — draw a reference curve of a child’s body measurement against age — and shipped different answers.

The CDC’s 2000 US growth charts took the polynomial route. The report’s table “Summary of curve smoothing procedures” records, for weight-for-age, that “All results from the 3-parameter linear model and LWR were combined and refit at midpoints of age intervals from birth to 20 years using a 10-parameter polynomial regression model for boys and 9-parameter polynomial regression model for girls.” Different charts in the same report were smoothed by different procedures again: a 10-parameter nonlinear model for stature-for-age, a 5-parameter polynomial for weight-for-length, a 4-parameter polynomial for BMI-for-age. One published instrument, four fitting recipes, chosen chart by chart.

The WHO’s 2007 growth reference for 5-19 year olds took the other route. It used the Box-Cox power exponential distribution with curve smoothing by cubic splines. The 1977 NCHS reference sample of 22,917 children was merged with the under-fives cross-sectional sample, giving a final fitting sample of 30,907 observations.

That is the trade-off in shipped parameter counts. Nine or ten globally coupled coefficients stretched across birth to 20 years, where moving one changes the curve everywhere. Against a local basis whose flexibility is set by knot placement and a penalty, and whose adjustments stay near where they are made.

Smooth terms can reveal useful structure. They remain easier to inspect than an unrestricted surface.

One such shape has been measured at a scale that makes the cost of a straight line explicit. In 2015 The Lancet published an analysis of 74,225,200 deaths in 384 locations across 13 countries. Gasparrini and colleagues attributed 7.71% of that mortality to non-optimum temperature — 7.29% to cold against 0.42% to heat. Zhao and colleagues reproduced the asymmetry globally in The Lancet Planetary Health in 2021: 5,083,173 deaths a year, 9.43% of all deaths, 8.52% cold and 0.91% heat. A single slope through temperature would collapse a curve whose two arms differ by a factor of about seventeen. It would report the average of a rise and a fall as something near nothing.

The basis carrying that curve is not left implicit either. The Methods section of the 2021 study spells it out: “cb is the cross-basis function of two natural cubic splines for temperature during 0–21 lag days, with three internal knots placed at the tenth, 75th, and 90th temperature percentiles”. Knots placed at percentiles rather than at round temperatures is flexibility spent where the observations are. The lag window is the admission that today’s cold is tomorrow’s death.

Too little penalization lets the curve chase sample noise. Too much returns a nearly straight effect. Knots, basis dimension, penalty strength and boundary constraints all shape the result.

Which of those knobs actually decides the picture has been surveyed rather than guessed. A 2019 review of spline function procedures in R screened 519 packages and kept 229 as relevant. Its conclusion: “Most differences can be attributed to the choice of hyper-parameters rather than the basis used.” The long argument about which basis family to use is mostly not the argument that determines the fitted curve.

The selection criterion, by contrast, can fail outright. Simon Wood reported in 2011 that REML and ML smoothness selection avoids the occasional severe undersmoothing failures to which GCV and AIC are prone. So “too little penalization” is not only a dial set carelessly by a human. It is a documented failure mode of two criteria in common use, and it has a known alternative.

Use cross-validation and stability plots, then inspect whether small data changes create new wiggles.

Start with domain expectations and earn each degree of flexibility. That sequence is not merely a matter of taste. Two drug regulators have formally qualified a workflow built on it. The EMA’s CHMP adopted its qualification opinion of MCP-Mod on 23 January 2014, and FDA’s CDER issued a fit-for-purpose determination letter for the same methodology on 26 May 2016. Both concluded it uses the available data better than traditional pairwise comparisons.

Under that methodology the plausible shapes of the dose-response curve are written down as a candidate set of models before the trial runs, tested against the data, and only then fitted. The CHMP gave the reason: “Considering dose in its proper functional form, i.e. as continuous rather than a qualitative, ordered categorical variable also offers advantages in terms of maximising the use of the available information through modelling and by allowing the interpolation of information across the dose range.”

Steps two and three below — select candidate smooths, choose basis and constraints — are exactly that pre-specification. The regulatory version differs from the everyday one in a single respect worth copying. The shape hypotheses are committed to before anyone has seen the outcome.

Suppose temperature matters only under high pressure. Two separate smooth terms cannot represent the conditional effect. A tensor-product smooth or an engineered interaction can model the surface, but interpretation and data requirements become more demanding.

Use residual slices and domain mechanisms to decide which interactions deserve complexity.

At boundaries and beyond observed ranges, basis functions can continue according to mathematical constraints rather than empirical evidence. Some systems use linear tails, clipping, monotonic constraints, or fallback rules. Show the data density beneath every effect plot, and define behavior for unsupported values.

For a global polynomial the damage can be read straight off the weights. Gelman and Imbens made the case in 2017, in the Journal of Business & Economic Statistics: high-order polynomials should not be used in regression discontinuity designs. A polynomial regression estimate is a weighted average of the outcomes. The values of those weights, they observe, “have nothing to do with the actual shape” of the underlying function. They follow from the estimator alone, and can be inspected before the outcome data is seen.

Their example is Matsudaira’s summer-school data: 68,798 students, with a forcing variable ranging over [−199, 168]. They tabulate the normalised weight given to the individual at the extreme value 168, against an average weight of 1. It is −9.1 at first order, then 34.4, −32.8, 16.7, −7.4 and 2.6 as the degree rises to six. They tie it to Runge’s phenomenon. Given a set of N pairs (xi, yi) on a compact interval [a, b], the N − 1th-order polynomial through all of them “becomes increasingly erratic, as the number of points increases, close to the boundary of the interval”.

A curve drawn confidently past the last dense region is an average dominated by the rows nobody has.