Cure Models (Mixture and Non-Mixture)

What cure models are and why they matter

Standard parametric survival models — Weibull, log-normal, Gompertz, and their relatives — assume that every patient will eventually experience the event if followed long enough. This assumption is wrong in a growing class of clinical settings: checkpoint-inhibitor immuno-oncology, certain haematological malignancies, childhood acute lymphoblastic leukaemia, and some chronic-disease applications where a subgroup achieves durable biological remission or immunological control and will never experience the event regardless of follow-up length.

When this two-population structure exists, forcing a standard model produces a misspecified extrapolation — a survival curve that asymptotes to zero and systematically underestimates mean survival in the treated arm.

Cure models, also called bounded survival models or long-term survivor models, directly parameterize this structure by estimating both the size of the cured fraction and the event-time distribution among those who remain at risk.

The cure phenomenon: when the Kaplan-Meier plateau tells a story

The canonical diagnostic is a sustained plateau in the Kaplan-Meier (KM) estimator at long follow-up. If, after sufficient events have accumulated, the KM curve stabilizes at a level materially above zero and remains stable rather than continuing to decline, this is evidence that a non-negligible fraction of the population may have escaped the event permanently.

In immuno-oncology this pattern was first observed consistently in melanoma and non-small-cell lung cancer trials evaluating PD-1/PD-L1 checkpoint inhibitors: five-year OS rates of 20–30% in populations where median OS with prior standard of care was under 12 months.

The plateau is the fingerprint of the durable-responder subgroup.

Three conditions strengthen the case for a genuine cure plateau rather than an artifact:

1. the plateau is stable across multiple consecutive data cuts;
2. surviving patients are followed well beyond the median event time; and
3. background all-cause mortality from age and comorbidity begins to visibly erode the plateau at very long follow-up — consistent with relative-survival cure models discussed below.

When any of these conditions fails, treat the apparent plateau with caution. The existing Kaplan-Meier entry covers the KM estimator mechanics; this entry focuses on what to do once a plateau is diagnosed.

The mixture cure model: S(t) = pi + (1 − pi) × S_u(t)

The mixture cure model, introduced by Boag (1949) and formalized by Farewell (1982), decomposes overall population survival into two components. Let pi denote the cure fraction — the proportion of patients who will never experience the event — and let S_u(t) denote the conditional survival function for the uncured (susceptible) subgroup. Overall survival at time t is:

S(t) = pi + (1 − pi) × S_u(t)

As time approaches infinity, S_u(t) approaches 0 (every uncured patient eventually has the event), so S(t) approaches pi. The long-run survival plateau equals the cure fraction exactly.

The model has two distinct parameter sets with two distinct interpretations. The cure probability pi — and any covariates predicting it — is typically estimated through a logistic regression on baseline patient characteristics: logit(pi_i) = gamma_0 + gamma_1 x_i. The latency survival S_u(t | z_i) is modeled separately with its own covariate vector z_i and a chosen parametric family (Weibull is the default; log-normal, log-logistic, or Royston-Parmar splines are alternatives).

This two-component structure means a treatment can increase pi (cure more patients) without changing the speed of uncured progression, change the latency without altering the cure fraction, or do both simultaneously — a clinical decomposition no single-hazard survival model can provide.

The log-likelihood for a mixture cure model with exponential latency is:

For observed events: log[(1 − pi) × f_u(t)] For censored observations: log[pi + (1 − pi) × S_u(t)]

The non-mixture (promotion-time) cure model

An alternative parameterization replaces the mixture structure with a competing-cause framework. The promotion-time model (Yakovlev and Tsodikov 1996; Chen, Ibrahim, and Sinha 1999) supposes each individual has a random number N of competing promotion events (e.g., metastatic foci) drawn from a Poisson distribution; the event occurs at the minimum of N independent promotion times. When N = 0, the individual is cured. The survival function is S(t) = exp(−theta × F(t)), where F(t) is the CDF of a promotion-time distribution and theta is the Poisson mean.

The cure fraction is exp(−theta).

