Stepwise Dynamic Competing Risks Model

• The stepwise dynamic competing risks model is a framework that updates risk at each transition as new event, time-scale, or longitudinal data become available.
• It incorporates various formulations, including Markov, semi-Markov, and non-Markov models, to handle sequential non-terminal and terminal events.
• The approach leverages landmark and joint modeling techniques to refine state occupation probabilities and cumulative incidence estimates in real time.

Searching arXiv for the cited paper and closely related work to support the article. arXivSearchTool.search({"query":"id:(Llopis-Cardona et al., 2021) OR id:(Liu et al., 2019) OR id:(Carollo et al., 2024) OR id:(Shen et al., 2023) OR id:(Niekerk et al., 2019) OR id:(Alvares et al., 2024) OR id:(Tong et al., 2021) OR id:(Allignol et al., 2013) OR id:(Wu et al., 2019) OR id:(Li et al., 2023) OR id:(Yu et al., 2023) OR id:(Alberge et al., 2024) OR id:(Albu et al., 2024) OR id:(Choi et al., 2020) OR id:(Rytgaard et al., 2021) OR id:(Zhang et al., 2018) OR id:(Nemchenko et al., 2018)","max_results":16}) The “Stepwise Dynamic Competing Risks Model” is an Editor’s term for a class of competing-risks formulations in which risk is updated as new event information, time scales, or longitudinal data accrue. In its simplest form, competing risks models describe one initial state with mutually exclusive terminal events, and observation stops at the first event. In more complex settings, especially when non-terminal events may occur sequentially, multi-state models provide the necessary extension by representing individuals as moving across several states over time with transition-specific hazards and probabilities (Llopis-Cardona et al., 2021). This framework includes Markov, semi-Markov, and non-Markov formulations; landmark proportional subdistribution hazards models; joint models with longitudinal biomarkers; recurrent-event constructions; and two-time-scale hazard surfaces, all of which operationalize dynamic updating of cumulative incidence and state occupation probabilities under competing events (Liu et al., 2019).

1. Conceptual scope and model class

Classical competing risks models are designed for settings with one initial state and mutually exclusive terminal events. Observation stops at the first event, and the occurrence of one event precludes observation of the others. They are therefore adequate when the events of interest are terminal and the process ends at the first observed event (Llopis-Cardona et al., 2021). By contrast, when non-terminal events such as disease onset, refracture, or illness can be followed by further events such as death, a multi-state representation is required because classical competing risks treat the first event as terminal by design and cannot represent subsequent transitions or hazards conditional on time since the non-terminal event (Llopis-Cardona et al., 2021).

Within this broader class, the defining “stepwise” feature is that risk changes at each transition as new information accrues. Multi-state models implement this by defining transition-specific hazards that can depend on prior events, current state, and time since entry into the current state. In a Markov model, the transition intensity depends only on current time and current state; in a semi-Markov model, it depends on current state and duration since entry into that state; and in a non-Markov model, it depends on richer history such as number of prior recurrences or covariate trajectory (Llopis-Cardona et al., 2021). This suggests that the phrase denotes a modeling strategy rather than a single standardized estimator.

Landmark formulations provide a complementary stepwise mechanism. At a landmark time $s$, the model conditions on being event-free at $s$ and updates the conditional cumulative incidence function using covariates observed up to $s$ (Liu et al., 2019). Dynamic prediction in this sense is sequential re-estimation or reconditioning over a grid of landmark times. Related developments use variable-length longitudinal histories, such as EEG summaries, to update competing-risk predictions as new information becomes available (Shen et al., 2023).

2. State-space representations and hazard formulations

A standard multi-state formulation begins with an explicit state graph. A typical example uses states $1=\text{healthy/initial}$, $2=\text{event A}$, $3=\text{event B}$, and $4=\text{death}$, with allowed transitions $1\to2$, $1\to3$, $1\to4$, $s$0, and $s$1, and optionally $s$2 or $s$3 if switching between non-terminal states is possible (Llopis-Cardona et al., 2021). In this setup, death is absorbing, whereas the non-terminal states are transient.

The general transition intensity for $s$4 at time $s$5 given history $s$6 and covariates $s$7 is

