AdvancedIndependently reviewed, not yet spot-checked by a human

Restricted cubic splines and the dose-response curve

How to drop the linearity assumption: where the knots go, why a zero-inflated clinical variable collides with the default knot placement, how to choose the reference point, why the confidence band pinches to a point there, and why the p value for the overall association and the p value for non-linearity must be reported separately.

Filling in something this site keeps recommending

Three pages already tell you to “use a restricted cubic spline instead”: linear regression where it discusses non-linearity, regression diagnostics where it reads residual plots, and the Cox model where it checks the assumption on continuous covariates. This page is where those three pointers land.

Start with the assumption being dropped. Put a continuous variable straight into a model:

\log h(t) = \log h_0(t) + \beta \cdot x + \cdots

The term $\beta \cdot x$ makes a strong claim: every one-unit increase in $x$ matters exactly as much as every other one. Going from 0 to 1 positive lymph node and going from 15 to 16 are forced to contribute the same number to the log hazard. Almost no clinician believes that, and yet it is the default, and nothing warns you about it.

A restricted cubic spline (RCS — in R the usual implementation is splines::ns(), a natural spline) cuts the range of $x$ at a few knots, fits a cubic piece between them, requires the pieces to join smoothly, and constrains the two tails to be straight. That constraint is what “restricted” and “natural” mean: without it, the sparsest parts of the data — the two ends — are exactly where the curve flails hardest.

The example on this page

survival::rotterdam is a breast cancer cohort from the Rotterdam tumour bank: 2982 patients and 1713 recurrence-free survival events (recurrence or death, whichever came first). The exposure is nodes, the number of positive axillary lymph nodes.

Three reasons for that choice, all of them properties of the data rather than conveniences:

Its relation to relapse really is non-linear — steep over the first few nodes, then flattening. That is precisely the shape a straight line cannot draw.
1436 of the 2982 patients have zero positive nodes (48.2%). That zero inflation breaks the textbook rule for placing knots, which is the next section.
Zero is a reference point clinicians recognise (node-negative disease), so the curve has an anchor the reader already understands.

library(survival)
library(splines)
data(cancer, package = "survival")

rot <- rotterdam
rot$rfs     <- pmax(rot$recur, rot$death)
rot$rfstime <- ifelse(rot$recur == 1, rot$rtime, rot$dtime)

# Linear version: one parameter, equal spacing forced
fit_lin <- coxph(Surv(rfstime, rfs) ~ nodes + age + size + meno, data = rot)

# Spline version, knots stated explicitly (why not the default: next section)
basis   <- ns(rot$nodes, knots = c(1, 4, 9), Boundary.knots = c(0, 20))
fit_spl <- coxph(Surv(rfstime, rfs) ~ basis + age + size + meno, data = rot)

# No exposure at all, for the test of the overall association
fit_non <- coxph(Surv(rfstime, rfs) ~ age + size + meno, data = rot)

anova(fit_non, fit_spl)   # overall association: is there a relationship?
anova(fit_lin, fit_spl)   # non-linearity: is that relationship a straight line?

Verified with R 4.6.0, survival 3.8.6 and splines 4.6.0. ns() ships with R itself; nothing extra to install.

import numpy as np
import pandas as pd
import statsmodels.api as sm
from patsy import dmatrix

# statsmodels' get_rdataset("rotterdam", "survival") fails: the Rdatasets index
# files it under survival's cancer bundle, so there is no standalone entry.
# Reading the CSV directly is the least trouble.
URL = "https://vincentarelbundock.github.io/Rdatasets/csv/survival/rotterdam.csv"
rot = pd.read_csv(URL)
rot["rfs"] = np.maximum(rot["recur"], rot["death"])
rot["rfstime"] = np.where(rot["recur"] == 1, rot["rtime"], rot["dtime"])

# constraints="center" is not optional. The unconstrained cr() basis sums to
# exactly 1 in every row, and a Cox model has no intercept, so that direction is
# not identifiable -- without it the nodes rows come back with nan for every
# standard error, and the whole inference half of the table is blank.
# lower_bound / upper_bound match R's Boundary.knots = c(0, 20).
SPL = 'cr(nodes, knots=[1, 4, 9], lower_bound=0, upper_bound=20, constraints="center")'

X = dmatrix(SPL + " + age + C(size) + meno",
            rot, return_type="dataframe").drop(columns="Intercept")