The non-mixture model is biologically motivated for cancer applications where the latent number of metastatic foci drives the event time; it also admits a proportional-hazards reparameterization on the promotion-time scale, making it tractable in standard software. For practical HTA applications with limited follow-up, mixture and non-mixture models are often empirically indistinguishable — model selection criteria and biological justification should guide the choice. In most NICE submissions, the mixture cure model is the more common framing.

Identifiability and the credibility threshold

Identifiability is the central methodological challenge of cure models. A cure fraction can only be credibly estimated if observed follow-up substantially exceeds the time when uncured events effectively cease — that is, if S_u(t) has approached zero in the observed window. If many patients remain at risk in the long tail, the model cannot reliably distinguish a true cure fraction from a very slow uncured survival function.

The practical implication for HTA submissions is significant: claiming a cure fraction of 20% from a two-year trial with a barely flattened KM tail is not credible; claiming 30% from a five-year trial with a stable plateau, zero late events, and background-mortality anchoring is considerably more defensible.

Formal identifiability tools include the Maller–Zhou (1992) sufficient-follow-up test, which assesses whether the largest observed event time is sufficiently close to the largest follow-up time to rule out an artifact plateau. Analysts should also inspect the smoothed hazard function: in the cured subgroup, the hazard should converge toward the background (population life-table) hazard rather than toward zero. The cure fraction estimate should remain stable across consecutive data cuts; a cure fraction that rises substantially from one interim data cut to the next is a sign of insufficient follow-up.

Relative-survival cure models and background-mortality anchoring

Standard mixture cure models estimated from all-cause survival do not account for background mortality. At long follow-up, even "cured" patients begin to die from unrelated causes, causing the observed KM plateau to gradually erode. Relative-survival cure models address this by modeling excess mortality relative to expected population survival, allowing the cure fraction to be interpreted as the proportion of patients who experience no excess disease-specific mortality — their mortality profile matches the general population's.

This is the approach recommended when follow-up exceeds five years or when older, comorbid cohorts are analyzed where background mortality is non-trivial.

The related concept on relative net survival covers the relative-survival framework and population mortality table linkage in detail.

HTA stakes: the cure fraction dominates lifetime QALY extrapolations

In health technology assessment submissions where the comparator shows a declining KM curve and the treatment arm shows a plateau, the cure fraction is frequently the single most influential assumption in the economic model. NICE Technical Support Document 21 (2021) on flexible methods for survival extrapolation specifically addresses cure models: it recommends using them only when genuine clinical evidence for long-term survivors exists, requires sensitivity analyses across alternative cure fractions, and flags the cure assumption as a focus area for ERG/EAG technical scrutiny.

The stakes are high: the difference between a 20% and a 30% assumed cure fraction in a lifetime cost-effectiveness model can shift the ICER by tens of thousands of pounds per QALY — enough to change the decision from approved to rejected.

Analysts must pre-specify the cure-model scenario in the analysis plan, report it alongside the full standard parametric candidate model set, and provide both clinical and statistical justification for any assumed cure fraction.

Covariates on cure probability versus latency separately

One of the most clinically informative features of the mixture cure model is the ability to specify separate covariates for the cure fraction and the latency distribution. This means an analyst can ask: does biomarker B predict who is cured (enters the logistic regression on pi), and does dose intensity predict how quickly uncured patients progress (enters the latency hazard)?

Standard single-hazard survival regression forces both mechanisms into a single coefficient set and cannot answer either question. Implementations in R (flexsurvcure, smcure) and the custom PROC NLMIXED approach in SAS both support separate covariate vectors.

RWE-specific caution: apparent plateaus as data-cutoff artifacts

In real-world evidence settings — claims databases, EHR cohorts, and registry studies — observed plateaus in rwPFS or rwOS frequently reflect data structure rather than biology. Three artifacts create false plateaus:

1. administrative censoring pile-up, where a large fraction of patients are censored at a common database end date, producing an abrupt flat segment;
2. outcome ascertainment gaps, where deaths are under-captured in claims data so that deceased patients appear as long-term survivors; and
3. short follow-up with sparse late events, where the KM estimator lacks events to show a continuing decline.

Before fitting a cure model to RWE data, analysts should:

