Splines and Flexible Functional Forms

In plain language

Splines and flexible functional forms are ways to model a curved (nonlinear) relationship between a continuous measurement — like BMI or age — and an outcome in a regression, without forcing the relationship to follow an arbitrary straight line or invented cut-off. The most common tool is the restricted cubic spline, which fits a smooth, piecewise curve through the data with a small number of anchor points called knots placed at percentiles of the exposure. The result is reported as a dose-response picture showing predicted risk across the full range of the exposure, not as raw regression numbers — because individual spline coefficients have no meaningful interpretation on their own. These methods are the correct alternative both to grouping patients into BMI categories (which throws away information and invents thresholds) and to assuming a straight-line relationship that may be J-shaped, U-shaped, or otherwise curved.

The problem with categorizing continuous exposures

One of the most persistent analytic errors in epidemiology and HEOR is converting a continuous covariate into categories: age bands, BMI ≥ 30 for "obese," creatinine above a clinical threshold. The appeal is transparency — indicator variables are simple to include and easy to present in tables. The cost is severe. First, within-category variation is entirely discarded: a BMI of 30.1 and a BMI of 39.9 are treated as identical, despite very different health risk profiles. Second, the cut-point is invented, not derived from biology: BMI ≥ 30 is a clinical and regulatory convention; placing the cut at 28 or 32 would produce a different result with no principled basis for the choice. Third, residual confounding increases because a coarser representation of a confounder leaves more of its effect uncontrolled. The correct solution — always — is to model continuous covariates continuously, and when the shape of the relationship is uncertain, let the data determine that shape rather than imposing an arbitrary one.

When the linear assumption fails

Entering a continuous variable as a single linear term is an improvement over categorization, but it imposes a constant slope: every one-unit increase in the covariate is assumed to carry the same increment in the outcome on the model scale. This is implausible for many exposures. Body-mass index has a well-documented J-shaped or U-shaped association with all-cause mortality: both underweight and severely obese individuals face elevated risk while the normal-weight range carries the lowest risk. A linear term misrepresents this as a simple upward (or downward) trend. Similar non-linearity appears in the alcohol-health relationship, in environmental exposures with threshold or saturation effects, and in pharmacological dose-response curves. Analysts can diagnose linearity failure by plotting component-plus-residual (partial residual) plots against the covariate, overlaying a lowess smoother, and comparing AIC/BIC for linear versus nonlinear specifications. A clear departure from a straight partial-residual line is grounds to fit a flexible functional form.

Restricted cubic splines: mechanics for practitioners

A spline is a piecewise polynomial that joins smoothly at breakpoints called knots. In a cubic spline the function value, first derivative, and second derivative all agree at each knot, producing a smooth continuous curve. Restricted cubic splines (RCS) — also called natural cubic splines — add a further constraint: outside the two outermost (boundary) knots, the function is forced to be linear. This boundary-linearity constraint prevents the explosive extrapolation that unrestricted cubic splines can produce in the tails, and it consumes fewer degrees of freedom: a k-knot RCS uses k − 2 degrees of freedom. With 3 knots the spline adds only 1 df beyond a linear term; with 5 knots it adds 3 df.

The conventional knot-count guidance (Harrell, Regression Modeling Strategies, 2015): - 3 knots — small samples or strong prior knowledge of a simple curve; 1 df added. - 4 knots — the practical default for most HEOR analyses; 2 df added. - 5 knots — large datasets (n > 1,000 events) with complex shapes; 3 df added.

Knot placement follows default percentile rules: for 3 knots, the 10th, 50th, and 90th percentiles of the covariate distribution; for 5 knots, the 5th, 27.5th, 50th, 72.5th, and 95th percentiles. Simulation studies consistently show that results are much more sensitive to the number of knots than to their exact positions within a few percentile points, so the default rules are robust in practice. Analysts need not search for optimal knot locations.

The dose-response curve as the primary deliverable