# ties="efron" matches R too: coxph() defaults to Efron, PHReg to Breslow.
fit = sm.PHReg(rot["rfstime"], X, status=rot["rfs"], ties="efron").fit()
print(fit.summary())

# The non-linearity test has to be assembled by hand: twice the difference
# in log-likelihood between the two models.
X_lin = dmatrix("nodes + age + C(size) + meno",
                rot, return_type="dataframe").drop(columns="Intercept")
fit_lin = sm.PHReg(rot["rfstime"], X_lin, status=rot["rfs"], ties="efron").fit()
lrt = 2 * (fit.llf - fit_lin.llf)   # 99.686 on 3 df

patsy's cr() is a natural cubic regression spline whose knot parameterisation differs from R's ns(), so the individual coefficients will not match, and they need not: the two are different bases for the same space. Align the knots, the boundary and the handling of ties and you are fitting the same model -- the non-linearity LRT above is 99.686 on 3 df, agreeing with R to six decimal places. constraints="center" is required rather than stylistic: without it the SE, t, p and confidence interval for every nodes row come back as nan. Every number on this page comes from the R above.

Where the knots go, and what a zero-inflated variable does to them

The most-cited rule comes from Harrell: three knots at the 10th, 50th and 90th percentiles; four at 5/35/65/95; five at 5/27.5/50/72.5/95. It is sound advice — put the knots where the data are dense, so every segment has enough people holding it up.

But it is advice about variables with a continuous distribution. Here is how nodes is distributed:

Percentile	5%	25%	50%	75%	95%
`nodes`	0	0	1	4	12

The two lowest percentiles are both zero, because nearly half the cohort (1436 patients, 48.2%) has no positive nodes. Placing knots by percentile therefore stacks several of them on one value — and that value happens to be the lower boundary. R responds by moving the knots itself, and says so:

What ns(nodes, df = 3) actually does here

The message R prints is shoving 'interior' knots matching boundary knots to inside (that is the English-locale wording; an R running in another language prints a translation, so the exact string will differ).

ns(x, df = 3) needs two interior knots and puts them at the tertiles — the 33.3rd and 66.7th percentiles, not the quartiles. In this data the lower tertile is itself 0, identical to the lower boundary, so R shoves it aside: the interior knots end up at 0.25, 2 and the boundary knots at 0, 34, the upper one dragged all the way out to the maximum of the data. The knots you thought were at the tertiles are not there; the upper boundary has been pulled away by a handful of extreme patients, and only 22 patients (0.74%) have more than 20 nodes at all.

This is a warning, not an error, and it appears once. In an R Markdown document or a pipeline it scrolls past easily, after which every figure and every p value is computed from the moved knots while the Methods section says “knots were placed at the tertiles”.

So this page states them explicitly: ns(nodes, knots = c(1, 4, 9), Boundary.knots = c(0, 20)). The interior knots are at 1, 4, 9 nodes — positions that carry clinical meaning for staging and that also sit where the data are dense — and the boundaries at 0 and 20. The upper boundary is deliberately not the maximum: beyond it there are too few patients, and anything drawn out there is extrapolation rather than estimation.

The reference point: the curve is always a comparison

The spline coefficients themselves cannot be read — basis1 through basis4 correspond to no clinical quantity. What can be read is the curve: holding the other covariates fixed, move $x$ from a reference value to some other value, and the risk multiplies by this much.

Every dose-response curve therefore needs a reference point chosen first. This page uses nodes = 0, node-negative disease, because clinicians recognise it. The median or the 10th percentile are the other common choices; which one you pick does not change the shape of the curve, only its vertical position, and it changes none of the p values.

Why the confidence band pinches to a point at the reference

Readers meeting an RCS figure in a paper often assume the band collapsing to zero width at the reference is a drawing error. It is correct, and the opposite would be the error.

The curve does not plot “the absolute risk at this $x$ ”. It plots “the risk at this $x$ relative to the reference”. The reference compared with itself is a ratio of exactly 1, with no uncertainty at all — it is a definition, not an estimate.

Technically: subtract the spline basis evaluated at the reference from the basis at each grid point to get $d(x)$ , then take the variance as $d(x)^\top V d(x)$ , where $V$ is the covariance matrix of the spline coefficients. At the reference $d = 0$ , so the variance is 0.