1. plot the censoring distribution separately from the event distribution to detect pile-up;
2. restrict to patients enrolled early enough to have genuinely long follow-up;
3. verify death ascertainment against a linked mortality source (NDI, Social Security Death Index, or vital statistics) where available; and
4. confirm the plateau does not coincide with the administrative data-cut date.

Pros, cons, and trade-offs

Pros of mixture cure models: directly parameterize the two-population structure present in immuno-oncology durable response; produce a long-run survival plateau that standard models cannot represent; decompose treatment effects into cure-fraction gain versus latency improvement; generate clinically interpretable quantities for HTA narrative; align with NICE TSD 21 guidance for innovative immuno-oncology submissions; latency component is flexible (Weibull, spline, log-normal).

Cons: identifiability requires mature follow-up and a genuinely stable plateau — applying cure models to early or immature data inflates the estimated cure fraction dramatically and produces unreliable extrapolations; the two-component structure doubles the number of parameters and inflates uncertainty in probabilistic sensitivity analysis; cure-fraction estimates are sensitive to the choice of latency distribution in the sparse tail; regulatory and HTA reviewers subject cure assumptions to intensive scrutiny.

Key trade-offs versus standard parametric models: standard Weibull or log-normal models are simpler, identifiable under any follow-up length, and easier to use in PSA — but they cannot represent a bounded long-run survival and will systematically underestimate mean survival in populations with genuine durable responders. The cure model gains biological fidelity at the cost of identifiability requirements and two-component complexity. Report both sets of models; do not present the cure model as the single primary analysis unless identifiability is well established.

When NOT to use — and when cure models are actively misleading

Do not use cure models when:

1. follow-up is immature (less than three to five years for most oncology indications, or the KM curve has not clearly plateaued);
2. the apparent plateau coincides with or is explained by data-cutoff censoring pile-up;
3. the Maller-Zhou identifiability test indicates insufficient follow-up to distinguish a true cure fraction from a very slow uncured tail;
4. the indication has no biological basis for permanent event prevention — most progressive neurodegenerative diseases, most cardiovascular endpoints in uncontrolled disease;
5. the cure fraction estimate changes substantially between consecutive data cuts, indicating the model is tracking a moving artifact rather than a stable biological phenomenon.

A cure model with wide confidence intervals on pi, or whose pi estimate requires a logistic regression with only a handful of long-term survivors, should not be the primary analysis model — report it as a sensitivity scenario alongside the standard parametric candidate set. Actively misleading use: claiming a cure fraction in RWE data whose apparent plateau coincides exactly with the administrative database end date is a common error that can grossly overestimate long-term efficacy.

Interpreting the output

In the worked example, a mixture cure model fit to a synthetic 50-patient immunotherapy trial estimates pi = 0.30 (30% cure fraction) and an exponential latency with S_u(24) = 0.50 (half of all uncured patients remain event-free at 24 months).

(1) Formal interpretation. The overall survival at 24 months is S(24) = 0.30 + 0.70 * 0.50 = 0.30 + 0.35 = 0.65. Sixty-five percent of the population is expected to be event-free at 24 months: 30 percentage points because they are in the permanently cured subgroup, and 35 percentage points because they are in the uncured subgroup (70% of total) but have not yet experienced the event (S_u(24) = 0.50 means 50% of uncured individuals survive to 24 months).

As time approaches infinity, S(t) approaches pi = 0.30 — the long-run plateau.

The cure fraction estimate and its 95% CI drive the long-run QALY extrapolation; the latency parameter drives near-term QALY accumulation. The two components are separately interpreted and separately interrogated by HTA technical reviewers.

(2) Practical interpretation. Of 100 patients starting treatment, the model estimates 30 will never progress or die from this cancer regardless of follow-up. The remaining 70 remain at risk, and at 24 months approximately half of them (35 patients) are still event-free, giving 65 event-free patients in total.

For a lifetime cost-effectiveness model, the 30% cure fraction means the treated arm shows a non-zero long-run survival plateau — a qualitatively different curve shape from any standard parametric model asymptoting to zero. If the cure fraction is raised from 30% to 40% in a deterministic sensitivity analysis, the entire tail of the survival curve lifts by 10 percentage points, typically adding a material number of life-years gained and potentially shifting the ICER by tens of thousands of pounds per QALY.

