Survival Extrapolation for HTA Using RWE

In plain language

Survival extrapolation is how health economists turn a short stretch of trial data into a full lifetime survival curve so they can estimate the total life-years a treatment provides. A clinical trial might only follow patients for two to three years, but a cost-effectiveness model needs to project what happens over an entire lifetime — sometimes 30 years. To do that, analysts fit a mathematical curve (called a parametric survival model) to the trial's observed data and then extend it far beyond what was actually seen. The critical catch is that several different curve shapes may all fit the observed data equally well, yet each one predicts a very different number of life-years in the unobserved future — and that difference can change whether a treatment looks cost-effective.

Survival extrapolation for HTA

is the step that turns a finite stretch of observed follow-up (a trial, registry, or claims cohort with, say, 18-36 months of data) into the lifetime survival curve a cost-effectiveness model requires. Health technology assessment is almost always conducted over a lifetime horizon, but time-to-event data are administratively censored at the data cut, so the analyst must fit a model to the observed Kaplan-Meier (KM) curve and project the hazard forward for years or decades. Because the incremental cost-effectiveness ratio (ICER) is dominated by the area under the survival curve (mean survival / life-years gained), and because the bulk of that area frequently lies beyond the observed data, the choice of extrapolation model is often the single most influential and most contested assumption in the entire appraisal. NICE DSU TSD 14 (Latimer) codified the standard workflow; TSD 21 extended it to flexible parametric and externally-informed approaches.

Core conceptual distinction

Survival extrapolation is not "survival analysis." Within-data estimands such as the hazard ratio, restricted mean survival time (RMST) over the observed window, or a Kaplan-Meier point estimate are interpolations governed by observed events. Extrapolation is an out-of-sample projection: the quantity of interest is mean (often discounted) survival to a lifetime horizon, which depends entirely on the assumed hazard shape in the unobserved tail where there are zero events to discipline the fit. Two models that are visually indistinguishable over the trial period - and have nearly identical AIC/BIC - can differ by years in projected mean survival because their tail hazards diverge (e.g., a log-normal with a long, decreasing hazard tail vs a Weibull with a monotone increasing hazard). The deliverable is therefore not a "best-fitting model" but a defensible projected hazard, justified by in-sample fit (AIC/BIC, visual, residuals) AND external plausibility (smoothed hazard shape, general-population mortality as a floor on the hazard, clinical expectation, and registry/long-term external data where available). The candidate set is the six standard parametric distributions (exponential, Weibull, Gompertz, log-normal, log-logistic, generalized gamma) plus flexible alternatives - Royston-Parmar restricted cubic splines, and mixture / non-mixture cure models when a plateau implies long-term survivors.

Pros, cons, and trade-offs