The common way to get this wrong is to take the standard errors from predict(type = "terms", se.fit = TRUE), draw the band from them, and shift the point estimate to the reference. Those standard errors are computed around the term’s own mean; shifting the estimate does not shift the uncertainty with it, so the band at the reference has a non-zero width — which amounts to claiming that a patient compared with himself has a hazard ratio somewhere in a 95% interval.

Dose-response curve. The horizontal axis is the number of positive axillary lymph nodes from 0 to 20; the vertical axis is the hazard ratio relative to 0 nodes on a logarithmic scale. The curve starts at 1 on the left and climbs steeply over the first few nodes — 1.23 at 1 node, 2.09 at 4, 3.39 at 9 — then clearly flattens, peaking near 3.84 at about 15 nodes and easing back to 3.55 at 20. A pale blue 95% confidence band wraps the curve, widening at the right where patients are few and closing to a single point at the reference on the left. A dashed horizontal line marks a hazard ratio of 1; a dotted vertical line and an open circle mark the reference at 0 nodes. Along the bottom runs a grey distribution strip: one thin vertical bar per integer node count, bar height being the number of patients. The scale is set by the largest non-zero group (1 node, 367 patients), so the bar at 0 is clipped and is drawn the same height as the bar at 1 even though it stands for 1436 patients, 48.2% of the cohort; both facts are labelled inside the figure. — The dose-response curve. The grey bars along the bottom are the distribution of nodes (bar height is the patient count, with the zero bar clipped — it is really 1436 patients), a reminder that the right-hand part of the curve rests on very few people, which is also where the band widens.Plotting script `figures/scripts/B2-05-splines.R`

Two p values, two questions, reported separately

This is what reviewers catch most often, and it is written into this site’s coverage requirements. A spline model supports two different likelihood ratio tests:

Test	Models compared	Question	χ²	df	p
Overall association (spline)	spline vs no `nodes` at all	Is `nodes` related to the outcome?	320.17	4	< 0.001
Overall association (linear)	linear `nodes` vs no `nodes`	Same question, straight lines only	220.49	1	< 0.001
Non-linearity	spline vs linear	Is that relationship a straight line?	99.69	3	< 0.001

The two questions are independent, and all four combinations occur:

Overall significant, non-linearity not — a relationship is there, and no departure from a straight line was detected. A linear term will do.
Both significant (the case on this page) — a relationship is there, and a straight line cannot describe it.
Overall not significant, non-linearity significant — uncommon but real. In a U-shaped relationship both ends rise while the middle falls, and a single linear coefficient averages them into nothing. Reporting only the linear model’s p value misses it entirely.
Neither significant — no association was detected in this study. That is not the same as saying there is none; it means this data did not provide enough evidence to detect one.

What the linear model actually gets wrong

Overlay the straight line the linear model was forced to draw, and the cost becomes visible.

Two lines on the same dose-response axes. The blue solid line is the spline curve with its pale blue 95% confidence band, climbing steeply and then flattening at around 15 nodes. The red dashed line is the straight line implied by the linear model, which on a logarithmic scale rises steadily from 1 at the reference to 4.37 at 20 nodes. The two separate early and are furthest apart in the middle of the range, near 7.6 nodes, where the straight line is only 0.57 times the curve; a short vertical segment on the figure marks the gap at 4 nodes, where the spline gives 2.09 and the line 1.34. The straight line keeps rising, crosses the spline curve near 20 nodes, and ends above the spline point estimate while still inside the confidence band. — Same data, same covariates. The only difference is whether nodes gets one parameter or 4. The straight line runs below the curve over most of the range and above it at the far end, where the fewest patients are.Plotting script `figures/scripts/B2-05-splines.R`

Positive nodes	Spline HR (95% CI)	HR the linear model claims
0	1 (reference)	1.00
1	1.23 (1.10–1.37)	1.08
2	1.49 (1.30–1.72)	1.16
4	2.09 (1.83–2.38)	1.34
9	3.39 (2.95–3.89)	1.94
15	3.84 (3.23–4.56)	3.02
20	3.55 (2.80–4.50)	4.37

Three things are worth stopping on:

