Fitting Transition-Specific Cox Models with mstate

From Expanded Histories to Transition-Specific Models

After multi-state histories have been expanded, a subject may contribute several rows: one for each transition-specific risk interval in which that subject is at risk. The next task is to describe how covariates relate to the intensity of each transition. A transition-specific Cox model allows both the baseline transition intensity and, when needed, the covariate effects to differ by transition.

This is not a model for a single time until “any event.” Each transition has its own event indicator and risk set. A subject who moves from one state to another may enter a later transition risk set, while a subject who has not yet entered the relevant state is not at risk for that transition.

In the expanded data, the usual counting-process fields are Tstart, Tstop, and status, together with a transition identifier such as trans. A row represents a subject at risk for one transition over \((\text{Tstart},\text{Tstop}]\); status is 1 if that transition occurs at the interval end and 0 otherwise. The data also retain subject identifiers and covariates. Check the actual field names and coding in the object produced by msprep before fitting.

The model assumes that, within each transition, the covariate effect is proportional over the selected time scale unless time-varying effects are specified. The risk intervals and event codes must correctly represent the observed histories; censoring must be conditionally non-informative for the transition intensity being modeled; and subjects must be independent for the ordinary likelihood-based analysis. Covariates must be measured and encoded appropriately, including their values over each interval if they change over time. A Markov assumption is not required just to fit these conditional transition-intensity models; assumptions needed for particular prediction targets are a separate issue.

Specifying Separate Effects in coxph

In coxph, Surv(Tstart, Tstop, status) identifies the counting-process response. Including strata(trans) gives each transition its own baseline hazard. That alone does not give each transition its own covariate effects: an un-interacted covariate in the formula has a common coefficient across transitions.

For transition-specific effects, mstate::expand.covs creates transition-indexed covariate columns. For example, expanding a binary treatment variable creates a treatment column for each transition. Including those columns together with strata(trans) permits a separate treatment coefficient for each transition and a separate baseline hazard for each transition. Use the names actually produced by expand.covs; they depend on its options and the transition indexing.

The Breslow tie method is commonly used in this workflow and is appropriate to specify explicitly when the fitted object will be passed to mstate::msfit. If there are repeated records per subject, the rows are not independent observations. For robust subject-level inference, a subject identifier can be supplied as a cluster in coxph, provided the downstream procedures support the resulting variance specification. Clustering changes the variance estimate, not the fitted coefficient estimates. Confirm the model and variance structure required by the analysis rather than adding a cluster term mechanically.

Why the Transition-Specific Likelihood Works

For transition \(r\), let \(Y_{ir}(t)\) indicate whether subject \(i\) is at risk for that transition at time \(t\), and let \(N_{ir}(t)\) count occurrences of that transition. Under the transition-specific proportional-hazards model, the contribution to the partial likelihood is formed from event times for that transition and the subjects in its risk set. With no tied event times, its factor for event subject \(i\) at time \(t\) is the subject’s relative intensity divided by the sum of relative intensities over the risk set.

Here \(\mathcal{D}_r\) is the set of observed event times for transition \(r\), and \(i(t)\) is the subject experiencing that transition at time \(t\). Separate baseline hazards cancel from this partial likelihood. When each transition has its own coefficients, the transition-specific likelihood factors can be maximized with their own covariate effects; the stratified Cox representation implements distinct baselines while the expanded covariate columns encode distinct slopes.

For tied events, the Breslow approximation at time \(t\) with \(d_r(t)\) events uses the denominator raised to the event count. The corresponding increment of the Breslow cumulative baseline hazard is \(\Delta\widehat H_{0r}(t)=d_r(t)/\sum_{j:Y_{jr}(t)=1}\exp\{\hat\beta_r^{\mathsf T}x_j(t)\}\). This makes explicit why a wrongly coded risk set or event count affects both the coefficient fit and the estimated baseline hazard.

Obtaining Transition Cumulative Hazards with mstate

The fitted coefficients summarize relative transition intensities. For later multi-state prediction, an absolute transition-hazard scale is also needed. The mstate::msfit function uses a fitted Cox model, a specified covariate profile, and the transition matrix to organize estimated cumulative transition hazards. It returns transition-indexed hazard information; it does not itself turn that information into state-occupation probabilities.

When transition effects use expanded columns, supply a profile consistent with the expansion. For a treatment profile with value 1, each transition-indexed treatment column should have value 1 for that profile; for value 0, each should have value 0. Include any other covariates required by the fitted formula, and use the same transition matrix that defines the transition indexing in the expanded data. Check that the supplied profile has all model variables and that the matrix and model use the same transition labels.

Diagnostics, Alternatives, and Limitations

A successful software call is not evidence that the model is correctly specified. Before interpretation, check that transition identifiers match the transition matrix, that each event is assigned to the correct transition, that the time intervals and event status are sensible, and that covariates do not change unexpectedly across rows. Inspect event counts by transition and the number of subjects at risk; sparse transitions can yield unstable or non-estimable coefficients. Look for warnings about convergence, singularity, or infinite estimates, and inspect coefficient standard errors and confidence intervals for signs of weak information.

The proportional-hazards assumption should be considered separately for each transition and covariate. Schoenfeld-residual checks and graphical or time-interaction assessments can reveal time-varying effects. A strong violation means that a constant transition-specific hazard ratio may be inadequate; a transition-specific time-varying effect or another model may be needed. Also examine whether covariate coding and missing-data handling are consistent across the expanded rows. One subject’s original covariate value should not accidentally be counted as a different value merely because the history was expanded.

The main alternative structures include a shared covariate effect across transitions, transition-stratified models with common slopes, and separate model fits for each transition. A common slope is more parsimonious but imposes equality of effects; separate fits are straightforward to inspect but require care to keep risk sets and covariate handling aligned. Choosing among these structures is a modeling decision, not a software default. Transition-specific models also do not automatically solve sparse-data problems, unmeasured confounding, informative censoring, or misspecification of the time scale.