Fitting an Illness-Death Model from Clinical Follow-Up Data

From Clinical Follow-Up to an Illness-Death Model

Clinical follow-up often records whether a patient develops an illness, dies, or is censored before either event. These records contain a history through states, not merely a single time to an endpoint. An illness-death model makes that history explicit: patients begin in the healthy state, may develop illness, and may die either before or after illness.

The usual illness-death model has three states: state 0 is alive without the illness, state 1 is alive after illness, and state 2 is dead. State 2 is absorbing. The permitted transitions are \(0\to1\), \(0\to2\), and \(1\to2\). A patient cannot make a \(0\to1\) transition after death, and a death after illness is a \(1\to2\) transition, not a second \(0\to2\) event.

Use a common time origin, such as randomization or cohort entry, and retain the actual event and censoring times. If illness occurs at time \(t_I\), a patient who remains under observation after illness contributes to the \(1\to2\) risk set starting at \(t_I\). In the clock-forward formulation used here, this interval is recorded on the time-since-origin scale, from \(t_I\) to death or censoring. A clock-reset or semi-Markov analysis uses a different time scale and requires its own risk-set construction.

Constructing Transition Records

A convenient counting-process representation has one record per patient per transition-specific risk interval. Each record contains a transition label, start and stop times, an event indicator, and the relevant covariates. The event indicator is 1 only if the event for that transition occurred at the stop time. Exiting a transition’s risk set because another transition occurred is recorded as a non-event for that transition.

In software, the same logic can be used to expand one patient-level record into transition-specific records, as in the earlier material on multi-state data expansion. Before fitting anything, check that intervals have valid start and stop times, that no patient is at risk for an outgoing transition after leaving its source state, and that each observed event agrees with the patient’s recorded history. Store baseline covariates consistently across a patient’s rows; time-varying covariates require values appropriate to each risk interval.

Transition-Specific Cox Models

For transition \(g\to h\), a Cox model specifies the transition intensity at time \(t\), conditional on being in state \(g\) just before \(t\). With baseline covariate vector \(x\), write

Each transition has its own baseline hazard \(\lambda_{0,gh}(t)\). Allowing a distinct coefficient vector \(\beta_{gh}\) means, for example, that a treatment can have different associations with illness, death before illness, and death after illness. A separate Cox fit for each transition is a direct implementation: subset to the transition’s risk intervals, use its event indicator, and fit a Cox model with start-stop time. Alternatively, a stratified fit can use transition-specific baseline hazards and transition-specific covariate terms.

The estimates above are illustrative model output, not a calculation from the four-patient table. That table is too small to support meaningful fits. In an actual analysis, check proportional hazards for each transition, including whether the covariate effect changes over time. Check event counts and risk-set sizes separately: a substantial cohort may still have few observed \(1\to2\) events. Also investigate convergence, influential observations, tied event times, and whether event dates and state histories are internally consistent.

Assumptions for Fitting and Predicting

The Cox fits require correctly identified event times, transition labels, and risk sets. Conditional on modeled covariates, censoring should be independent of the transition process in the relevant risk sets; informative loss to follow-up can distort both estimated hazards and state probabilities. The proportional-hazards form requires the covariate hazard ratio for each transition to be constant over the modeled time scale unless time-varying effects are explicitly included. Covariates used for a clinical profile must be measured and coded on the same scale as in model fitting.

A further assumption concerns the state-probability calculation. The standard forward calculation from transition intensities treats the process as Markov conditional on the modeled covariates and time scale: given the current state and covariates, the next transition intensity does not additionally depend on the earlier path. In an illness-death setting this can be a substantive assumption. For example, if death after illness depends on time since illness beyond what is represented by current state and covariates, a clock-forward Markov model may not adequately describe the process. The semi-Markov material addresses clock-reset risk sets; do not silently substitute one time scale for another.