- Standard six parametric distributions (TSD 14) vs flexible parametric splines (TSD 21): the six closed-form distributions are transparent, fast, easy to do probabilistic sensitivity analysis (PSA) on (multivariate-normal draws of the parameters), and familiar to reviewers. Cost: each imposes a rigid, usually monotone or single-turning-point hazard, so when the true hazard is non-monotone (e.g., early treatment-related mortality then a plateau) none fits well and the choice among poor fits drives the answer. Royston-Parmar splines flex to the observed hazard and can be anchored to external data, but they are less constrained in the tail and can extrapolate implausibly if knots sit near the data boundary - more knots improve in-sample fit while worsening tail behavior. Prefer the six when the smoothed hazard is plausibly monotone and follow-up is mature; prefer splines or externally-anchored models when the hazard is complex or the tail is sparse. - Independent per-arm fits vs a single joint model with a treatment covariate (proportional hazards / accelerated failure time): fitting each arm separately maximizes within-arm fit but lets the two extrapolations cross or diverge implausibly and ignores the question of whether the treatment effect persists. A joint model (or a relative-effect-on-a-reference-curve approach) enforces a coherent relationship but bakes in a PH/AFT assumption that may be false in the tail. Prefer independent fits only when proportionality is clearly violated AND you separately justify the long-term effect; otherwise model the treatment effect explicitly so you can test waning. - Naive extrapolation vs general-population mortality flooring / relative-survival framing: ignoring background mortality lets a fitted curve project survival probabilities that exceed the age-matched general population - clinically impossible and a frequent ERG/EAG critique. Capping the all-cause hazard at the general-population life-table hazard (or modeling excess/relative survival) is more defensible but requires linkage to national life tables and an assumption about excess hazard in the tail. Prefer flooring/relative survival in any chronic or oncology model run to a lifetime horizon, especially in older populations. - vs RMST or within-trial analysis: RMST avoids the extrapolation problem entirely by restricting to the observed horizon, but it cannot answer the lifetime question HTA demands. Use RMST as a face-validity anchor (the model's restricted mean over the observed window should match the empirical RMST), not as a substitute.

When to use

Any cost-utility or cost-effectiveness analysis with a lifetime (or long, e.g., >5-year) horizon where survival or time-to-progression is a model input and follow-up is shorter than the horizon - i.e., essentially every oncology and most chronic-disease submissions to NICE, CDA-AMC (formerly CADTH), PBAC, IQWiG, ICER. Use it to populate partitioned-survival models (overall and progression-free survival curves) and to derive transition probabilities or dwell times for Markov / state-transition models. It is also the right tool when RWE (registry or claims) provides longer-term follow-up that can validate or directly inform the trial-based tail.

When NOT to use - and when it is actively misleading or dangerous

- The horizon does not exceed observed follow-up. If the decision horizon is fully covered by data (e.g., an acute condition resolved within the trial), extrapolation adds assumption-driven uncertainty for no benefit - report the empirical KM / RMST instead. - Few events / immature data with no external anchor. Extrapolating from a curve with little tail information (a handful of late events, wide late KM confidence bands) produces projections driven almost entirely by the parametric form, not the data; different equally-plausible distributions can imply mean survival differences large enough to flip the ICER across the willingness-to-pay threshold. Without external long-term data or expert elicitation this is conjecture dressed as analysis. - A plateau is fit with a non-cure distribution (or vice versa). Forcing a standard distribution onto immunotherapy-style data with a long-term survivor fraction will either understate the plateau (monotone increasing hazard) or invent immortality; conversely, fitting a cure model when the plateau is an artifact of administrative censoring (everyone simply ran out of follow-up at the data cut) fabricates cured patients. The apparent plateau at the right edge of a KM curve is frequently censoring, not biology. - Tail survival exceeds the general population. If the projected all-cause survival is better than the age/sex-matched life table, the model is impossible and must be floored - presenting it unfloored is misleading. - Assuming a treatment effect persists for life without evidence. Extrapolating a within-trial hazard ratio indefinitely (no waning) is a common and consequential optimistic bias; the persistence of the effect beyond the data must be an explicit, tested scenario, not a silent default.

Data-source operational depth

- Trial / IPD (the TSD 14 base case): cleanest hazards but shortest follow-up; the tail is almost entirely unobserved, so the danger is over-reading an apparent late plateau that is really administrative censoring at the data cut. Always overlay the number-at-risk and a smoothed hazard plot before choosing a distribution; a hazard estimated from <10 events is noise. - Registry: the natural source of long-term external data to anchor or validate the trial tail, but registries are themselves administratively censored at extraction and often have differential loss to follow-up (sicker patients drop out or die undocumented), which can create a spurious late survival improvement (the survivors look healthy because the unwell are missing). Confirm completeness of vital status (link to a death index) before trusting registry tails. - Claims: can give long real-world follow-up but the failure modes bite hard for extrapolation. Medicare Advantage (MA)-only person-time lacks fee-for-service (FFS) claims, so death and late events are differentially unobserved and the tail survival is artifactually inflated unless MA-only spans are excluded or vital status is linked. Mortality in elderly claims is dominated by competing, non-disease causes, so a naive all-cause extrapolation conflates the disease hazard with background mortality - frame as relative survival or floor the hazard at the general-population life table. Immortal time inflates early survival if the index date is pinned at a landmark (e.g., a second prescription) rather than at treatment initiation. Outcome (death) capture lags and claim reversals distort the most recent months - hold out the final incomplete quarters from the fit. - Linked claims-registry-vital records: the ideal substrate (registry severity + claims duration + reliable mortality from the death index), but linkage selects the linkable subset and introduces order/fill/service-date discrepancies that must be reconciled before time-zero and event dates are set.

Worked claims example

Question: lifetime overall survival for a new first-line therapy vs standard of care in metastatic NSCLC, to populate a partitioned-survival cost-utility model run to a 30-year (lifetime) horizon, using a commercial + Medicare FFS claims database linked to the National Death Index, where observed follow-up is ~30 months. (1) Cohort and time zero: adults with a metastatic NSCLC diagnosis and a first qualifying systemic-therapy fill (`person_id`, `fill_date`, `days_supply`, `dx` codes); index_date = that first fill; require continuous A/B/D FFS enrollment for a 365-day baseline and exclude MA-only person-time so death and late events are observable. (2) Outcome and follow-up: death from the linked death index; censor at disenrollment (excluding MA switches handled as above), end of data, and the data-cut date; drop the final two incomplete quarters to avoid reversal/lag artifacts. (3) Plot KM with number-at-risk and a smoothed hazard; the hazard is non-monotone (early peak, then declining) so the standard six fit poorly. (4) Fit the six parametric distributions plus a 2-knot Royston-Parmar spline per arm; rank by AIC/BIC, overlay each fitted curve on the KM, and inspect projected hazards to 30 years. (5) Floor the all-cause hazard at the age/sex-matched US life-table hazard so projected survival never exceeds the general population. (6) Select the generalized gamma (best AIC and a clinically plausible declining-then-leveling hazard), confirm its restricted mean over 30 months matches the empirical RMST as a face-validity check, and compute mean discounted life-years per arm. (7) PSA: draw the fitted parameter vectors from their multivariate-normal sampling distribution (Cholesky of the covariance matrix), recompute lifetime survival per draw, and propagate to the ICER. (8) Scenario analyses on the distribution choice (Weibull, log-normal, spline), a treatment-effect-waning scenario (hazard ratio reverts to 1 after 5 years), and inclusion vs exclusion of the registry-derived long-term tail.

Interpreting the output

Six standard parametric distributions fitted to 30-month trial data produce projected mean survival ranging from 2.1 years (exponential) to 6.2 years (log-normal) — a three-fold spread from models with nearly identical AIC (within 4 points).

Formal interpretation. Each parametric family implies a distinct hazard shape in the unobserved tail beyond the data window, and that tail — not the observed 30 months — dominates the lifetime mean survival and therefore the ICER. AIC and BIC measure in-sample fit and carry no information about which extrapolated tail is correct; they cannot discriminate among the six distributions beyond the last observed event time. The log-normal and generalized-gamma families allow non-monotone hazards and project substantially longer survival, while the exponential assumes constant hazard and produces the shortest estimate. Selection among families must therefore be justified by external biological plausibility, registry or epidemiological long-term mortality data, or clinical expert elicitation — not goodness-of-fit statistics alone.

Practical interpretation. A three-fold difference in projected mean survival translates directly to a three-fold difference in estimated life-years gained and potentially shifts the cost-effectiveness estimate from below to well above the payer willingness-to-pay threshold. Health technology assessment bodies (NICE, ICER) require structural uncertainty sensitivity analyses across all plausible families, with the base-case selection defended in writing. The choice of extrapolation model is typically the single largest driver of ICER uncertainty in immature survival data.

Worked example

Scenario

A randomized trial of a new lung cancer therapy followed patients for 30 months. The trial team recorded whether each patient had died and when. Now a health economist needs to populate a cost-effectiveness model that runs to a 30-year lifetime horizon. The observed data only cover the first 30 months, so a parametric survival model must be fit to those data and then extrapolated forward for the remaining 27.5 years. Three common distributions — exponential, Weibull, and log-normal — all fit the 30-month observed data reasonably well, but they make very different assumptions about how the hazard behaves in the unobserved tail.

Dataset

Summary of observed 30-month trial data for the treatment arm (used to fit each parametric model). The three fitted distributions are then projected forward to 30 years.

Steps

• All three models are fit to the same 30-month trial data; their AIC scores are close (408-412), so no model clearly wins on in-sample fit alone.

• Each model is then projected forward to 30 years using its mathematical formula — this is the extrapolation step, where the curves diverge sharply because each assumes a different shape for the hazard in the unobserved tail.

• The exponential model assumes the hazard (risk of dying) stays constant forever, which is pessimistic for a cancer therapy; it projects only 2.1 years of mean survival.

• The Weibull model allows the hazard to change over time in a single direction; it projects 3.4 years of mean survival.

• The log-normal model assumes the hazard first rises then falls (a long, slow declining tail), which is optimistic; it projects 6.2 years of mean survival.

• The analyst must choose — or present all three as scenarios — because the cost-effectiveness ratio depends heavily on this choice: the difference between 2.1 and 6.2 life-years is large enough to flip the conclusion about whether the treatment is cost-effective.

Result

Projected mean survival ranges from 2.1 years (exponential) to 3.4 years (Weibull) to 6.2 years (log-normal) — a three-fold spread from the most pessimistic to the most optimistic model, all fit to identical 30-month data. The choice of extrapolation model, not the trial itself, drives the final cost-effectiveness answer.

Runnable example

import numpy as np
import pandas as pd
from lifelines import (WeibullFitter, LogNormalFitter, LogLogisticFitter,
                       GeneralizedGammaFitter, ExponentialFitter)

HORIZON_YEARS = 30.0   # lifetime horizon for the HTA model
DISCOUNT = 0.035       # annual discount rate on life-years (NICE reference case)
GRID = np.linspace(0, HORIZON_YEARS, int(HORIZON_YEARS * 12) + 1)  # monthly grid to the horizon (30y -> 361 points)

# lifelines supplies five of the six TSD-14 standard distributions; Gompertz has no lifelines
# fitter and is fit via flexsurv (R) instead (see the R implementation below).
CANDIDATES = {
    "exponential": ExponentialFitter,
    "weibull":     WeibullFitter,
    "lognormal":   LogNormalFitter,
    "loglogistic": LogLogisticFitter,
    "gengamma":    GeneralizedGammaFitter,
}

def discounted_life_years(times, surv, rate):
    # Trapezoidal area under S(t) with continuous discounting -> discounted mean survival.
    disc = np.exp(-rate * times)
    integrand = surv * disc
    return np.trapz(integrand, times)

def fit_and_extrapolate(df_arm: pd.DataFrame) -> pd.DataFrame:
    rows = []
    for name, Fitter in CANDIDATES.items():
        f = Fitter()
        f.fit(df_arm["time_years"], event_observed=df_arm["event"])
        surv = f.survival_function_at_times(GRID).values  # projected S(t) to the horizon
        rows.append({
            "model": name,
            "AIC": f.AIC_,
            "BIC": getattr(f, "BIC_", np.nan),
            "mean_LY": np.trapz(surv, GRID),               # undiscounted lifetime mean survival
            "disc_LY": discounted_life_years(GRID, surv, DISCOUNT),
        })
    return pd.DataFrame(rows).sort_values("AIC").reset_index(drop=True)

# Per-arm fits; compare disc_LY across distributions to expose tail-driven divergence.
results = {arm: fit_and_extrapolate(g) for arm, g in surv.groupby("arm")}
for arm, tbl in results.items():
    print(f"\n=== {arm} ===")
    print(tbl.to_string(index=False))
# Selection rule: lowest AIC/BIC AND a plausible projected hazard AND survival floored at the
# general-population life table (apply min(model_hazard, lifetable_hazard) before integrating).

library(flexsurv)
library(dplyr)

HORIZON  <- 30      # lifetime horizon (years)
DISCOUNT <- 0.035   # annual discount rate (NICE reference case)
GRID     <- seq(0, HORIZON, by = 1/12)

dists <- c("exp", "weibull", "gompertz", "lnorm", "llogis", "gengamma")

fit_one_arm <- function(d) {
  fits <- lapply(dists, function(dd)
    flexsurvreg(Surv(time_years, event) ~ 1, data = d, dist = dd))
  names(fits) <- dists
  # Royston-Parmar spline (2 internal knots) on the log-cumulative-hazard scale.
  fits[["spline_k2"]] <- flexsurvspline(Surv(time_years, event) ~ 1, data = d,
                                        k = 2, scale = "hazard")

  aic <- sapply(fits, AIC)

  # Mean (discounted) life-years = trapezoidal area under projected S(t) * discount factor.
  disc <- exp(-DISCOUNT * GRID)
  ly <- sapply(fits, function(f) {
    S <- summary(f, t = GRID, type = "survival", ci = FALSE)[[1]]$est
    c(mean_LY = sum(diff(GRID) * (head(S, -1) + tail(S, -1)) / 2),
      disc_LY = sum(diff(GRID) * (head(S * disc, -1) + tail(S * disc, -1)) / 2))
  })
  data.frame(model = names(fits), AIC = aic,
             mean_LY = ly["mean_LY", ], disc_LY = ly["disc_LY", ],
             row.names = NULL) %>% arrange(AIC)
}

results <- surv %>% group_split(arm) %>%
  setNames(sort(unique(surv$arm))) %>% lapply(fit_one_arm)
print(results)

# Selection: lowest AIC/BIC AND plausible smoothed hazard AND survival never exceeding the matched
# general population. For relative survival / flooring, add bhazard = d$bhaz to flexsurvreg() and use
# a cure model (flexsurvcure) when a genuine plateau is supported by mature follow-up.
# PSA: draw parameters from normboot.flexsurvreg(fit, B = 1000) and recompute disc_LY per draw.