Cubic Spline

Cubic Spline in Survival Modelling and Economic Evaluation

A cubic spline joins short cubic polynomials end to end at chosen points called knots, and the joins are made smooth so that the curve and its first two derivatives agree wherever two pieces meet. In health economics its best-known use is inside survival models, where it lets a fitted survival or hazard curve bend to follow trial data that no single standard distribution matches. The choice matters because a curve fitted to a few years of follow-up is usually carried over a lifetime horizon, and the behaviour of the spline beyond its last knot then shapes mean survival, life years and the cost per QALY. Spline-based models are among the flexible options the NICE manual names for extrapolating survival, and NICE Decision Support Unit Technical Support Document 21 (TSD 21) describes how they work and where they mislead.

How cubic pieces are joined at knots

Between two neighbouring knots the curve is an ordinary cubic polynomial. A different cubic is used on the next interval. With knots $k_1<k_2<\dots<k_m$, the piece on the interval from $k_j$ to $k_{j+1}$ can be written as

$$s_j(x)=a_j+b_j(x-k_j)+c_j(x-k_j)^2+d_j(x-k_j)^3$$

where $x$ is the input variable (in survival work usually log time), $k_j$ is the knot at the left end of the interval and $a_j$, $b_j$, $c_j$ and $d_j$ are the coefficients of that piece.

In TSD 21's description the function has continuous first and second derivatives at the knots, so the gradient and the rate of change of the gradient agree where pieces meet. Only the third derivative may jump at a knot, which is why a cubic spline can change shape at each knot without a visible corner. More knots therefore allow a more complex shape.

Counting parameters shows what the constraints remove. Each cubic piece has four coefficients and each interior knot imposes three equality conditions, so a cubic spline with $K$ interior knots and $K+1$ pieces has $p=4(K+1)-3K=K+4$ free parameters, the count the R bs function returns ($K+3$ basis columns plus an intercept).

Natural and restricted cubic splines

An unrestricted cubic spline can swing sharply beyond its outermost knots, because each end piece is a full cubic with no data on its outer side. A natural cubic spline adds the constraint that the function is linear beyond the boundary knots, which is how the R splines documentation defines its natural spline basis. In regression the same object is usually called a restricted cubic spline. Durrleman and Simon described restricted cubic spline regression for medical data, comparing it with non-parametric methods for the relationship between age and survival in the Stanford Heart Transplant data.

Treating the two boundary points as extra knots gives $K+6$ parameters, and linearity in both tails removes four of them, so a natural spline with $K$ interior knots has $K+2$ parameters including the intercept. The basis used in survival software, documented for the R flexsurv package, writes the spline as

$$s(x;\boldsymbol{\gamma})=\gamma_0+\gamma_1 x+\gamma_2 v_1(x)+\dots+\gamma_{m+1}v_m(x), \qquad v_j(x)=(x-k_j)_+^3-\lambda_j(x-k_{\min})_+^3-(1-\lambda_j)(x-k_{\max})_+^3, \qquad \lambda_j=\frac{k_{\max}-k_j}{k_{\max}-k_{\min}}$$

where $k_{\min}$ and $k_{\max}$ are the boundary knots, $k_1$ to $k_m$ are the interior knots, $\gamma_0$ to $\gamma_{m+1}$ are the coefficients and $(u)_+$ equals $u$ when $u$ is positive and 0 otherwise. Below $k_{\min}$ every basis function is zero. Above $k_{\max}$ the cubic and quadratic terms in each $v_j$ cancel, so the function is a straight line at both ends. The restricted cubic spline page covers regression uses, interpolation covers splines drawn exactly through given points, and curve fitting covers fitting curves to noisy data in general.

Royston-Parmar splines on the log cumulative hazard scale

The main health economic use of cubic splines comes from Royston and Parmar, who extended the Weibull and log-logistic models by using natural cubic splines to smooth the baseline log cumulative hazard and the log cumulative odds. TSD 21 starts from the Weibull model. Its log cumulative hazard is a straight line in log time,

$$\ln H(t)=\ln\lambda+\gamma\ln t$$

where $H(t)$ is the cumulative hazard at time $t$, $\lambda$ is the Weibull scale parameter and $\gamma$ is its shape parameter. A flexible parametric model replaces the straight line with a restricted cubic spline in log time,

$$\ln H(t\mid\mathbf{z})=s(\ln t;\boldsymbol{\gamma})+\boldsymbol{\beta}^{\mathsf T}\mathbf{z}$$