For covariate profile \(x\), the fitted transition intensities form a matrix \(A(t\mid x)\). Its off-diagonal entries are the relevant intensities, and each diagonal entry is the negative sum of the other entries in its row. Thus the row for state 0 has entries for \(0\to1\) and \(0\to2\), while the row for state 1 has an entry for \(1\to2\); state 2 has no outgoing intensity. Starting in state 0, the state-probability row vector is obtained by advancing the initial vector through these transition intensities. For fitted Cox models, the transition-specific baseline cumulative hazards and profile-specific relative hazards supply the increments used in the product-integral calculation.

Derivation and a State-Probability Calculation

A short derivation shows why the illness probability depends on both illness and death hazards. Suppose the three transition intensities are constant, with rates \(a=\lambda_{01}\), \(b=\lambda_{02}\), and \(c=\lambda_{12}\), and patients start in state 0. Let \(p_j(t)\) be the probability of occupying state \(j\) at time \(t\). Over a short interval, state 0 loses probability through either outgoing transition, so \(p_0'(t)=-(a+b)p_0(t)\). State 1 gains probability from illness and loses it through death, so \(p_1'(t)=a p_0(t)-c p_1(t)\), with \(p_0(0)=1\) and \(p_1(0)=0\). Solving the first equation gives \(p_0(t)=e^{-(a+b)t}\). Substituting this into the second equation and integrating gives, when \(a+b\ne c\),

The integral accounts for every possible illness time \(s\): the patient must remain in state 0 until \(s\), enter state 1 at that time, and then remain alive in state 1 through \(t\). The death-state probability is the remaining probability, \(p_2(t)=1-p_0(t)-p_1(t)\). This constant-rate calculation is a transparent illustration; Cox-based predictions generally use estimated baseline hazard increments and a product integral rather than assuming constant rates.

With Cox fits, prediction for a profile \(x\) follows the same state-flow logic, but each transition hazard is multiplied by \(\exp(\hat\beta_{gh}^{\mathsf T}x)\). The fitted cumulative hazards are combined over event-time increments to obtain the transition-probability matrix from time 0 to the prediction horizon. The resulting row corresponding to the initial state gives the predicted probabilities of occupying each state. State probabilities are mutually exclusive and exhaustive, so their sum should be one, subject only to numerical rounding.

Diagnostics, Alternatives, and Limitations

• Audit the histories. A death after recorded illness belongs to \(1\to2\); a death before illness belongs to \(0\to2\). Check that no interval begins after its patient has left the source state.
• Check transition-specific model fit. Assess proportional hazards and inspect event counts for each transition. A satisfactory diagnostic for \(0\to1\) does not establish adequate fit for \(1\to2\).
• Check the time scale and Markov assumption. If post-illness risk depends on time since illness, a clock-forward model may miss that history dependence. Consider whether a clock-reset formulation is substantively appropriate.
• Check predictions. Probabilities must lie between 0 and 1 and sum to 1 at each horizon. Compare predictions with observed state occupancy where appropriate; disagreement can reveal model misspecification, but agreement alone does not prove the model is correct.

Transition-specific Cox regression is flexible about each baseline hazard and gives covariate effects without choosing a parametric hazard shape. Aalen-Johansen estimation provides nonparametric state-occupancy estimates for the observed cohort, but does not by itself give the same covariate-specific regression model. Parametric transition models impose more structure on the hazard shapes and can be useful for extrapolation, but their predictions can be sensitive to that structure. No approach repairs incorrect event histories or informative censoring, and a covariate-specific prediction remains conditional on the fitted model and its assumptions.

Common Mistake

A frequent error is to treat the fitted \(0\to1\) hazard ratio as the relative risk of being ill at a fixed time. The hazard ratio compares instantaneous illness transition intensities among patients still in state 0; it does not account by itself for deaths before illness, deaths after illness, or the timing of transitions. State occupancy requires combining all applicable transition hazards. Another error is to leave a patient in the \(0\to2\) risk set after illness: that patient is no longer in state 0 and cannot experience a \(0\to2\) transition thereafter.