$s$8

A Cox-type Markov specification is

$s$9

whereas the semi-Markov extension replaces calendar time by time since entry into the current state,

$s$0

An illness-death example is the post-illness death hazard $s$1, which depends on elapsed time since illness onset (Llopis-Cardona et al., 2021).

For terminal competing risks, the canonical quantities are the cause-specific hazard

$s$2

overall survival

$s$3

and the cumulative incidence function

$s$4

Fine–Gray models instead target the cumulative incidence directly through a subdistribution hazard that keeps failures from other causes in the risk set with modified weights (Llopis-Cardona et al., 2021). This is appropriate when outcomes are terminal competing events and the goal is direct modeling of the cumulative incidence, but it is not designed for sequential non-terminal events or path probabilities after a non-terminal event (Llopis-Cardona et al., 2021).

A related generalization replaces a single time axis by two jointly evolving time scales. In the formulation with age at diagnosis $s$5 and time since diagnosis $s$6, the cause-specific hazard surface is $s$7, overall survival is

$s$8

with $s$9, and the cumulative incidence is

$1=\text{healthy/initial}$0

On a regular Lexis grid, these hazards can be treated as piecewise constant within grid cells and then smoothed by tensor-product $1=\text{healthy/initial}$1-splines (Carollo et al., 2024).

3. Estimation, inference, and dynamic prediction

Nonparametric estimation in multi-state competing-risks settings is centered on the Aalen–Johansen estimator,

$1=\text{healthy/initial}$2

where $1=\text{healthy/initial}$3 contains the Nelson–Aalen increments for each observed transition type at time $1=\text{healthy/initial}$4. Under independent censoring, this yields consistent estimates of state occupation probabilities and path probabilities (Llopis-Cardona et al., 2021). Transition-specific Cox models fit separate hazards for each allowed transition using transition-specific risk sets, optionally with time-dependent covariates (Llopis-Cardona et al., 2021).

Likelihood-based formulations generalize this to full multi-state paths. Under a multiplicative intensity model, the subject-level likelihood combines hazard contributions at observed transitions with an exponential term involving the integrated sum of outgoing hazards from the current state (Llopis-Cardona et al., 2021). Bayesian versions place priors on baseline hazards and regression coefficients and use posterior samples to compute transition probabilities and cumulative incidence functions, for example with MCMC or INLA (Llopis-Cardona et al., 2021, Niekerk et al., 2019).

Dynamic prediction is usually expressed through state occupation probabilities

$1=\text{healthy/initial}$5

or through path-specific probabilities. For an illness-death structure, the probability of remaining in the initial state is

$1=\text{healthy/initial}$6

and post-event survival in the illness state may depend on time since the illness transition via a semi-Markov hazard (Llopis-Cardona et al., 2021). After observing a non-terminal event at time $1=\text{healthy/initial}$7, subsequent hazards are updated to reflect time since event and new covariate values, so that prediction becomes explicitly conditional on the realized path (Llopis-Cardona et al., 2021).

Landmark Fine–Gray models offer a different route to dynamic prediction. At landmark time $1=\text{healthy/initial}$8, the subdistribution hazard for cause $1=\text{healthy/initial}$9 is modeled on residual time $2=\text{event A}$0,

$2=\text{event A}$1

leading directly to the conditional cumulative incidence over $2=\text{event A}$2 (Liu et al., 2019). The “landmark PSH supermodel” stacks landmark datasets and models the landmark dependence of coefficients as smooth functions of $2=\text{event A}$3, enabling one-step estimation across multiple landmarks (Liu et al., 2019). A related landmark subdistribution hazard approach with multiple longitudinal biomarkers uses local polynomial smoothing to estimate $2=\text{event A}$4 without explicitly imputing biomarker values at the landmark or modeling their trajectories (Wu et al., 2019).

4. Extensions with longitudinal markers, recurrent events, and alternative learning objectives