This is a critical communication point: individual coefficients from a spline model are not interpretable in isolation. A 3-knot RCS produces two regression parameters (one linear term, one nonlinear component); neither has a stand-alone clinical meaning. The correct output is the predicted outcome value plotted as a function of the continuous covariate across its observed range, with a reference value chosen (e.g., BMI = 25 kg/m²) and the outcome expressed relative to that reference. For a logistic regression this is the log-odds ratio (or exponentiated OR) with a 95% confidence interval band. For a Cox model it is the log-HR versus the reference. This dose-response curve is the primary deliverable: it shows the shape, direction, and uncertainty of the exposure-outcome relationship in a single figure that decision-makers can interpret without statistical expertise. Reports that present only the numerical spline coefficients are providing uninterpretable output.

Fractional polynomials as a rival

Fractional polynomials (FP) model nonlinearity by selecting the best-fitting powers of X from the closed set {−2, −1, −0.5, 0, 0.5, 1, 2, 3} (where power 0 means log(X)), with a model-selection algorithm (MFP) that chooses one or two power terms. FPs avoid knots entirely and produce a smooth, closed-form transformation. They handle extreme-tail behaviour well when the correct transformation happens to be one of the available powers (e.g., log for right-skewed exposures); RCS is more flexible for complex interior shapes (multiple inflection points or humps) and for exposure ranges spanning several orders of magnitude. In practice both methods are occasionally compared as a sensitivity check; systematic comparisons in simulation show RCS and FPs perform similarly across a broad range of true shapes with neither dominating uniformly.

Royston-Parmar flexible parametric survival models

In survival analysis, Royston-Parmar (RP) models place a restricted cubic spline on the log cumulative hazard (for a proportional-hazards parameterisation) or log cumulative odds (for proportional odds). This produces a fully parametric model that can capture complex, non-monotone hazard shapes — early-rise-then-plateau, hump-shaped — that the six standard parametric distributions (exponential, Weibull, Gompertz, log-normal, log-logistic, generalised gamma) cannot accommodate. In an etiological study the primary output is the dose-response curve of the log hazard ratio as a function of a continuous covariate modelled with an RCS term, not the baseline hazard spline coefficients. For the HTA survival extrapolation context — where RP splines are used to extrapolate beyond observed follow-up and their knot positions are a key sensitivity parameter — see the dedicated entry on survival-extrapolation-hta-rwe. This entry focuses on their use as flexible functional forms within a regression framework, not on the extrapolation problem specifically.

Penalized splines and generalised additive models

Penalized splines (P-splines) and thin-plate splines in a generalised additive model (GAM) automate smoothness selection: instead of pre-specifying the number of knots, the analyst specifies a penalty on the second derivative of the smooth function, and the effective degrees of freedom are chosen by restricted maximum likelihood (REML). The mgcv package in R implements this via `gam(y ~ s(x), family = ...)`. GAMs are well-suited for exploratory analysis of many continuous predictors simultaneously or when knot selection is genuinely uncertain. Trade-offs: the estimated effective df is data-adaptive and the inference is approximate (confidence intervals rely on Bayesian posterior arguments); the smooth cannot easily be transported to other software as a fixed analytical formula; and the automated smoothness parameter can over-smooth in sparse tails. For a confirmatory analysis with a specific continuous exposure, a pre-specified RCS with a defended number of knots is more transparent and easier to pre-register and validate.

Overfitting guardrails

Every additional knot costs one degree of freedom. The standard rule of thumb: require at least 10–15 outcome events (EPV) per degree of freedom added by the spline. For a 5-knot RCS on one covariate (3 df) in a binary outcome study, this implies at least 30–45 events for that spline term before accounting for other covariates in the model. When EPV is low, reduce the knot count to 3 or to a linear term. Internal validation using bootstrap optimism correction (Harrell's `validate()` in R's rms package) quantifies how much the apparent C-statistic or calibration slope overfits the training data. AIC and BIC across models with different knot counts provide a complementary selection criterion.

Pros, cons, and trade-offs

Restricted cubic splines: Pros — flexible, can represent any smooth shape; boundary linearity prevents tail explosion; default percentile knot rules are robust; fully integrated with standard regression models and their inference machinery (Wald tests, bootstrap CIs); widely implemented in R (rms, splines), Python (patsy, formulaic), and SAS (EFFECT statement). Cons — individual coefficients are uninterpretable; requires adequate EPV; linear outside boundary knots may be wrong for structural extrapolation; does not automatically select the optimal knot count. Prefer when the functional form is uncertain and EPV ≥ 10 per spline df.