where $s$ is the spline above evaluated at $x=\ln t$, $\mathbf{z}$ holds covariates such as treatment arm and $\boldsymbol{\beta}$ holds their coefficients. Survival follows as $S(t)=\exp\{-H(t)\}$. With no interior knots the spline is a straight line and the model is the Weibull model. TSD 21 notes that covariate effects in this form are log hazard ratios, so exponentiated coefficients are hazard ratios. Logistic and probit transformations of survival give spline extensions of the log-logistic and log-normal models in the same way. The Royston-Parmar model page sets out the family in full.

Proportional hazards can be relaxed by interacting treatment with a second spline function of log time. TSD 21 observes that this second spline often needs fewer knots than the baseline, because the difference between two hazard functions is usually simpler than the hazard itself. Models can also be fitted on the log hazard scale, which needs numerical integration during estimation but may make assumptions about a treatment effect after follow-up easier to judge.

Knots and degrees of freedom in survival software

Unless the analyst supplies knot locations, software applies a default rule. TSD 21 reports that current implementations put the boundary knots at the minimum and maximum uncensored log event times and spread interior knots evenly over their distribution, so a five-knot model has knots at the 0th, 25th, 50th, 75th and 100th centiles. In the R flexsurv package the number of interior knots equals the degrees of freedom minus one in the notation of the Stata command stpm, and zero interior knots gives a Weibull, log-logistic or log-normal model depending on the scale.

Within follow-up the choice is forgiving, and TSD 21 cites evidence that predicted survival inside the observed period is very insensitive to the number and location of knots once there are enough to capture the shape, while asking for attention to the number of events in trial data, since many spline parameters may be impractical when events are few. The Akaike and Bayesian information criteria are often used as an informal guide to the number of knots, but TSD 21 adds that they say little about how well a model extrapolates, a point developed on the Akaike information criterion page.

Worked example with one interior knot

An illustrative model with hypothetical numbers shows how the basis is built and why the tail behaves as it does. Time is measured in months and the knots sit on the log time scale at $k_{\min}=0$, $k_1=1$ and $k_{\max}=3$, which correspond to about 1, 2.7 and 20.1 months. The last knot stands for the end of trial follow-up, and the coefficients are $\gamma_0=-3$, $\gamma_1=0.8$ and $\gamma_2=0.05$.

1. Basis weight. The single interior knot gives $\lambda_1=(3-1)/(3-0)=0.6667$, so the basis function is $v_1(x)=(x-1)_+^3-0.6667x_+^3-0.3333(x-3)_+^3$.
2. Survival at 12 months. Here $x=\ln 12=2.4849$. The basis value is $v_1=1.4849^3-0.6667\times15.3436=3.2741-10.2291=-6.9550$, so the spline is $s=-3+1.9879-0.3478=-1.3599$. The cumulative hazard is $e^{-1.3599}=0.2567$ and survival is $e^{-0.2567}=0.7736$, so about 77% of patients are alive.
3. Survival at the last knot. At $x=3$ the basis value is $v_1=8-18=-10$, so the spline is $s=-3+2.4-0.5=-1.1$, the cumulative hazard is $e^{-1.1}=0.3329$ and survival is 0.7169.
4. The tail beyond the last knot. For $x$ above 3 the derivative $v_1'(x)=3(x-1)^2-2x^2-(x-3)^2$ equals $-6$ at every point. The slope of the log cumulative hazard is therefore fixed at $s'=0.8-0.3=0.5$, and the extrapolated tail is a Weibull curve with shape 0.5, $H(t)=0.3329(t/20.09)^{0.5}$. At 60 months the cumulative hazard is $H=0.3329\times1.7284=0.5754$ and survival is 0.5625.
5. The hazard. Differentiating $H(t)=\exp\{s(\ln t)\}$ gives $h(t)=H(t)\,s'(\ln t)/t$. The monthly hazard is $h=0.3329\times0.5/20.09=0.0083$ at the last knot and $h=0.5754\times0.5/60=0.0048$ at 60 months.

Inside follow-up the spline lets the slope fall from 0.8 at the first boundary knot to 0.5 at the last, a bend that a single Weibull curve cannot make. Beyond the last knot the model has no further flexibility, so the hazard keeps falling in proportion to $t^{-0.5}$ for as long as the model runs. In a cohort that is ageing, that assumption would need checking against general population mortality before the curve could inform mean survival.

What a spline tail assumes when extrapolated

