Fitting Parametric Transition Models with flexsurv

Why Fit Parametric Transition Models?

A multi-state prediction depends on the transition hazards that move individuals between states. The Aalen–Johansen estimator and Cox-based approaches provide useful nonparametric or semiparametric predictions; parametric models instead specify a survival-time distribution for each transition. This gives smooth hazards and can support predictions at times where data are sparse, but it makes the distributional assumptions consequential.

We use a simple illness-death process: state 1 is illness-free, state 2 is illness, and state 3 is death. Transitions \(1\to2\) and \(1\to3\) can occur from state 1; transition \(2\to3\) can occur from state 2. State 3 is absorbing. The focus here is fitting separate Weibull or log-normal models to these transition-specific risk intervals and combining their hazards under a clock-forward Markov model.

For each transition, construct records describing time at risk, the transition event indicator, and any covariates. For a transition out of state 1, a different first transition ends that transition’s risk interval and is treated as a censoring event in that transition-specific cause-specific hazard model. For \(2\to3\), a person enters the risk set on entering state 2. In the clock-forward model used here, the time scale remains time since the common origin; the interval therefore starts at the time state 2 was entered. Correct entry and exit times are essential.

Fitting Weibull and Log-Normal Models

The R package flexsurv fits parametric survival distributions by maximum likelihood. In flexsurvreg, a Weibull model can be specified with dist="weibull" and a log-normal model with dist="lnorm". A counting-process survival response, Surv(Tstart,Tstop,status), represents an interval beginning at Tstart and ending at Tstop; status is 1 if the transition of interest occurred at its end and 0 otherwise. For a \(2\to3\) model on the clock-forward time scale, this form represents delayed entry at state-2 entry. Check the installed package’s documentation and confirm that its survival-response handling matches the data structure being supplied.

For example, with a multi-state long data set called ms, where trans identifies the transition, status marks the transition of interest, and Tstart and Tstop are its risk-interval endpoints, fits can be requested with flexsurvreg(Surv(Tstart,Tstop,status) ~ 1, data=subset(ms, trans=="1-2"), dist="weibull") and the same formula with dist="lnorm". Repeat for each transition. Add covariates only when the model is intended to represent their effects on that transition.

The separate cause-specific likelihood contribution for an observed interval ending at \(t_i\), beginning at \(a_i\), with event indicator \(\delta_i\), is proportional to the hazard at the event time if the event occurred, multiplied by survival over the interval. For a specified transition hazard \(h(t)\), this contribution is

The full transition-specific likelihood is the product over its risk intervals. Other transitions terminate an interval for this cause-specific model, but they are not treated as evidence that the person remained in the original state indefinitely. This likelihood-based fitting relies on correct transition coding, appropriate risk intervals, independent subjects, and censoring that is non-informative for the transition hazard conditional on modeled covariates and the observed risk history. It also assumes the chosen parametric family adequately represents the transition hazard.

The Weibull parameterization used here has shape \(k>0\) and scale \(\lambda>0\), with survival \(S(t)=\exp\{-(t/\lambda)^k\}\). Its hazard is \(h(t)=(k/\lambda)(t/\lambda)^{k-1}\). Thus \(k=1\) gives a constant hazard, \(k>1\) an increasing hazard, and \(k<1\) a decreasing hazard. For a log-normal model, \(\log T\sim N(\mu,\sigma^2)\), so \(S(t)=1-\Phi((\log t-\mu)/\sigma)\) for \(t>0\). Its hazard is

where \(\phi\) and \(\Phi\) are the standard normal density and distribution functions. The log-normal model is an accelerated failure-time model; its hazard need not be proportional over time and can rise and later fall.

From Transition Hazards to Occupancy Predictions

For the three-state illness-death process, arrange the transition hazards in the intensity matrix \(Q(t)\). Its off-diagonal entries are the allowed transition hazards and each diagonal entry is the negative sum of the other entries in that row. For example,

If \(P(u,t)\) is the matrix of transition probabilities from each starting state at time \(u\) to each state at time \(t\), its rows evolve according to \(dP(u,t)/dt=P(u,t)Q(t)\), with \(P(u,u)\) equal to the identity matrix. Because Weibull and log-normal hazards usually vary with time, this system generally requires numerical integration. Evaluate the fitted transition hazards at successive times, form \(Q(t)\), and integrate the matrix equation. Starting from state 1, the first row of \(P(0,t)\) gives the predicted probabilities of being illness-free, ill, or dead at time \(t\).

A useful derivation explains why the rows of the intensity matrix sum to zero. Over a short interval \(\Delta t\), the probability of a transition \(r\to s\) is \(q_{rs}(t)\Delta t+o(\Delta t)\). The probability of remaining in state \(r\) is therefore \(1-\sum_{s\ne r}q_{rs}(t)\Delta t+o(\Delta t)\). Subtracting the identity matrix, dividing by \(\Delta t\), and taking the limit gives the diagonal entry \(q_{rr}(t)=-\sum_{s\ne r}q_{rs}(t)\). This construction preserves total probability when the differential equations are integrated.

Model Checks, Comparisons, and Limitations

Before interpreting predictions, examine each transition fit separately. Check that the event count and person-time are plausible, that risk intervals begin and end at the recorded state changes, and that no transition is assigned to someone who was not in its origin state. Compare fitted survival or cumulative hazard patterns with nonparametric summaries for the corresponding risk set. Inspect fitted hazards over the observed time range, and assess whether their shape is credible where data are sparse. Convergence, extreme parameter estimates, or fitted hazards that change implausibly near the edge of follow-up are warnings, not merely software details.

The approaches answer related but distinct questions. A cause-specific Cox model leaves the baseline hazard unspecified and estimates covariate effects under a proportional-hazards structure; a parametric model specifies the baseline hazard family and can yield smooth extrapolations. The Aalen–Johansen estimator estimates occupancy without imposing these parametric shapes, but does not by itself provide the same distribution-based extrapolation. Weibull is especially restrictive when the true hazard is non-monotone; log-normal allows a non-monotone hazard but still imposes a particular distributional form. Neither model should be selected solely because it yields a convenient prediction.

This construction also has substantive limits. Predictions require that all relevant transitions are represented, that the fitted transition hazards are applicable to the target covariate profile, and that the Markov assumption is adequate. If future risk depends on time spent in the current state or other past history beyond what is modeled, the clock-forward Markov intensity matrix may not be adequate; the clock-reset semi-Markov approach is a different model. Extrapolation beyond observed follow-up is especially sensitive to parametric assumptions. Finally, independently fitted transition models do not make uncertainty in the assembled probabilities disappear: uncertainty in each fitted model should be propagated, for example by simulation or bootstrap, when intervals for predictions are required.