Fractional polynomials: Pros — closed-form transformation; no knot selection; can represent extreme-tail behaviour well with standard powers. Cons — restricted to the available power library; may fail for complex interior shapes; model search inflates effective degrees of freedom. Prefer when prior theory suggests a specific power transformation and the exposure range is wide.

Linear term: Pros — interpretable slope and CI; minimum df; analytically and computationally trivial; always valid as a starting point. Cons — systematically wrong when the true shape is curved; misses J-shaped or U-shaped patterns; can produce residual confounding bias when the covariate is an important confounder with nonlinear effects. Prefer when EPV is too small for flexible forms, or when prior evidence strongly supports linearity.

Categorization: Do not prefer categorization over any of the above for continuous exposures or confounders in a regression model. The only legitimate use of categories after fitting is for clinical communication — presenting the fitted curve evaluated at clinically recognisable BMI bands, for example — not as a substitute for modelling the continuous relationship.

When to use

Use splines and flexible functional forms whenever: (1) a continuous covariate is entered as an exposure or as a confounder in a regression model and the functional form is not firmly established by prior literature; (2) the exposure range is wide enough that a constant-slope assumption is non-trivial; (3) EPV ≥ 10 per spline df. Key applications in HEOR include: age when used as a confounder (rather than a study-design stratum); BMI, blood pressure, creatinine, HbA1c in outcome models; cumulative drug dose in dose-response analyses; biomarker-outcome relationships in registry or EHR studies. In propensity score models, RCS for continuous confounders is particularly important: a linear propensity score term for a nonlinearly confounded variable leaves systematic residual confounding after matching or weighting.

When NOT to use

Do not use splines when: (1) EPV is very low (< 10 per df) — reduce to a linear term and validate; (2) the model output must be a portable, closed-form risk score for bedside use — RCS requires software to evaluate and cannot be simplified to a lookup table; (3) the continuous covariate takes only a handful of distinct integer values (0, 1, 2, 3 prior hospitalisations) — use indicator variables instead; (4) extrapolation beyond the observed range is the scientific goal — the boundary linearity of RCS reduces but does not eliminate extrapolation risk, and structural extrapolation should use explicitly motivated parametric forms (Weibull, Royston-Parmar with external anchoring); (5) the association is already established as linear in multiple prior studies — the added complexity buys nothing and costs interpretability.

Interpreting the output

In the worked example, logistic regression of 5-year cardiovascular mortality on BMI is fit with a linear term (AIC = 2450) and a 3-knot restricted cubic spline (AIC = 2410). The spline model produces two regression coefficients; neither is reported individually. Instead, the dose-response curve is plotted: log-odds ratio of 5-year mortality versus BMI, evaluated across the observed BMI range with BMI = 28 kg/m² (the median knot and reference value) set to log-OR = 0.

(1) Formal interpretation. The individual parameters β₁ and β₂ of the 3-knot RCS basis are not estimands with clinical meaning; they are coordinates of a basis representation. The estimand is the conditional log-odds ratio of 5-year CV mortality at each BMI level relative to the reference BMI = 28, adjusting for all other model covariates, evaluated from the fitted spline function. The AIC difference of 40 units (AIC linear 2450 minus AIC spline 2410 = 40) provides strong evidence against the linear constraint: by convention a ΔAIC > 10 is considered strong evidence that the additional parameters improve fit beyond chance. The spline model uses 1 extra degree of freedom (k − 2 = 3 − 2 = 1) relative to the linear model; the improvement of 40 AIC units for 1 df is substantial.

(2) Practical interpretation. The dose-response curve reveals a J-shaped BMI-mortality relationship: predicted risk is lowest near the median (BMI ≈ 28) and rises at both the underweight extreme (BMI < 20) and the severely obese extreme (BMI > 38). A linear model forced this into a single slope, systematically under-estimating risk at both extremes. For a payer or HTA body, the shape matters: a policy targeting exclusively "BMI ≥ 30" misses elevated risk among the underweight, while an intervention designed from the correctly specified curve would target both tails.

Worked example

Scenario