Decision diagram

flowchart LR
  subgraph Population["All patients at time zero"]
    direction TB
    C["Cured subgroup<br/>proportion = pi<br/>Never experience the event<br/>S_cured(t) = 1 for all t"]
    U["Uncured subgroup<br/>proportion = 1 - pi<br/>Will eventually experience event<br/>S_u(t) → 0 as t → ∞"]
  end
  Population --> S["Overall survival<br/>S(t) = pi + (1 − pi) × S_u(t)<br/>Long-run plateau = pi"]
  S --> HTA["HTA extrapolation<br/>Lifetime horizon QALY<br/>dominated by pi"]
  S --> ID["Identifiability check<br/>KM plateau stable across<br/>data cuts? Maller-Zhou test?"]

Worked example

Scenario

A biostatistician is analyzing overall survival from a 50-patient single-arm trial of a PD-1 checkpoint inhibitor in metastatic melanoma. The Kaplan-Meier curve shows a stable plateau at approximately 30% beginning around month 30, with no events after that point across two consecutive data cuts. A mixture cure model is fitted using maximum likelihood: the model estimates a cure fraction pi = 0.30 and an exponential latency distribution for uncured patients with rate chosen so that the median uncured survival is 24 months, yielding S_u(24) = 0.50. The task is to compute the overall survival probability at 24 months and confirm the long-run plateau.

Dataset

Representative survival data for 10 of the 50 trial patients (synthetic). event_observed = 1 means death was recorded; event_observed = 0 means the patient was alive at last contact (last data cut). Patients still event-free beyond 30 months are the long-term survivors whose sustained follow-up provides the evidence for the cure fraction.

Steps

Result

S(24) = pi + (1 - pi) * S_u(24) = 0.30 + 0.70 * 0.50 = 0.30 + 0.35 = 0.65. At 24 months, 65% of patients are expected to be event-free: 30% because they are permanently cured and 35% because they are in the uncured group but have not yet experienced the event. The long-run plateau is pi = 0.30.

Trade-offs

Runnable example

import numpy as np
from scipy.optimize import minimize

def mixture_cure_loglik(params, t, event):
    """Negative log-likelihood for mixture cure model with exponential latency.

    params[0] = logit(pi)  -- cure fraction on log-odds scale (unconstrained)
    params[1] = log(lam)   -- exponential rate on log scale (ensures lam > 0)
    t      : array of observed times (event or censoring)
    event  : array of 0/1 indicators (1 = event observed, 0 = censored)
    """
    logit_pi, log_lam = params
    pi  = 1.0 / (1.0 + np.exp(-logit_pi))   # cure fraction in (0, 1)
    lam = np.exp(log_lam)                     # positive exponential rate
    S_u = np.exp(-lam * t)                    # S_u(t) = exp(-lambda * t)
    f_u = lam * S_u                           # density: f_u(t) = lambda * exp(-lambda * t)
    # Log-contributions per observation
    ll_event  = np.log((1.0 - pi) * f_u + 1e-15)  # observed event
    ll_censor = np.log(pi + (1.0 - pi) * S_u + 1e-15)  # censored
    ll = np.sum(event * ll_event + (1.0 - event) * ll_censor)
    return -ll  # scipy.minimize minimises, so return negative LL

# ── Synthetic data: 35 uncured (events) + 15 cured (censored at 60 months) ──
np.random.seed(42)
n_cured, n_uncured = 15, 35
t_uncured = np.random.exponential(scale=24, size=n_uncured)  # median ~17 months
t_cured   = np.full(n_cured, 60.0)                           # censored at 60 months
t_all     = np.concatenate([t_uncured, t_cured])
event_all = np.concatenate([np.ones(n_uncured), np.zeros(n_cured)])

# ── Fit via Nelder-Mead ──
init   = [0.0, np.log(1.0 / 24.0)]     # pi = 0.5, rate ≈ 1/24
result = minimize(mixture_cure_loglik, init,
                  args=(t_all, event_all), method="Nelder-Mead",
                  options={"xatol": 1e-8, "fatol": 1e-8, "maxiter": 10000})