The worst understatement is in the middle of the range, not at the start. The ratio of the line to the curve bottoms out near 7.6 nodes, where the linear model gives an HR only 0.57 times the spline’s; at 1 node the ratio is 0.88, a far smaller gap. The 4-node gap marked on the figure (spline 2.09, line 1.34) is where most patients are, not where the discrepancy is largest. Across the whole stretch from 1 to 15 nodes, the straight line lies outside the spline’s confidence band.
The far tail is overstated instead. At 20 nodes the linear model extrapolates to 4.37 while the spline gives 3.55 (2.80–4.50). The line is back inside the interval here, but only because so few patients remain that the interval has grown too wide to separate the two — a wide interval is not evidence that two models agree.
The slope eases off sharply after the first few nodes. From 0 to 4 nodes the HR climbs from 1 to 2.09; from 4 to 9 it reaches 3.39; and from 9 all the way to 15 it only reaches 3.84, after which it stops rising. Clinically that makes sense: past a certain nodal burden, a few more nodes carry little extra information. A linear term has no structural way to express that flattening — it can only keep climbing.

How this page connects

Wondering why “just cut it into quartiles” is usually worse — see categorising a continuous variable. The short version: categorising trades power for the linearity assumption, whereas a spline drops the assumption without paying that price.
What to do after a residual plot shows curvature — back to regression diagnostics.
The continuous-covariate assumption in a Cox model — see the Cox model. Note that this assumption (log hazard linear in the covariate) and the proportional hazards assumption are two different things; fixing the first with a spline does nothing for the second.
Splines also appear on calibration, in a different role: there the spline draws a flexible calibration curve of predicted against observed risk, and the exposure is the model’s own linear predictor.

Common misuses

Misuse	Why it is wrong
Using `ns(x, df = 3)` without reading the warning	With a zero-inflated or highly discrete variable the knots get moved, and the Methods section no longer matches the model
Reporting a single p value	“Is there a relationship” and “is it a straight line” are two questions, and reviewers ask
Calling the relationship linear because the non-linearity test was not significant	Only “no departure from linearity was detected” is supportable
Not stating the reference point	The vertical axis of the curve is then undefined and nobody can reproduce it
A confidence band with non-zero width at the reference	The standard errors are wrong, usually from `type = "terms"`
Drawing the curve out to the maximum of the data	The tail is extrapolation, its shape set by a handful of people
Using a high-order polynomial instead	Polynomials are global; an outlier at one end swings the other end
Adjusting knots until the curve looks right	That is choosing the model from the result; p values and intervals stop meaning anything
Interpreting the spline coefficients themselves	The basis functions carry no clinical meaning; only the curve can be read
Claiming the Cox assumptions are satisfied because the spline fixed the shape	Proportional hazards is a separate assumption and needs its own check

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B2-05-splines.R

Read the figure

The answer comes from the same statistical output that produced this page's figures, not from a number typed in beside them.

Entered linearly, each additional node carries a hazard ratio of 1.077. What does the spline's hazard ratio for nodes = 4 against nodes = 0 show?

Show the answer and why

Correct answer: 2.09 - well above what the linear model gives at the same point, so linearity understates this stretch

The spline gives 2.09 at nodes = 4 where the linear model gives only 1.34 - forcing one constant multiplier across the whole range understates the low end badly. 1.08 is that per-node multiplier: a slope rather than a contrast at a point, so it is not the same quantity at all. What makes this matter is that nearly half the patients have zero nodes, so the steepest part of the curve sits exactly where the people are.

Nonlinearity is tested by a likelihood ratio test. Which model does it compare the spline against, and what are the degrees of freedom?

Show the answer and why

Correct answer: Against the model with nodes entered linearly, 3 degrees of freedom

The nonlinearity test asks whether linearity is sufficient, so the comparison model is the linear one and the difference is 3 parameters. Comparing against the empty model gives 4 degrees of freedom, but that test asks whether nodes are associated with the outcome at all - a test that will almost certainly pass, so reading it as evidence of nonlinearity answers a different question with a foregone conclusion. 1 degree of freedom is the linear model against the empty one. Choosing the wrong comparison model is the standard misuse of this family of tests.

What happens if ns(nodes, df = 3) is used and R places the knots itself?

Show the answer and why

Correct answer: The first interior knot lands at 0.25 - knots go on quantiles, and with nearly half the sample at zero they bunch together

Automatic placement puts knots on quantiles, and with nearly half the patients at zero the quantiles collapse together: the first interior knot lands at 0.25, and R even warns that it has pushed interior knots away from the boundary. 4.00 is the upper-quartile node count rather than a knot, and 1.00 is the knot this page specifies by hand. Quantile placement is a reasonable default for a spread-out variable and fails on a pile of zeros - which is when the knots need choosing deliberately, or the zeros need thinking about separately.

Chapters that use this method

Sources and licences

This page is original writing