Royston-Parmar Model

Royston-Parmar Model for Survival Extrapolation

The Royston-Parmar model, sometimes called the RP spline model, is a parametric survival model used in health economic evaluation to fit trial survival data whose hazard changes shape over follow-up. It writes a transformation of the survival curve, most often the log of the cumulative hazard, as a smooth curve in log time, so the hazard can turn and change shape in ways that a single Weibull curve, whose hazard only rises or only falls, cannot follow. Patrick Royston and Mahesh Parmar published it in Statistics in Medicine in 2002, and its fitted curves supply the lifetime survival that a partitioned survival model or state transition model needs. Its weak point for decision making is the tail, because beyond the last knot the curve is a straight line on the transformed scale and the lifetime estimate rests on that line. Cure and competing-risk versions belong to the wider flexible parametric model family.

Origin as an extension of the Weibull and log-logistic models

Royston and Parmar noted that censored survival data were almost always analysed with Cox proportional hazards regression, although parametric models can sometimes handle non-proportional hazards simply and make the hazard much easier to visualise. They proposed extensions of the Weibull and log-logistic models in which natural cubic splines smooth the baseline log cumulative hazard and the log cumulative odds of failure, with further extensions for non-proportional covariate effects and a test of whether the chosen scale suits the covariate effects.

NICE Decision Support Unit Technical Support Document 21 (TSD21) states that flexible parametric models were developed because standard parametric models often could not capture the hazard shapes seen in applied studies. A natural cubic spline is the same object as a restricted cubic spline, joined smoothly at knots and constrained to be linear outside the boundary knots.

The model on the log cumulative hazard scale

The starting point is the Weibull model. TSD21 shows that a Weibull survival curve, transformed to the log cumulative hazard scale, 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 the shape parameter. The Royston-Parmar model replaces this straight line with a spline and adds covariates.

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

where $z$ is a vector of covariates such as treatment arm, $\beta$ holds their coefficients and $s(\ln t;\gamma)$ is a natural cubic spline in log time with coefficients $\gamma$. TSD21 notes that this is a proportional hazards model, so exponentiating a coefficient in $\beta$ gives a hazard ratio. The parameterisation used by Royston and Parmar, as set out in the documentation of the R package flexsurv, is

