Multi-state models

Competing risks describe a single transition from an initial state to one of several end states. Multi-state models generalize this to processes where subjects move through several states over time, possibly passing through intermediate ones. The classic example is the illness-death model: a subject starts healthy, may become ill, and may die either directly or after becoming ill. Greenwood estimates the transition and occupancy probabilities for such models with the Aalen-Johansen estimator. This page shows how to prepare multi-state data and interpret the results.

flowchart LR
    H[Healthy] --> I[Ill]
    H --> D[Death]
    I --> D

Figure 1: The illness-death model. Subjects can move from health to illness, and reach death from either state.

Preparing multi-state data

A multi-state analysis needs one row per period a subject spends in a state. Each row records the interval (start, stop], the state occupied during it, and the state transitioned to at the end, or a censoring marker if none. A subject who becomes ill and then dies contributes two rows: one for the healthy period ending in the illness transition, and one for the ill period ending in death.

We build this from mgus2, treating progression as the intermediate “pcm” state and death as the absorbing state.

import greenwood as gw

# Load the mgus2 dataset
mg = gw.load_dataset("mgus2", backend="polars")

# Build multi-state intervals: each subject gets one or two rows
start, stop, state, event = [], [], [], []
for i in range(len(mg)):
    pt, ft = mg["ptime"][i], mg["futime"][i]
    progressed, died = mg["pstat"][i] == 1, mg["death"][i] == 1
    if progressed and pt < ft:
        start += [0, pt]
        stop += [pt, ft]
        state += ["mgus", "pcm"]
        event += ["pcm", "death" if died else None]
    else:
        start += [0]
        stop += [ft]
        state += ["mgus"]
        event += ["death" if died else ("pcm" if progressed else None)]

# Keep only positive-length intervals
rows = [(a, b, s, e) for a, b, s, e in zip(start, stop, state, event) if b > a]
start, stop, state, event = map(list, zip(*rows))

The loop appends to four parallel lists, so a single interval is spread across start[i], stop[i], state[i], and event[i] at the same index. After the loop we drop any zero-length intervals and unpack the surviving rows back into the four lists. Counting the entries tells us how many intervals we ended up with, which is more than the number of subjects because progressing subjects contribute two rows each.

# Total number of intervals (more than subjects due to progressors)
len(start)

Reading the four lists position by position is easier if we zip them back together and look at the first several intervals as (start, stop, state, event) tuples.

# Preview the first six intervals as (start, stop, state, event) tuples
list(zip(start, stop, state, event))[:6]

[(0, 30, 'mgus', 'death'),
 (0, 25, 'mgus', 'death'),
 (0, 46, 'mgus', 'death'),
 (0, 92, 'mgus', 'death'),
 (0, 8, 'mgus', 'death'),
 (0, 4, 'mgus', 'death')]

Each tuple is one interval (start, stop]: the subject occupies state from just after start up to and including stop, and at stop transitions to event. A None in the event position marks a censored interval, one that ends without a transition. A subject who only ever holds mgus produces a single row, while a subject who progresses produces two consecutive rows, the second of which has state pcm, exactly the illness-death pattern the model needs.

Estimating transition and occupancy probabilities

The MultiState estimator forms the Aalen-Johansen product of the per-time transition matrices and reports, at each time, the probability of occupying each state. This is the multi-state generalization of the survival curve.

# Rebuild the intervals from scratch for a self-contained example
mg = gw.load_dataset("mgus2", backend="polars")

start, stop, state, event = [], [], [], []
for i in range(len(mg)):
    pt, ft = mg["ptime"][i], mg["futime"][i]
    progressed, died = mg["pstat"][i] == 1, mg["death"][i] == 1
    if progressed and pt < ft:
        start += [0, pt]
        stop += [pt, ft]
        state += ["mgus", "pcm"]
        event += ["pcm", "death" if died else None]
    else:
        start += [0]
        stop += [ft]
        state += ["mgus"]
        event += ["death" if died else ("pcm" if progressed else None)]

# Keep only positive-length intervals
rows = [(a, b, s, e) for a, b, s, e in zip(start, stop, state, event) if b > a]
start, stop, state, event = map(list, zip(*rows))

# Fit the multi-state Aalen-Johansen estimator from lists
ms = gw.MultiState().fit(
    start=start, stop=stop, state=state, event=event, states=("mgus", "pcm", "death")
)

When the intervals already live in a data frame, name the columns and pass the frame as data= instead. This works for any backend.

import polars as pl

# Package the same intervals as a Polars DataFrame
intervals = pl.DataFrame({"tstart": start, "tstop": stop, "from": state, "to": event})

# Fit using column names and data= instead of raw lists
ms = gw.MultiState().fit(
    start="tstart",
    stop="tstop",
    state="from",
    event="to",
    states=("mgus", "pcm", "death"),
    data=intervals,
)

The fitted estimator holds an occupancy curve for each of the three states. In its data frame the columns after time give the occupancy probability of mgus, pcm, and death respectively, and the preview reports the full number of time points.

# Occupancy probabilities over time for each state
ms.to_frame(format="polars")

	time f64	mgus f64	pcm f64	death f64
PolarsRows276Columns4
0	1	0.969653179191	0	0.0303468208092
1	2	0.947960714332	0.00144616432392	0.0505931213442
2	3	0.937114481902	0.00144616432392	0.0614393537736
3	4	0.924822085149	0.00216924648589	0.073008668365
4	5	0.916868181368	0.00289232864785	0.0802394899846
5	6	0.909637359748	0.00361541080981	0.0867472294423
6	7	0.905298866776	0.00361541080981	0.091085722414
7	8	0.895175716509	0.00433849297177	0.10048579052
273	376	0.054500072561	0.0272500362805	0.918249891158
274	394	0.054500072561	0.0272500362805	0.918249891158
275	424	0	0.0272500362805	0.972749963719

Each column after time is the probability that a subject is in that state at that time. The probabilities sum to one across states, because every subject is always in exactly one state. Early on, almost everyone is in the initial mgus state; over time, mass flows into pcm and, predominantly, death.

You can read the occupancy probabilities at specific times with predict.

# Read occupancy probabilities at specific time points
ms.predict(times=[60, 120, 240], format="polars")

	time f64	mgus f64	pcm f64	death f64
PolarsRows3Columns4
0	60	0.645529276758	0.0160062741857	0.338464449057
1	120	0.404460127907	0.0120267339061	0.583513138187
2	240	0.176158307922	0.0114890821233	0.812352609955

The full transition-probability matrices, giving the probability of moving from each state to each other state over the interval from time zero, are available as the transition_ attribute for users who need them.

Occupancy versus transition probabilities

Occupancy probabilities answer “what is the chance a subject is in state k at time t?”. The transition matrix answers the more detailed “given a subject in state j at the start, what is the chance they are in state k at time t?”. The occupancy probabilities are the initial-state distribution multiplied through the transition matrix.

Competing risks is a special case

A competing-risks model is a multi-state model with one initial state and several absorbing states and no intermediate transitions. If that is your situation, use AalenJohansen directly, which also gives standard errors for the cumulative incidence.

Next steps

You can now estimate occupancy probabilities for general multi-state processes.

Competing risks covers the important special case with standard errors and the Fine-Gray regression model.
Survival data explains the counting-process response and the categorical event (first category = censoring) that underlie multi-state data.
The Quick start links every topic together if you want a map of the whole package.