Joint modeling extends stepwise competing-risks prediction by linking event risk to evolving biomarkers. In one formulation for multiple myeloma, two biomarkers—M-spike and involved free light chains—are modeled by bi-exponential mixed-effects trajectories within each line of therapy, while competing risks of death and transition to the next line of therapy are described by Weibull cause-specific hazards whose linear predictors depend on the log-transformed bi-exponential parameters (Alvares et al., 2024). Predictions are updated at line-of-therapy initiation and within-line landmarks by integrating over the joint posterior of model parameters and random effects (Alvares et al., 2024). The paper reports that corrected two-stage estimation achieved similar predictive performance to joint estimation with materially reduced computation time, approximately a $2=\text{event A}$5 reduction, from $2=\text{event A}$6 hours to $2=\text{event A}$7 hours (Alvares et al., 2024).

A related backward joint model, crBJM, factorizes the likelihood into the distribution of competing-risks data and the distribution of longitudinal data given the competing-risks data. This yields dynamic risk prediction through one-dimensional integrals over latent future event times and event types, while avoiding the high-dimensional random-effects integration typical of shared-random-effects joint models with multiple biomarkers (Li et al., 2023). The model also allows prediction of future longitudinal trajectories conditional on remaining at risk at a future time (Li et al., 2023).

Recurrent-event settings require additional state augmentation. In a joint dynamic model for recurrent competing risks, longitudinal marker states, and health-status states, all intensities are updated multiplicatively as functions of the current health state, marker state, covariates, past event counts, and effective ages (Tong et al., 2021). Here the stepwise dynamics arise because every event or transition updates the relevant states or virtual ages, which in turn instantaneously modifies all intensities (Tong et al., 2021). In non-Markov illness-death settings, modified competing-risks estimators can be used to estimate transition probabilities nonparametrically when the post-illness hazard depends on time since illness rather than only current calendar time (Allignol et al., 2013).

Machine-learning variants target the same prediction problem with different objectives. Dynamic-DeepHit and a modified DDRSA estimate patient-level cumulative incidence functions on a discretized time grid from variable-length EEG histories under three competing risks: awakening, death not due to withdrawal of life-sustaining therapy, and withdrawal of life-sustaining therapy (Shen et al., 2023). In a cohort of $2=\text{event A}$8 post-cardiac-arrest coma patients, the classical Fine–Gray model using static features and the latest hour’s EEG summary was highly competitive with Dynamic-DeepHit; for awakening at landmark $2=\text{event A}$9 hours and horizon $3=\text{event B}$0 hours, the reported time-dependent c-indices were $3=\text{event B}$1 for Fine–Gray and $3=\text{event B}$2 for Dynamic-DeepHit (Shen et al., 2023). Other approaches directly optimize proper censoring-adjusted multiclass scores for cumulative incidence, as in MultiIncidence, or optimize pairwise time-dependent concordance surrogates, as in SSPN (Alberge et al., 2024, Nemchenko et al., 2018).

5. Applied illustrations and empirical patterns

A frequently cited applied example is recurrent hip fracture. The state space comprises discharge after index hip fracture ($3=\text{event B}$3), recurrent hip fracture ($3=\text{event B}$4), and death ($3=\text{event B}$5), with transitions $3=\text{event B}$6, $3=\text{event B}$7, and $3=\text{event B}$8 (Llopis-Cardona et al., 2021). A competing-risks model for $3=\text{event B}$9 and $4=\text{death}$0 and an illness-death model for $4=\text{death}$1, $4=\text{death}$2, and $4=\text{death}$3 were fitted using Bayesian Cox-type hazards with Weibull baselines and covariates age and sex (Llopis-Cardona et al., 2021). For the common transitions $4=\text{death}$4 and $4=\text{death}$5, both models yielded nearly identical cumulative incidences, whereas the illness-death model additionally provided state occupation probabilities and estimates of death after refracture (Llopis-Cardona et al., 2021). One-year death after refracture was reported as about $4=\text{death}$6 at age $4=\text{death}$7 and $4=\text{death}$8 at age $4=\text{death}$9 for men, and $1\to2$0 and $1\to2$1 for women (Llopis-Cardona et al., 2021).