A cohort study examines the association between body mass index (BMI) and 5-year cardiovascular mortality in 800 adults. Three specifications of BMI are compared in a logistic regression: (1) a binary obesity indicator (BMI ≥ 30 kg/m²), (2) BMI as a continuous linear term, and (3) BMI modeled with a 3-knot restricted cubic spline. Knot positions for the spline are determined from a 10-person pilot subsample. The analyst uses AIC to compare model fit and then plots the dose-response curve from the winning spline model to visualize the shape of the BMI-mortality relationship.

Dataset

Top panel: the 10-person pilot subsample sorted by BMI, used to place knots at the 10th, 50th, and 90th percentiles. Bottom panel: full-cohort AIC comparison across the three model specifications.

Steps

• Sort the 10-person pilot BMI values in ascending order: 18, 20, 23, 25, 27, 29, 31, 34, 38, 42 (n = 10 values).

• Knot 1 at the 10th percentile: in a sorted list of 10, this is the 1st value = 18 kg/m².

• Knot 2 at the 50th percentile (median): average of the 5th and 6th sorted values = (27 + 29) / 2 = 28 kg/m².

• Knot 3 at the 90th percentile: the 9th value = 38 kg/m². Three knots placed at 18, 28, 38.

• The binary model collapses BMI into obese vs not; every BMI difference within a category is discarded. AIC = 2480.

• The linear model uses one slope for the entire BMI range. AIC = 2450. AIC improvement over binary: 2480 - 2450 = 30.

• The 3-knot spline adds 1 extra degree of freedom (k - 2 = 3 - 2 = 1) on the full cohort with knots at 18, 28, 38. AIC = 2410. AIC improvement over linear: 2450 - 2410 = 40.

• A difference of 40 AIC units strongly favors the spline (rule of thumb: ΔAIC > 10 is strong evidence). The dose-response curve is plotted with BMI = 28 as the reference value. Individual spline coefficients are not reported because they have no stand-alone interpretation.

Result

Knot positions: 18, 28, 38 kg/m². Knot 2 = (27 + 29) / 2 = 28. AIC linear = 2450, AIC spline = 2410; AIC difference = 2450 - 2410 = 40 in favor of the spline. The dose-response curve reveals a J-shaped BMI-mortality association invisible in the linear or binary models.

Runnable example

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt

# ── Simulate J-shaped BMI-mortality data (nadir near BMI 26) ──
rng = np.random.default_rng(42)
n = 800
bmi = rng.uniform(15, 50, n)
true_log_odds = 0.004 * (bmi - 26) ** 2 - 1.8   # U-shaped truth
prob_true = 1 / (1 + np.exp(-true_log_odds))
event = rng.binomial(1, prob_true, n)
df = pd.DataFrame({"bmi": bmi, "event": event})

# ── Knot positions at 10th / 50th / 90th percentiles of observed BMI ──
p10, p50, p90 = (float(np.percentile(bmi, p)) for p in [10, 50, 90])
print(f"Knots at 10th/50th/90th pctiles: {p10:.1f} / {p50:.1f} / {p90:.1f} kg/m²")

# ── Fit linear and 3-knot RCS models ──
fit_lin = smf.logit("event ~ bmi", data=df).fit(disp=0)

# patsy cr(): natural cubic spline = RCS with boundary linearity.
# One interior knot (at p50) with boundary knots at p10, p90 gives a 3-knot spline
# and adds exactly 1 df beyond the linear term (k - 2 = 3 - 2 = 1).
rcs_formula = (
    f"event ~ cr(bmi, knots=[{p50:.2f}], "
    f"lower_bound={p10:.2f}, upper_bound={p90:.2f})"
)
fit_rcs = smf.logit(rcs_formula, data=df).fit(disp=0)

print(f"\nAIC  linear : {fit_lin.aic:.1f}")
print(f"AIC  spline : {fit_rcs.aic:.1f}")
print(f"ΔAIC (linear - spline): {fit_lin.aic - fit_rcs.aic:.1f}")
print("\nWARNING: spline coefficients below have NO stand-alone clinical interpretation.")
print("Report the dose-response curve (plot), not these coefficients.")
print(fit_rcs.summary2().tables[1][["Coef.", "Std.Err.", "P>|z|"]])