$$s(x;\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$$

where $x$ is log time, $m$ is the number of internal knots, $k_1,\dots,k_m$ are the internal knots, $k_{\min}$ and $k_{\max}$ are the boundary knots, $\lambda_j = (k_{\max}-k_j)/(k_{\max}-k_{\min})$, and $(u)_+$ equals $u$ when $u$ is positive and zero otherwise. With $m = 0$ only $\gamma_0$ and $\gamma_1$ remain and the model is a Weibull model with shape $\gamma_1$. Writing $\eta(t)$ for the right-hand side, survival and hazard follow directly.

$$S(t \mid z) = \exp\{-\exp(\eta(t))\}, \qquad h(t \mid z) = \frac{1}{t}\,\frac{ds}{dx}\,\exp(\eta(t))$$

where $S(t \mid z)$ is the probability of being event-free beyond $t$, $h(t \mid z)$ is the hazard, and $ds/dx$ is the slope of the spline at $x = \ln t$. The hazard formula holds when covariate effects do not vary with time. Parameters are estimated by maximum likelihood.

Scales, knots and time-dependent effects

The log cumulative hazard is one of three scales in common use. Each fixes what a time-constant covariate effect means and which distribution the model reduces to without internal knots, and NICE DSU Technical Support Document 14 (TSD14) likens these models to generalised linear models with link functions.

Knot choice sets how flexible the curve is. TSD21 states that current software usually places knots evenly along the distribution of uncensored log event times, with the boundary knots at the shortest and longest of those times, so that five knots sit at the 0th, 25th, 50th, 75th and 100th centiles. Predicted survival within follow-up is very insensitive to the number and position of knots, provided there are enough, while the Akaike and Bayesian information criteria give an informal guide to the number of knots that says little about extrapolation.

Proportional hazards can be relaxed by interacting a covariate with a second spline in log time, which gives a hazard ratio that changes over follow-up. TSD21 notes that the second spline can have fewer knots than the baseline, since the shape of the hazard is usually more complex than the deviation between two hazards. The flexsurv documentation reports that, in the breast cancer data of the 2002 paper, a proportional hazards spline with one internal knot gave log hazard ratios of 0.835 and 1.612 for two prognostic groups, against 0.840 and 1.618 from a Cox model.

Worked example with one internal knot

The numbers in this example are illustrative. A cancer trial has five years of follow-up, with uncensored event times running between 0.5 and 5 years and a median of 2 years. A Royston-Parmar model on the log cumulative hazard scale is fitted with boundary knots at $\ln 0.5 = -0.693$ and $\ln 5 = 1.609$ and one internal knot at $\ln 2 = 0.693$. The fitted coefficients are $\gamma_0 = -1.8$, $\gamma_1 = 1.4$ and $\gamma_2 = 0.18$, with a treatment coefficient of $\beta = -0.35$.

1. Knot weight. The weight for the internal knot is

$$\lambda_1 = \frac{1.609 - 0.693}{1.609 + 0.693} = \frac{0.916}{2.302} = 0.3979$$

where the numerator is $k_{\max} - k_1$ and the denominator is $k_{\max} - k_{\min}$.

2. Basis function at 3 years. Here $x = \ln 3 = 1.0986$, which lies between the internal knot and the upper boundary knot, so the last term of $v_1$ is zero.

$$v_1 = 0.4055^3 - 0.3979 \times 1.7918^3 = 0.0667 - 2.2889 = -2.2222$$

where 0.4055 is $x - k_1$ and 1.7918 is $x - k_{\min}$.

3. Survival for the comparator arm at 3 years. The spline value is

$$s = -1.8 + 1.4 \times 1.0986 + 0.18 \times (-2.2222) = -0.6620$$

where the three terms are $\gamma_0$, $\gamma_1 x$ and $\gamma_2 v_1$. The cumulative hazard is $H = e^{-0.6620} = 0.5158$ and survival is $S = e^{-0.5158} = 0.597$.

4. Hazard at 3 years. Differentiating $v_1$ gives its slope of $-3.3391$, so the spline slope and the hazard are

$$1.4 + 0.18 \times (-3.3391) = 0.799, \qquad 0.799 \times 0.5158 / 3 = 0.137$$

where the second product is the spline slope times $H$ divided by $t$, a hazard of 0.137 per year.

5. Treatment arm. The hazard ratio is $e^{-0.35} = 0.705$, so treated survival at 3 years is $S = 0.597^{0.705} = 0.695$.

Repeating steps 2 to 5 at other times gives these curves. The hazard rises to a peak near one year and then falls, a shape no single Weibull curve can produce.

6. What drives the extrapolation. Beyond the last knot at 5 years the cubic and squared terms of $v_1$ cancel, and the spline becomes a straight line in log time with slope

$$1.4 + 0.18 \times (-3.811) = 0.714$$

where $-3.811$ is the constant slope of $v_1$ beyond $k_{\max}$. After 5 years the comparator hazard follows a Weibull shape with shape parameter 0.714 and keeps falling. A hypothetical tail that kept these curves up to 5 years and then continued with a slope of 1.0 would give 20-year survival of 0.049 rather than 0.132. By numerical integration, the comparator's restricted mean survival time over 20 years is 7.24 years with the first tail and 6.20 years with the second, a difference of 1.04 life years from two curves that are the same within the trial.

Extrapolation beyond the last knot

TSD21 states that projection beyond the data is set by the linear term beyond the final knot. Without external data a Royston-Parmar extrapolation rests entirely on this linearity, and with many knots it may follow the trend at the end of follow-up, where events are few. A good fit within the data can still produce a plateau lasting many years, with a hazard below that of the general population.

The TSD21 simulation study bears this out. Flexible parametric models had very low bias for restricted mean survival at 3 years across all scenarios once they had at least 3 degrees of freedom. Lifetime mean survival was another matter, with high bias in some scenarios when survival was long or unobserved heterogeneity large, and no less bias than standard models when a cure fraction was present. Fitting the model as a relative survival model, so that background mortality is part of the hazard, reduced bias appreciably in some scenarios. TSD21 recommends that background mortality and other relevant external information should always be considered for incorporation.

Use in NICE technology appraisal

TSD14, published in 2011 and last updated in March 2013, said that Royston and Parmar's spline-based models were potentially very useful, although none had yet been used in NICE appraisals, and advised considering them when standard models fit poorly or extrapolate implausibly. By 2020 TSD21 could cite appraisals that used them. In TA498 the evidence review group used a spline model for treatment duration, and in TA517 the company used spline models for progression-free and overall survival.

TSD21 also records a misreading. In TA483 the company, evidence review group and committee all read a model with 2 knots as implying 3 heterogeneous subgroups of patients, whereas such a model only describes a complex hazard. The NICE technology appraisal manual (PMG36) now names spline-based approaches among the more flexible options for extrapolating time-to-event data. Section 4.6.25 asks for the clinical plausibility of the hazard functions and the uncertainty in the extrapolated portion to be assessed, and section 4.6.24 points to TSD21 when the proportional hazards assumption does not hold in some of the studies being synthesised.

How it differs from neighbouring survival models

The Royston-Parmar model sits between the Cox model, which leaves the baseline hazard unspecified, and the standard distributions, which fix its shape. What matters for a cost-effectiveness model is how each behaves after follow-up ends. The Cox model gives no estimates beyond the last event and needs a separate baseline curve. The Royston-Parmar model follows its final linear segment, whose slope is set by late events and knot choice, while a cure model adds a cured fraction that requires a credible cure assumption.

A spline model in general can bend any regression relationship, such as cost against age. The Royston-Parmar model is the specific use of a natural cubic spline in log time on a survival scale.

Misreadings that recur in appraisals

Most errors come from treating a good fit inside follow-up as evidence about the tail. Four recur.

• Reading the knots as subgroups of patients, as in TA483, when a knot is only a joining point of the curve.
• Choosing the number of knots by the lowest information criterion and then treating that model as the best extrapolation.
• Quoting one hazard ratio from a model with a time-dependent treatment effect, whose hazard ratio changes with time.
• Assuming the linear tail is safe, when on the log cumulative hazard scale it is a Weibull-shaped hazard that can fall below general population mortality.

Sources

• Jackson CH. flexsurv: a platform for parametric survival modeling in R. Journal of Statistical Software. 2016;70(8). https://doi.org/10.18637/jss.v070.i08
• Latimer NR. NICE DSU Technical Support Document 14: Survival analysis for economic evaluations alongside clinical trials, extrapolation with patient-level data. Sheffield: Decision Support Unit, ScHARR, University of Sheffield; June 2011, last updated March 2013. https://www.sheffield.ac.uk/sites/default/files/2022-02/TSD14-Survival-analysis.updated-March-2013.v2.pdf
• National Institute for Health and Care Excellence. NICE technology appraisal and highly specialised technologies guidance: the manual (PMG36), sections 4.6.24 and 4.6.25. London: NICE; published 31 January 2022, last updated 13 October 2026. https://www.nice.org.uk/process/pmg36/chapter/economic-evaluation-2
• 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; 23 January 2020. https://www.sheffield.ac.uk/sites/default/files/2022-02/TSD21-Flex-Surv-TSD-21_Final_alt_text.pdf