Standard parametric survival models, such as the Weibull or Gompertz distributions, are often insufficient for capturing complex, non-monotonic hazard functions. When hazards rise and fall, or exhibit other intricate patterns, a more flexible approach is required. This is where flexible parametric spline models, notably Royston-Parmar models (RPMs), become indispensable. RPMs utilise restricted cubic splines to model a transformation of the survival function, typically on the log-cumulative hazard scale.

The core idea is to model a function $g(S(t,z))$ as a spline function of log time, s(x | oldsymbol{eta}, oldsymbol{\gamma}, k), where $x = \ln(t)$. The choice of transformation, $g(.)$, is critical. Modelling on the log-cumulative hazard scale, $g(S(t)) = \ln(-\ln(S(t)))$, is the most common approach. This directly ensures that the estimated survival function $S(t)$ is always between 0 and 1, and the hazard function is always positive. Alternatively, one can model on the log-hazard scale, $g(S(t))=\ln(h(t))$. While this can be more intuitive for interpreting covariate effects on the hazard itself, it does not automatically constrain the survival function to be valid. The log-cumulative hazard approach is generally more robust for extrapolation.

The flexibility of these models is governed by the number and position of knots in the spline function. Knots are the points where the polynomial pieces of the spline join. Using too few knots can lead to underfitting, failing to capture the true shape of the hazard. Conversely, too many knots can lead to overfitting, where the model captures random noise in the data rather than the underlying signal. This results in unstable and biologically implausible extrapolations.

Knot placement is typically based on the quantiles of the log-event times, ensuring the model has flexibility where most events occur. For example, a common strategy is to place knots at the minimum and maximum log-event times, with internal knots at specific centiles (e.g., 33rd and 67th). The degrees of freedom for the spline function are equal to the number of knots minus one.

For transformative therapies where a subset of patients may be considered 'cured', standard survival models are inadequate because they assume all patients remain at risk of the event. Cure models explicitly account for this by modelling a 'cure fraction' — a proportion of the population that is not expected to experience the event. There are two primary classes of cure models: Mixture Cure Models (MCM) and Non-Mixture Cure Models (NMCM).

MCMs assume the population is a mix of two sub-populations: a 'cured' group with a survival probability of 1, and a 'susceptible' group whose time to event is described by a standard survival function, $S_u(t)$. The model estimates both the cure fraction, $\pi$, and the parameters of the survival function for the susceptible group.

NMCMs, or bounded hazard models, provide an alternative. Instead of assuming a mixture, they model the survival of the whole population using a function that eventually reaches a non-zero plateau. This plateau represents the statistical cure fraction. The hazard function in an NMCM approaches zero over time, whereas in an MCM, the hazard for the susceptible group may not. These models are often used when the mechanism of 'cure' is thought to be a long-term remission rather than a complete biological cure from the outset.

A major challenge in health technology assessment, particularly for submissions to bodies like NICE, is ensuring that long-term survival extrapolations are clinically and biologically plausible. A common pitfall is a model predicting that patients with a specific condition will live longer than the general population, which is usually unrealistic. This requires incorporating background mortality.

The NICE DSU Technical Support Document 21 provides guidance on this. The recommended approach involves using external data, typically national life tables, to represent general population mortality. The hazard rate for the modelled population, $h_p(t)$, should be constrained such that it does not fall below the hazard rate of the matched general population, $h_g(t)$.

In practice, this often means that after a certain point in time, the model's hazard rate is fixed to the general population rate. For a relative survival model, this is the point where the relative survival function plateaus at a constant value, implying the excess hazard has fallen to zero. For other models, it may involve a piecewise approach where the extrapolated hazard, $h(t)$, is defined as $max(h_p(t), h_g(t))$. This ensures the model respects the fundamental biological constraint imposed by background mortality, leading to more credible long-term cost-effectiveness estimates.

This integration of external mortality data is crucial for validating the long-term behaviour of any complex survival model, from flexible splines to cure models.