Survival analysis bridge (phase 1)

A validation exhibit: the same reference model reproduced with example-local functions, before heormodel.survival existed.

Open In Colab

This page is a validation exhibit, not the current way to bring a survival curve into heormodel. It reproduces the reference model below with example-local functions, as it was done before heormodel.survival existed. Use Survival curves instead, which covers the same model through the public module.

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 draws 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

For 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 result is 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 gives one parameter set per iteration, so run_psa propagates the estimation uncertainty 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% interval contains the analytic 4.927. As the trial size grows, the estimates converge to shape 1.2 and scale 6.0 and the interval narrows around that value.

Back to top