Survival analysis bridge

Open In Colab

This tutorial shows how to turn a fitted parametric survival curve into inputs the engines accept: sampled event times for MicrosimModel.continuous and per-cycle death probabilities for MarkovModel. The reference curve is Weibull, S(t) = exp(-(t / scale) ** shape) with shape 1.2 and scale 6.0 years, run as a two-state alive-and-dead model at a 3% discount rate. Full script: examples/survival_bridge.py.

Defining the survival curve

The survival function and its inverse are all the two engines need: one to integrate, one to sample from.

import numpy as np
import pandas as pd

SHAPE, SCALE, DISCOUNT = 1.2, 6.0, 0.03

def survival(t, shape, scale):
    return np.exp(-((t / scale) ** shape))  # Weibull survival S(t)

def sample_times(rng, size, shape, scale):
    u = rng.random(size)  # invert S(t) = u for the event time
    return scale * (-np.log(u)) ** (1 / shape)

The discounted life expectancy is the integral of exp(-rate t) S(t), the value both engines target: 4.927 years.

from scipy.integrate import quad
discounted_le = quad(lambda t: np.exp(-DISCOUNT * t) * survival(t, SHAPE, SCALE), 0, np.inf)[0]
round(discounted_le, 5)  # discounted life expectancy, the value both engines target
4.92709

Sampling event times for an individual model

MicrosimModel.continuous races one death time per person; event_times returns the time to each state, with inf where a move cannot happen.

from heormodel.models import MicrosimModel
from heormodel.run import run_psa

STATES, ARM = ("alive", "dead"), "Standard care"

def event_times(params, intervention, state, attrs, rng):
    times = np.full((len(state), 2), np.inf)  # columns: to alive, to dead
    alive = state == 0
    times[alive, 1] = sample_times(rng, int(alive.sum()), params["shape"], params["scale"])
    return times

def reward_rates(params, intervention, state, attrs):
    alive = (state == 0).astype(float)
    return np.zeros(len(state)), alive  # no cost; one life-year per year alive

continuous = MicrosimModel.continuous(
    states=STATES, event_times=event_times, state_reward_rates=reward_rates,
    interventions=[ARM], horizon=80.0, n_individuals=200_000,
    discount_rate=DISCOUNT, effect="lifeyears")

The sampled cohort recovers the analytic discounted life expectancy within Monte Carlo error.

draws = pd.DataFrame({"shape": [SHAPE], "scale": [SCALE]}, index=pd.RangeIndex(1, name="iteration"))
run_psa(continuous, draws, seed=1, sequential=True).outcomes.summary().loc[ARM, "lifeyears"]
np.float64(4.930517038770994)

Building per-cycle transition probabilities for a cohort model

The other consumer is MarkovModel: the per-cycle probability of dying is 1 - S(k+1) / S(k), stacked into an age-varying transition array.

from heormodel.models import CohortSpec, MarkovModel

N_CYCLES = 60

def transitions_and_rewards(params, intervention):
    surv = survival(np.arange(N_CYCLES + 1), params["shape"], params["scale"])
    death = 1 - surv[1:] / surv[:-1]  # conditional death probability each cycle
    transition = np.zeros((N_CYCLES, 2, 2))
    transition[:, 0, 0], transition[:, 0, 1] = 1 - death, death
    transition[:, 1, 1] = 1.0  # dead is absorbing
    return CohortSpec(transition, np.zeros(2), np.array([1.0, 0.0]))

cohort = MarkovModel(
    states=STATES, interventions=[ARM], transitions_and_rewards=transitions_and_rewards,
    n_cycles=N_CYCLES, initial_state="alive", discount_rate=DISCOUNT,
    cycle_correction="half_cycle", effect="lifeyears")
cohort.evaluate(draws).summary().loc[ARM, "lifeyears"]
np.float64(4.945314843927852)

The cohort lands within the half-cycle correction error of the continuous value, so a model author gets the same answer whichever engine they use.

Fitting the curve and carrying its uncertainty

A real analysis fits the curve to data first. Here a Weibull is fit by maximum likelihood to a right-censored sample, returning the shape, scale, and their log-scale covariance.

from scipy.optimize import minimize

def fit_weibull(times, events):
    def negative_log_likelihood(theta):
        shape, scale = np.exp(theta)  # optimize on the log scale so both stay positive
        return -np.sum(events * (np.log(shape / scale) + (shape - 1) * np.log(times / scale))
                       - (times / scale) ** shape)
    fit = minimize(negative_log_likelihood, np.log([1.0, times.mean()]), method="BFGS")
    return np.exp(fit.x), fit.hess_inv  # (shape, scale), log-scale covariance

rng = np.random.default_rng(20260714)
event_time = sample_times(rng, 300, SHAPE, SCALE)  # a 300-patient trial
observed = np.minimum(event_time, 12.0)            # administrative censoring at 12 years
(shape_hat, scale_hat), cov_log = fit_weibull(observed, (event_time <= 12.0).astype(float))
round(shape_hat, 3), round(scale_hat, 3)
(np.float64(1.203), np.float64(5.987))

Drawing shape and scale from the fit’s asymptotic distribution puts the estimation uncertainty on the iteration index, so run_psa propagates it like any other parameter.

log_draws = rng.multivariate_normal(np.log([shape_hat, scale_hat]), cov_log, size=1_000)
fitted = pd.DataFrame(np.exp(log_draws), columns=["shape", "scale"],
                      index=pd.RangeIndex(1, 1_001, name="iteration"))
life_years = run_psa(continuous, fitted, seed=2).outcomes.effects_wide()[ARM]
round(life_years.mean(), 3), np.round(np.percentile(life_years, [2.5, 97.5]), 3)
(np.float64(4.919), array([4.503, 5.34 ]))

The 95% credible interval contains the analytic 4.927, and as the trial size grows the estimates converge to shape 1.2 and scale 6.0 and the interval shrinks toward that value.

Back to top