logit_pi_hat, log_lam_hat = result.x
pi_hat  = 1.0 / (1.0 + np.exp(-logit_pi_hat))
lam_hat = np.exp(log_lam_hat)

print(f"Converged: {result.success}  |  Iterations: {result.nit}")
print(f"Estimated cure fraction pi: {pi_hat:.3f}")
print(f"Estimated exponential rate lambda: {lam_hat:.4f}")
print(f"Estimated median uncured survival: {np.log(2)/lam_hat:.1f} months")

# ── Compute S(t) at 24 months using the worked-example formula ──
t_eval = 24.0
S_u_24 = np.exp(-lam_hat * t_eval)
S_24   = pi_hat + (1.0 - pi_hat) * S_u_24
print(f"\nS_u(24 months) = {S_u_24:.3f}")
print(f"S(24 months)   = {pi_hat:.3f} + (1 - {pi_hat:.3f}) * {S_u_24:.3f} = {S_24:.3f}")
print(f"Long-run plateau (t -> infinity): S -> pi = {pi_hat:.3f}")

# ── Survival curve at selected time points ──
t_grid = np.array([6, 12, 18, 24, 36, 48, 60])
for t in t_grid:
    S_t = pi_hat + (1.0 - pi_hat) * np.exp(-lam_hat * t)
    print(f"  S({t:2d}) = {S_t:.3f}")

# install.packages(c("flexsurvcure", "smcure", "survival"))
library(flexsurvcure)
library(survival)

# ── Synthetic dataset mirroring the worked example ──
set.seed(42)
n_cured   <- 15L; n_uncured <- 35L
t_data    <- c(rexp(n_uncured, rate = 1/24),   # uncured: exponential, mean = 24 months
               rep(60, n_cured))                # cured:   censored at 60 months
event_data <- c(rep(1L, n_uncured), rep(0L, n_cured))
df <- data.frame(time = t_data, status = event_data)

# ── Mixture cure model: exponential latency, no covariates ──
# mixture = TRUE  → S(t) = pi + (1-pi)*S_u(t)  [mixture parameterisation]
# dist    = "exp" → exponential latency for uncured patients
fit_mcm <- flexsurvcure(Surv(time, status) ~ 1, data = df,
                         dist = "exp", mixture = TRUE)
print(fit_mcm)

# ── Extract cure fraction and latency rate ──
pi_hat  <- plogis(coef(fit_mcm)["theta"])    # theta is logit(pi) in flexsurvcure
lam_hat <- exp(coef(fit_mcm)["rate"])         # log-scale rate parameter

cat(sprintf("Estimated cure fraction pi : %.3f\n", pi_hat))
cat(sprintf("Latency rate lambda         : %.4f\n", lam_hat))
cat(sprintf("Median uncured survival     : %.1f months\n", log(2) / lam_hat))

# ── Compute S(24) using the worked-example formula ──
S_u_24 <- exp(-lam_hat * 24)
S_24   <- pi_hat + (1 - pi_hat) * S_u_24
cat(sprintf("\nS_u(24 months) = %.3f\n",  S_u_24))
cat(sprintf("S(24 months)   = %.3f + %.3f * %.3f = %.3f\n",
            pi_hat, (1 - pi_hat), S_u_24, S_24))
cat(sprintf("Long-run plateau S(inf) = pi = %.3f\n", pi_hat))

# ── Survival predictions at a grid of time points ──
t_seq <- seq(0, 72, by = 6)
pred  <- summary(fit_mcm, t = t_seq, type = "survival")[[1]]
print(pred[, c("time", "est", "lcl", "ucl")])

# ── Non-mixture (promotion-time) cure model comparison ──
fit_nmcm <- flexsurvcure(Surv(time, status) ~ 1, data = df,
                          dist = "exp", mixture = FALSE)
cat("\nNon-mixture model AIC:", AIC(fit_nmcm),
    "  Mixture model AIC:", AIC(fit_mcm), "\n")