In a two-time-scale competing-risks analysis of breast cancer mortality, hazards were modeled over age at diagnosis and time since diagnosis using two-dimensional $1\to2$2-splines and generalized linear array models (Carollo et al., 2024). In the SEER application, age at diagnosis was available with a last open-ended category and was ungrouped using a two-dimensional penalised composite link model before fitting the hazards (Carollo et al., 2024). The resulting hazard and cumulative incidence surfaces revealed patterns that depended jointly on age at diagnosis and time since diagnosis, including non-proportional chemotherapy effects across the two time scales (Carollo et al., 2024).

In multiple myeloma, the blockwise line-of-therapy joint model found that higher baseline M-spike and free light chains increased death hazard across lines of therapy, while baseline and growth of both biomarkers increased the hazard of transitioning to the next line of therapy (Alvares et al., 2024). Dynamic discrimination was satisfactory in earlier lines of therapy and deteriorated in later lines, where data were sparser and biomarker informativeness diminished (Alvares et al., 2024).

In hospital epidemiology, a dynamic random-forest comparison for central line-associated bloodstream infection treated discharge and death as competing events and evaluated daily landmark predictions for 7-day risk (Albu et al., 2024). The data included $1\to2$3 admissions, $1\to2$4 catheter episodes, $1\to2$5 CLABSI events, $1\to2$6 deaths, and $1\to2$7 discharges (Albu et al., 2024). Binary, multinomial, and competing-risks models performed similarly, with AUROC about $1\to2$8 for baseline predictions and about $1\to2$9 around day $1\to3$0, whereas survival models censoring competing events at their occurrence overestimated the CLABSI risk with Expected:Observed ratios between $1\to3$1 and $1\to3$2 (Albu et al., 2024).

6. Assumptions, software, and methodological boundaries

The standard assumptions recur across formulations: independent or non-informative right censoring, correct specification of the allowed transition structure and time scales, and an appropriate choice among Markov, semi-Markov, and non-Markov dependence structures (Llopis-Cardona et al., 2021). When duration in state matters, semi-Markov formulations are often favored; when full path history matters, non-Markov structures become necessary (Llopis-Cardona et al., 2021). This suggests that a major misconception is to treat all sequential event settings as adequately handled by first-event competing-risks models.

Model checking depends on the formulation. For transition-specific proportional hazards models, Schoenfeld residuals remain standard; for multi-state models, comparison of Aalen–Johansen estimates with model-based predictions is recommended; for Markov assumptions, dependence on time since entry should be examined; and in Bayesian implementations, posterior predictive checks for state occupation and transition counts are appropriate (Llopis-Cardona et al., 2021). In landmark models, calibration through observed-to-expected ratios, time-dependent Brier scores, and time-dependent AUC has been emphasized (Liu et al., 2019, Wu et al., 2019).

A concise summary of typical implementation choices is useful.

Setting	Typical formulation	Representative tools
Terminal first events	Cause-specific Cox or Fine–Gray	survival, cmprsk
Sequential non-terminal events	Multi-state Cox, parametric, or Bayesian model	mstate, msm, flexsurv, INLA
Joint longitudinal and competing risks	Joint model or backward joint model	rstan, Stan, R-INLA
Two time scales	Lexis-grid competing-risks with $1\to3$3-splines	TwoTimeScales

Frequentist and Bayesian software ecosystems are both established. For multi-state analysis, commonly listed tools are mstate, msm, flexsurv, and survival; for Fine–Gray, cmprsk; and for Bayesian latent Gaussian formulations, INLA (Llopis-Cardona et al., 2021). Two-time-scale competing-risks smoothing is implemented in the TwoTimeScales package (Carollo et al., 2024). Joint competing-risks models using latent Gaussian structures can be fitted in R-INLA, including non-Gaussian longitudinal data, spatial structures, time-dependent splines, and several association structures (Niekerk et al., 2019).

Overall, the stepwise dynamic competing-risks framework is best understood as a unifying statistical strategy: encode clinically or scientifically plausible event pathways, represent evolving risk through transition-specific or landmark-specific quantities, and update predicted cumulative incidence or state occupation probabilities whenever new event or biomarker information becomes available. Classical competing risks remain appropriate when events are terminal and first-event cumulative incidence is the only target; once non-terminal events, recurrent events, multiple time scales, or longitudinally updated prediction become central, the broader stepwise dynamic formulation is required (Llopis-Cardona et al., 2021).