TSD 21 states that projection beyond the data is dictated by the linear term beyond the final knot on the transformed scale, which may or may not be sensible. With many knots that line can follow the trend at the very end of follow-up, which may rest on few events, and changing the number of knots can alter the extrapolation even when the fit within the data barely moves. A good fit within the range of the data may still give implausible extrapolations, such as a long plateau whose implied hazard sits below that of the general population. In one TSD 21 example every fitted hazard was decreasing at the end of follow-up, which the authors say should cause alarm because the true hazard of an ageing cohort would be expected to start rising.

The NICE manual (PMG36, section 4.6.25) lists spline-based approaches among the more flexible options for extrapolation of time-to-event data. It asks that the clinical plausibility of each model's underlying hazard function be routinely assessed and that uncertainty in the extrapolated portion be explored. TSD 21 recommends plotting extrapolated survival and hazard functions together with the expected mortality rate for any extrapolation. Spline models can also carry background mortality in a relative survival framework, and in a flexible parametric cure model the placement of the final boundary knot sets the time from which cure is assumed.

Where cubic splines sit among survival curve options

A spline model is one of several ways to escape the fixed shapes of a standard parametric survival model. For complex hazards TSD 21 lists flexible parametric models with splines or fractional polynomials, models that impose a cure proportion, more general mixture models and conditional approaches such as piecewise models and landmarking. A spline describes the shape of one hazard function, whereas a mixture cure model treats the population as two groups, cured and uncured, with separate survival.

When splines are used in a partitioned survival model, the extrapolated curves for progression-free and overall survival set the time spent in each health state, so the spline assumptions pass straight into costs and QALYs.

Spline models in NICE appraisals and how they are misread

TSD 21 records spline models in NICE appraisals. In TA498 (lenvatinib with everolimus for previously treated advanced renal cell carcinoma) the evidence review group used a flexible parametric model with splines for treatment duration. In TA517 (avelumab for metastatic Merkel cell carcinoma) the company used spline models for progression-free and overall survival.

The same document reports that in TA483 the company, the evidence review group and the appraisal committee all read a model with 2 knots as implying three heterogeneous subgroups of patients. TSD 21 explains that such a model simply represents a complex hazard function and makes no direct assumption about subgroups. Three further misreadings are easy to make.

• Natural taken to mean straight. The spline is linear only beyond its boundary knots and only on the transformed scale. On the log cumulative hazard scale against log time, a straight tail is a Weibull shape, so the hazard can keep rising, or keep falling towards zero, for as long as the model runs.
• Coefficients read as effects. The spline coefficients depend on the basis chosen, and flexsurv offers a second, orthogonal natural cubic spline basis for the same model. A report needs the scale, the knot locations and the fitted curves, since coefficients from different bases are not comparable.
• Fit taken as evidence for the tail. A model can match the Kaplan-Meier curve closely and still extrapolate badly, because information criteria and visual fit say little about the period beyond follow-up.

Sources

• Durrleman S, Simon R. Flexible regression models with cubic splines. Statistics in Medicine. 1989;8(5):551-561. https://doi.org/10.1002/sim.4780080504
• Jackson C. flexsurvspline: flexible survival regression using the Royston/Parmar spline model. flexsurv R package documentation, version 2.3 (accessed 11 October 2026). https://search.r-project.org/CRAN/refmans/flexsurv/html/flexsurvspline.html
• National Institute for Health and Care Excellence. NICE technology appraisal and highly specialised technologies guidance: the manual (PMG36). Published 31 January 2022, last updated 13 October 2026. Section 4.6.25. https://www.nice.org.uk/process/pmg36/chapter/economic-evaluation-2
• R Core Team. ns and bs: natural cubic spline and B-spline basis functions. splines package documentation, R version 4.6.1 (accessed 11 October 2026). https://search.r-project.org/R/refmans/splines/html/ns.html
• Royston P, Parmar MKB. Flexible parametric proportional-hazards and proportional-odds models for censored survival data, with application to prognostic modelling and estimation of treatment effects. Statistics in Medicine. 2002;21(15):2175-2197. https://doi.org/10.1002/sim.1203
• Rutherford MJ, Lambert PC, Sweeting MJ, Pennington B, Crowther MJ, Abrams KR, Latimer NR. NICE DSU Technical Support Document 21: Flexible methods for survival analysis. Sheffield: Decision Support Unit, ScHARR, University of Sheffield; 2020 (updated March 2022). Executive summary and sections 2.1.3, 3.2, 3.6, 5.2.3 and 6.3. https://doi.org/10.15131/shef.data.33849118