# ── smcure: separate covariate formulas for cure fraction and latency ──
# library(smcure)
# Suppose we have a binary treatment indicator 'trt' and a continuous biomarker 'bm'
# smcure(Surv(time, status) ~ trt,            # covariates on latency (PH or AFT)
#        cureform = ~ trt + bm,               # separate logistic on pi
#        data = df, model = "ph", nboot = 200)
# The cureform argument controls which variables enter the logistic cure-fraction model;
# the main formula controls the latency hazard. Both can overlap with different subsets.

/* ── Create synthetic survival dataset ── */
data work.trial;
  call streaminit(42);
  do i = 1 to 50;
    if i <= 35 then do;   /* Uncured: Exponential(mean=24) */
      time   = rand('exponential') * 24;
      status = 1;
    end;
    else do;              /* Cured: censored at 60 months  */
      time   = 60;
      status = 0;
    end;
    drop i;
    output;
  end;
run;

/* ── Mixture cure model via PROC NLMIXED ──────────────────────────────
   Log-likelihood per observation:
     Event   (status=1): log[(1-pi) * lam * exp(-lam*t)]
     Censored(status=0): log[pi + (1-pi) * exp(-lam*t)]
   Parameters estimated on unconstrained scales:
     logit_pi: logit of cure fraction pi
     log_lam : log of exponential rate lambda                           */
proc nlmixed data=work.trial;
  /* Constrained-scale quantities */
  pi  = exp(logit_pi)  / (1 + exp(logit_pi));    /* cure fraction in (0,1)   */
  lam = exp(log_lam);                              /* exponential rate > 0     */
  S_u = exp(-lam * time);                          /* S_u(t) = exp(-lambda*t)  */
  f_u = lam * S_u;                                 /* f_u(t) = lambda*exp(...)  */
  S   = pi + (1 - pi) * S_u;                      /* overall mixture survival  */

  /* Log-likelihood contribution (GENERAL distribution in NLMIXED) */
  if status = 1 then ll = log((1 - pi) * f_u + 1e-12);   /* event    */
  else               ll = log(S + 1e-12);                  /* censored */
  model time ~ general(ll);

  /* Starting values (logit(0.3)=-0.847, log(1/24)=-3.178) */
  parms logit_pi = -0.847  log_lam = -3.178;

  /* Back-transformed estimates of clinical interest */
  estimate 'Cure fraction pi'
    exp(logit_pi) / (1 + exp(logit_pi));
  estimate 'Latency rate lambda'
    exp(log_lam);
  estimate 'Median uncured survival (months)'
    log(2) / exp(log_lam);
  estimate 'S_u(24 months)'
    exp(-exp(log_lam) * 24);
  estimate 'Overall S(24 months)  [= pi + (1-pi)*S_u(24)]'
    exp(logit_pi) / (1 + exp(logit_pi))
    + (1 - exp(logit_pi) / (1 + exp(logit_pi))) * exp(-exp(log_lam) * 24);
run;

/* ── NOTES ─────────────────────────────────────────────────────────────
   To add covariates on the cure fraction, replace logit_pi with a linear
   predictor: logit_pi = gamma0 + gamma1 * trt + ...; add the corresponding
   PARMS entries. Similarly, covariates on the latency rate: replace log_lam
   with a linear predictor.
   PROC NLMIXED uses quasi-Newton optimisation with Hessian-based standard
   errors; for complex models or sparse data, add TECHNIQUE=NRRIDG and
   verify convergence with the HESSIAN option.
   For the non-mixture (promotion-time) model, replace S with
   S = exp(-theta * (1 - S_u)) and the event density with the corresponding
   derivative; see Lambert (2007) Stata Journal for the non-mixture log-lik. */

Citations

1. [1]Othus M, Barlogie B, LeBlanc ML, Crowley JJ. Cure models as a useful statistical tool for analyzing survival. Clinical Cancer Research. 2012;18(14):3731-3736.
2. [2]Lambert PC. Modeling of the cure fraction in survival studies. The Stata Journal. 2007;7(3):351-375.

1. [3]Amdahl J. flexsurvcure: Flexible Parametric Cure Models. CRAN. 2017.