# ── Dose-response curve: log-OR vs reference BMI (= median) ──
ref_bmi = p50
bmi_grid = pd.DataFrame({"bmi": np.linspace(float(bmi.min()), float(bmi.max()), 200)})
ref_df   = pd.DataFrame({"bmi": [ref_bmi]})

def prob_to_logodds(p):
    p = np.clip(p, 1e-9, 1 - 1e-9)
    return np.log(p / (1 - p))

log_or = (prob_to_logodds(fit_rcs.predict(bmi_grid))
          - float(prob_to_logodds(fit_rcs.predict(ref_df))))

plt.figure(figsize=(7, 4))
plt.axhline(0, color="grey", lw=0.8, ls="--")
plt.axvline(ref_bmi, color="grey", lw=0.8, ls=":",
            label=f"Reference BMI = {ref_bmi:.1f}")
plt.plot(bmi_grid["bmi"], log_or, color="steelblue", lw=2,
         label="3-knot RCS log-OR")
plt.xlabel("BMI (kg/m²)")
plt.ylabel(f"Log-OR vs BMI = {ref_bmi:.1f} kg/m²")
plt.title("Dose-response: BMI and 5-year CV mortality")
plt.legend()
plt.tight_layout()
plt.savefig("bmi_spline_curve.png", dpi=120)
print(f"\nDose-response curve saved to bmi_spline_curve.png")
print("This curve is the primary deliverable. The spline coefficients are not reported.")

library(rms)

# ── Simulate the same J-shaped BMI-mortality data ──
set.seed(42)
n    <- 800
bmi  <- runif(n, 15, 50)
prob <- plogis(0.004 * (bmi - 26)^2 - 1.8)
df   <- data.frame(bmi = bmi, event = rbinom(n, 1, prob))

# ── rms requires datadist for Predict() and confidence intervals ──
dd <- datadist(df)
options(datadist = "dd")

# ── Fit linear vs 3-knot RCS logistic regression ──
fit_lin <- lrm(event ~ bmi,        data = df, x = TRUE, y = TRUE)
fit_rcs <- lrm(event ~ rcs(bmi, 3), data = df, x = TRUE, y = TRUE)
# rcs(bmi, 3): 3-knot RCS with knots at 10th/50th/90th percentiles by default.
# Knot positions are stored in attr(rcs(bmi, 3), "parms").

cat(sprintf("AIC  linear : %.1f\n", AIC(fit_lin)))
cat(sprintf("AIC  spline : %.1f\n", AIC(fit_rcs)))
cat(sprintf("ΔAIC (linear - spline): %.1f\n", AIC(fit_lin) - AIC(fit_rcs)))

# Wald test that the nonlinear spline term = 0 (tests whether spline improves on linear)
anova(fit_rcs)   # shows separate tests for linear vs nonlinear components

# ── Dose-response curve: log-OR vs reference BMI (= median, set by datadist) ──
# Predict() evaluates the fitted model on a grid and subtracts the reference value.
# The reference is the median of bmi in the datadist object (adjustable via ref.zero).
p <- Predict(fit_rcs, bmi = seq(15, 50, by = 0.5), fun = identity, ref.zero = TRUE)
# ref.zero = TRUE: log-OR is set to 0 at the reference value (datadist median).

plot(p,
     xlab = "BMI (kg/m²)",
     ylab = "Log-OR vs reference BMI (median)",
     main = "Dose-response: BMI and 5-year CV mortality (3-knot RCS)")
abline(h = 0, lty = 2, col = "grey")

# ── Alternative: splines::ns() for use outside rms (e.g., glm, coxph) ──
# ns() = natural cubic spline = RCS with boundary linearity.
# Knots must be specified explicitly (not automatic percentiles like rms::rcs).
k   <- quantile(df$bmi, c(0.10, 0.50, 0.90))
fit_ns <- glm(event ~ splines::ns(bmi, knots = k["50%"],
                                  Boundary.knots = c(k["10%"], k["90%"])),
              data = df, family = binomial)
cat(sprintf("AIC via ns(): %.1f\n", AIC(fit_ns)))
# WARNING: ns() coefficients are not interpretable individually.
# Use predict() on a grid and plot log-OR vs reference, as above.