The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
87,324 characters
Estimating Pathway Treatment Effects in the Presence of Intermediate Events with Multi-State Data
\def\spacingset#1{\renewcommand{\baselinestretch}
{#1}\small\normalsize} \spacingset{1}
\if11
{
\title{\bf Estimating Pathway Treatment Effects in the Presence of Intermediate Events with Multi-State Data}
\author{Yuhao Deng\thanks{
Yuhao Deng and Haoyu Wei contributed equally. The authors thank Kajsa Kvist from Novo Nordisk A/S (Denmark) for the insightful discussions on the LEADER Trial.}\hspace{.2cm}\\
Vaccine and Infectious Disease Division, Fred Hutchinson Cancer Center\\
Haoyu Wei$^*$ \\
Department of Economics, University of California San Diego \\
Donglin Zeng \\
Department of Biostatistics, University of Michigan \\
Rui Song \\
Amazon Inc. \\
Xiao-Hua Zhou \\
Department of Biostatistics and Beijing International Center for \\
Mathematical Research, Peking University
}
\date{\vspace{-20pt}}
\maketitle
} \fi
\if01
{
\bigskip
\bigskip
\bigskip
\begin{center}
{\LARGE\bf Estimating Pathway Treatment Effects in the Presence of Intermediate Events with Multi-State Data}
\end{center}
\medskip
} \fi
\bigskip
\begin{abstract}
During clinical trials evaluating a drug's effect on a survival endpoint, intermediate events often occur in addition to the primary event. The treatment can exert its effect on the primary endpoint along multiple pathways through intermediate events. Assumptions for identifying mediation effects, such as sequential ignorability in natural effects or the dismissible components condition in separable effects, fail because intermediate events act as treatment-induced confounding. To understand the effect along each pathway, we consider hypothetical interventions in transitions between event statuses to mimic the treatment mechanism. The hypothetical interventions adjust for effects through intermediate events and marginalize over unobserved treatment-induced confounding, if any. Based on the derived efficient influence functions for the counterfactual cumulative incidences under hypothetical interventions, we construct multiply robust and semiparametrically efficient estimators for pathway treatment effects. Our proposed framework enables the examination of treatment effects through each transition, on each event, and along each path. By analyzing data from the LEADER Trial, we find that liraglutide significantly reduces the risk of cardiovascular and microvascular events. The reduction in all-cause mortality is primarily mediated by its effects on expanded major adverse cardiovascular events.
\end{abstract}
\noindent
{\it Keywords:} Causal inference; Efficient influence function; Mediation; Multi-state model; Time-to-event; Transition.
\vfill
\newpage
\spacingset{1.8}
\section{Introduction} \label{sec:intro}
When evaluating the effect of a new drug on a survival endpoint, individuals may experience multiple intermediate events, thereby generating multiple transition pathways from baseline to the primary endpoint. The effects operating through each pathway provide insight into multiple questions, such as how the drug is effective and whether additional post-treatment interventions should be recommended. For example, to study the direct effect, researchers wish to isolate the effect of treatment itself on the primary endpoint, independent of any modification through intermediate events. By identifying the transitions associated with the largest treatment effects, doctors can provide intensive prompts to increase treatment adherence when patients are at risk of these transitions. By identifying the transition pathways through which treatment is ineffective or harmful, doctors can impose additional medications to prevent these transitions.
Treating the primary event and intermediate events as distinct states yields the so-called multi-state data \citep{hougaard1999multi, andersen2002multi, putter2007tutorial}. In particular, when the primary event is a terminal state and the intermediate events may or may not occur before it, this results in special cases called competing or semi-competing risks data, which have been extensively studied \citep{prentice1978analysis, kodell1980illness, fine1999proportional, fine2001semi, xu2010statistical, andersen2012competing, chen2012maximum}. For such data, mediation-type frameworks, such as natural effects and separable effects, have been proposed to estimate direct and indirect treatment effects on the primary event by counterfactually controlling for intermediate event risks, under no-unmeasured-confounding-like assumptions \citep{huang2021causal, weir2022counterfactual, martinussen2023estimation, deng2024direct, breum2024estimation}. However, these frameworks only considered at most one type of intermediate event. Since intermediate events act as treatment-induced confounders for subsequent transitions, natural effects are not identifiable due to violation of sequential ignorability if there are both semi-competing and competing events \citep{shpitser2016causal, miles2020semiparametric}. The separable effects framework posits that the original treatment consists of multiple components, each of which only directly influences a single event \citep{stensrud2021generalized, stensrud2022separable, robins2022interventionist}. Separable effects target the effect of each component, but are not informative about the effect through each transition pathway. Even if the effects of each component are theoretically well-defined, it is hard to identify biologically meaningful treatment components with isolated effects in a real trial. Moreover, unobserved treatment-induced confounding is unavoidable due to the complex interactions among events during disease progression, so the no-unmeasured-confounding-like assumptions underlying the identification of natural and separable effects cannot be satisfied.
To avoid these untestable or unverifiable assumptions, we consider hypothetical interventions associated with each transition. In particular, we envision a hypothetical scenario in which local interventions are imposed at each possible transition between events, while all other transitions remain unaffected. This intervention mimics the organic effect that shifts the risks of intermediate events \citep{lok2015defining}. Essentially, we cut one directed edge in the causal graph. Formally, we extend the randomized (stochastic) interventional effects framework, originally proposed for mediation analysis in longitudinal studies, to accommodate multi-state data \citep{vanderweele2017mediation, vansteelandt2017interventional, lin2017interventional, lin2017mediation, diaz2020causal, hejazi2022nonparametric, valeri2023multistate}. The randomized intervention envisions random draws of the intermediate event process from a specific distribution, either conditioning on or marginalizing over time-varying confounding \citep{deng2026randomized}. When unmeasured confounders are present, the randomized intervention remains valid in a hypothetical sense by marginalizing over them. Nevertheless, interpreting such an intervention is less straightforward: it mimics the generation of intermediate events under ``recanting twins'' of unmeasured confounders \citep{vo2026recanting}. Identifying treatment effects under hypothetical interventions does not require cross-world assumptions such as sequential ignorability as long as the hypothetical intervention is meaningful. The randomized interventional effects framework is powerful in studies with multiple mediators or multiple outcomes. For example, path-specific effects can be defined by considering sequences of randomized interventions \citep{vansteelandt2019Mediation, Diaz2022Causal, tai2023causal, Kormaksson2024Dynamic}.
There are several advantages to evaluating treatment effects through hypothetical interventions. First, by intervening in a single transition and comparing changes in cumulative incidence functions, one can determine whether the treatment affects this transition. Second, by intervening in the transitions into a state, one can identify which events are affected by the treatment. Third, by intervening in all the transitions along a path, one can estimate the path-specific treatment effect. Estimating transition-, event-, and path-specific effects generally concerns estimating counterfactual cumulative incidences under hypothetical interventions imposed on transitions. Rather than intervening in event times, the key to distinguishing pathway effects in multi-state data is to intervene in the instantaneous risk of transitions, as characterized by hazard functions. For example, a practitioner may consider exposing healthy individuals to a high-risk environment for cardiovascular events, thereby introducing a hypothetical transition-specific hazard. Next, the practitioner may expose individuals to a low-risk environment for microvascular events, thereby introducing another transition-specific hazard. The levels of transition-specific hazards in the hypothetical world can be specified as the observable transition hazards in the real world, regardless of treatment-induced confounding.
However, hypothetical interventions for multi-state data pose conceptual challenges for randomized interventions. First, interventions can only be performed for at-risk individuals who have not experienced the terminal event. Second, once an intermediate event occurs under the intervention, subsequent interventions should not cause this event to recur. Third, the interventions should be applied sequentially across all event-counting processes, not only to a single intermediate event, because the status of any event can modify the risk of another. In addition to the conceptual difficulty of defining the randomized intervention, the interaction among intermediate events makes estimation more complex than modeling competing risks. While cumulative incidence functions can be calculated from transition hazards using the Kolmogorov forward equation \citep{andersen2002multi}, such a plug-in estimator is inefficient and model-dependent. If the working models are misspecified, the resulting estimator may be biased, yielding misleading conclusions. Statistically efficient and robust estimation methods remain understudied.
In this paper, we define randomized interventional effects for multi-state data. We rigorously specify the sequential randomized interventions in event-counting processes that control transition hazards by marginalizing over potential treatment-induced confounding. We demonstrate the identifiability of the counterfactual cumulative incidence of each event under any given intervention. The total effect on an event can be decomposed into the sum of interventional effects arising from sequentially intervening at each transition. We derive the efficient influence function for the counterfactual cumulative incidence under semiparametric theory \citep{bickel1993efficient} and provide a semiparametrically efficient, multiply robust estimator that does not require modeling of treatment-induced confounding. We establish the asymptotic properties for the proposed estimator. We also propose hypothesis testing methods for treatment effects.
\section{Randomized Interventional Effects} \label{sec:frame}
\subsection{Notations and Estimands}
Suppose there are $K$ states, including initial (baseline), intermediate, and terminal states. The most typical multi-state data involve only one initial state (denoted by O), such as competing or semi-competing risks data. The transitions between states can be represented by a graph $(\mathcal{G}, \mathcal{E})$, where $\mathcal{G}$ is the set of states (nodes) and $\mathcal{E}$ is the set of one-step transitions (edges) between adjacent states. Let $\underline{\mathcal{G}}_k$ represent the set of states with one-step transitions from state $k$ to its adjacent following states. We denote the set of paths from state $k$ to state $g$ as $\mathcal{Q}_{kg}$, and the set of paths from initial states to state $g$ as $\mathcal{Q}_g$. A path $q \in \mathcal{Q}_{kg}$ is represented as an ordered sequence of states $q = (q(0), q(1), \ldots, q(l(q)))$, where $q(0) = k$, $q(l(q)) = g$, and $l(q)$ denotes the number of transitions along $q$. We assume that each state is visited at most once on any path. For example, Figure \ref{fig:model} shows the multi-state structure for our motivating data, where O denotes the initial state, E denotes expanded major adverse cardiovascular events, M denotes microvascular events, and D denotes death. For this multi-state structure, $\mathcal{Q}_{\text{D}} = \mathcal{Q}_{\text{OD}} = \{q_1,q_2,q_3,q_4,q_5\}$. For a single individual, it can either remain in the initial state O, or reach an event $g\in\{\text{E},\text{M},\text{D}\}$ along a single path in $\mathcal{Q}_{\text{O}g}$ at the end of the study.
\begin{figure}[!tb]
\centering
\includegraphics[width=0.8\textwidth]{Figure1.pdf}
\caption{Five potential paths from baseline (O) to death (D) with expanded major adverse cardiovascular events (E) and microvascular events (M) as two intermediate states.}
\label{fig:model}
\end{figure}
Transition between event statuses is the core of modeling multi-state data, and treatment effect evaluation should start from each transition. An interpretable causal estimand should satisfy two rules:
\begin{itemize}
\item The effect on a transition only matters on the population at risk of the transition, so the treatment should shift the transition hazard.
\item The effect on a transition is only present after conditioning on the history before the transition; since latent variables such as frailty are unobserved, the hazard should marginalize over latent variables.
\end{itemize}
Specifically, for each transition $(k,g) \in \mathcal{E}$, we generate $|\underline{\mathcal{G}}_{k}|$ hypothetical counting processes are according to the hazards of $\underline{\mathcal{G}}_{k}$, and the first jump among these counting processes determines the next state an individual will transition to. In the counterfactual world, there are $|\mathcal{E}|$ transitions that are subject to intervention. Let $\overline{a} = (a_{kg}: (k,g)\in\mathcal{E})$ denote the levels of interventions, where the entries of $\overline{a}$ may take different values. Here, $a_{kg} = 1$ indicates that the transition from $k$ to $g$ is intervened in by drawing a counting process from the natural level under treatment, and $a_{kg} = 0$ indicates that it is intervened in by drawing a counting process from the natural level under control. In this way, interventions are sequential from one transition to the next, resulting in a hypothetical transition trajectory.
Let $W \in \mathbb{R}^p$ be baseline covariates with support $\mathcal{W}$. The end of the study is set at $t^*$. In practice, $t^*$ should not be larger than the maximum follow-up time by design. Mathematically, we can characterize the above hypothetical transitions in terms of conditional hazard functions as described below. Let $X^{\overline{a}}(t)$ represent the current state at time $t$ under the interventions of level $\overline{a}$. We define the hypothetical transition trajectory up to time $t$ as $\overline{X}^{\overline{a}}(t) = (X^{\overline{a}}(s) \in \mathcal{G}: 0 \leq s \leq t)$. This hypothetical transition trajectory is generated in the following way. Let $X^a(t)$ be the potential state at time $t$ if the unit is assigned treatment $a$ at baseline, and let
\begin{equation*}
\mathrm{d} \Lambda_{kg}^{a}(t \mid W,\overline{X}^{a}(t)) := \mathrm{P} \{X^{a}(t + \mathrm{d} t)=g \mid X^a(t)=k, W,\overline{X}^{a}(t)\}
\end{equation*}
be the potential transition hazard. For the transition $(k,g) \in \mathcal{E}$, we intervene in the counting process of event $g$ according to the reference hazard associated with treatment $a_{kg} \in \{0,1\}$. Then the counting process of $g$ given history $(W,\overline{X}^{\overline{a}}(t))$ subject to $X^{\overline{a}}(t)=k$ under the randomized intervention is drawn from the hazard
\begin{equation} \label{int_haz}
\mathrm{d} \Lambda_{kg}^{\overline{a}}(t \mid W,\overline{X}^{\overline{a}}(t)) := \mathrm{P} \{X^{\overline{a}}(t + \mathrm{d} t)=g \mid X^{\overline{a}}(t)=k,W,\overline{X}^{\overline{a}}(t)\}.
\end{equation}
Such a randomized intervention postulates that the counterfactual hazards of transitions between two adjacent states are locally stable under intervention. Notably, the randomized intervention remains hypothetically achievable even in the presence of unmeasured confounding. In the presence of potential treatment-induced confounding $\overline{L}^{a}(t)$, the hazard \eqref{int_haz} marginalizes over treatment-induced confounding,
\begin{equation}\label{draw_uc}
\mathrm{d} \Lambda_{kg}^{\overline{a}}(t \mid W,\overline{X}^{\overline{a}}(t)) := \int_{\overline{l}(t)} \mathrm{d} \Lambda_{kg}^{a_{kg}}(t \mid W,\overline{l}(t),\overline{X}^{\overline{a}}(t)) \mathrm{P} (\overline{L}^{a_{kg}}(t)=\overline{l}(t) \mid W, \overline{X}^{\overline{a}}(t)) \mathrm{d} \overline{l}(t).
\end{equation}
The primary object of interest is the time to reach a specific state $j \in \mathcal{G}$ (typically a terminal event such as death), defined as $T_j^{\overline{a}} = \inf\{t \in [0, t^*]: X^{\overline{a}}(t) = j \}$ if state $j$ is on the trajectory $\overline{X}^{\overline{a}}(t^*)$, otherwise we set $T_j^{\overline{a}} = \infty$.
The primary estimand is the counterfactual cumulative incidence function (CIF) of state $j \in \mathcal{G}$, given by
\begin{equation}\label{target_par_all}
F_j^{\overline{a}}(t) := \mathrm{P} ( T_j^{\overline{a}} \leq t ) = \mathrm{P}(j \in \overline{X}^{\overline{a}}(t)), \ t \in [0, t^*].
\end{equation}
This estimand evaluates the probability that the event $j$ occurs by time $t$ under this hypothetical intervention. It remains well-defined even if $T^{\overline{a}}_j$ is infinite. Note that the paths from baseline to state $j$ are not unique. Let $F_j^{\overline{a}}(t;q)$ denote the sub-incidence (or path-specific cumulative incidence) of state $j$ along a specific path $q \in \mathcal{Q}_j$. Then, the counterfactual cumulative incidence function can be decomposed as
\begin{equation}\label{par_relationship}
F_j^{\overline{a}}(t) = \sum_{q \in \mathcal{Q}_{j}} F_j^{\overline{a}}(t; q).
\end{equation} The counterfactual CIF is the key to defining transition-, event-, and path-specific effects, as these effects are defined by comparing counterfactual CIFs under different hypothetical interventions.
To illustrate the interventions and estimands in our framework, we consider the multi-state data structure in Figure \ref{fig:model}. The level of interventions is indexed by a vector of seven dimensions,
\[
\overline{a} = (a_{\text{OE}}, a_{\text{OM}}, a_{\text{OD}}, a_{\text{EM}}, a_{\text{ME}}, a_{\text{ED}}, a_{\text{MD}}).
\]
For example, we can impose treatment to the transitions along the path $\text{O}\to\text{E}\to\text{D}$, while leaving all other transitions under control, that is, $\overline{a} = (1,0,0,0,0,1,0)$. Figure \ref{fig:dag_int} presents the causal graph of the counterfactual event-counting processes: $N_{\text{E}}^{\overline{a}}(t)$ for E, $N_{\text{M}}^{\overline{a}}(t)$ for M, and $N_{\text{D}}^{\overline{a}}(t)$ for D. Given the history up to time $t$, if the individual is in the initial state O, we draw $N_{\text{OE}}^{\overline{a}}(t)$ based on the transition hazard of O$\to$E under treated, $N_{\text{OM}}^{\overline{a}}(t)$ based on the transition hazard of O$\to$M under control, and $N_{\text{OD}}^{\overline{a}}(t)$ based on the transition hazard of O$\to$D under control; if the individual is in state E, we draw $N_{\text{EM}}^{\overline{a}}(t)$ based on the transition hazard of E$\to$M under control and $N_{\text{ED}}^{\overline{a}}(t)$ based on the transition hazard of E$\to$D under treated; if the individual is in state M, we draw $N_{\text{ME}}^{\overline{a}}(t)$ based on the transition hazard of M$\to$E under control and $N_{\text{MD}}^{\overline{a}}(t)$ based on the transition hazard of M$\to$D under control.
The event-counting processes for E, M, and D are $N_{\text{E}}^{\overline{a}}(t) = \max\{N_{\text{OE}}^{\overline{a}}(t),N_{\text{ME}}^{\overline{a}}(t)\}$, $N_{\text{M}}^{\overline{a}}(t) = \max\{N_{\text{OM}}^{\overline{a}}(t),N_{\text{EM}}^{\overline{a}}(t)\}$, and $N_{\text{D}}^{\overline{a}}(t) = \max\{N_{\text{OD}}^{\overline{a}}(t),N_{\text{ED}}^{\overline{a}}(t),N_{\text{MD}}^{\overline{a}}(t)\}$
by aggregating all transitions to this event, respectively.
The counterfactual cumulative incidence function of D is $F_{\text{D}}^{\overline{a}}(t) = \mathrm{P}\{N_{\text{D}}^{\overline{a}}(t)=1\}$.
\begin{figure}[!tb]
\centering
\includegraphics[width=0.5\textwidth]{Figure2.pdf}
\caption{Causal graphs of counterfactual counting processes. Solid arrows represent stochastic effects on the counting processes of E, M, and D. Dashed arrows represent deterministic effects, enforcing that $N_{\text{E}}^{\overline{a}}(t+s) = N_{\text{E}}^{\overline{a}}(t)$ and $N_{\text{M}}^{\overline{a}}(t+s) = N_{\text{M}}^{\overline{a}}(t)$ for all $s>0$ if $N_{\text{D}}^{\overline{a}}(t) = 1$ at some $t>0$.} \label{fig:dag_int}
\end{figure}
\subsection{Identifiability}
According to the Kolmogorov forward equation in multi-state models, the central goal of identifying the counterfactual cumulative incidence function is to identify transition hazards. In the hypothetical world, by definition, the intervention leads to transition hazards exactly identical to the potential hazards. Therefore, we only need assumptions to identify potential hazards. The following assumptions are commonly made in causal inference methods for survival data.
\begin{assumption}[Ignorability] \label{ass:ign}
$A \protect\mathpalette{\protect\independenT}{\perp} \overline{X}^{a}(t^*) \mid W$ for $a=0,1$.
\end{assumption}
Ignorability (also known as unconfoundedness or exchangeability) excludes unmeasured confounding between the treatment assignment and the potential outcomes. This assumption does not exclude treatment-induced confounding. To account for censoring, let $T_{\text{c}}^a$ be the potential censoring time and $\mathrm{d}N^a_{\text{c}}(t)$ be the increase (jump) of the potential counting process for censoring in $[t,t+dt)$ under treatment assignment $a$. We assume censoring is independent of potential outcomes given the history.
\begin{assumption}[Random censoring]\label{ass:ran}
$ N_{\text{c}}^a(t) \protect\mathpalette{\protect\independenT}{\perp} \overline{X}^a(t+dt) \mid W, \overline{X}^{a}(t), A=a$.
\end{assumption}
Additionally, we assume that the treatment probability and the uncensored probability at the end of the study are both positive.
\begin{assumption}[Positivity] \label{ass:pos}
There exists a positive constant $\epsilon > 0$ such that $\mathrm{P}(A=a \mid W) > \epsilon$, and $\mathrm{P} (T_{\text{c}}^A \geq t^*, A=a \mid W, \overline{X}^{\overline{a}}(t^*)) > \epsilon$ on the support of $(W, \overline{X}^{\overline{a}}(t^*))$.
\end{assumption}
Positivity has two implications. First, the propensity score (i.e., the probability of receiving treatment) is bounded away from 0 and 1, ensuring we have data from both the treated and control groups.
Second, at the end of the follow-up period $t^*$, a non-negligible proportion of subjects remain at risk and uncensored within the strata defined by $W$ and $\overline{X}^{\overline{a}}(t^*)$.
Let $\overline{X}(t) = ( X(s): 0 \leq s \leq t )$ be the observed trajectory from the treatment initiation to time $t$. The trajectory $X(\cdot)$ is defined so that $X(t) = g$ if the individual is in state $g \in \mathcal{G}$ at time $t$ for $0 \leq t \leq t^*$. If the individual is censored at $T_{\text{c}}$, we denote $X(t) = \text{c}$ for $t>T_{\text{c}}$. We let $\Delta_k = 1$ if state $k$ is visited before the end of the study. Let $T_k = \inf\{t\in[0,t^*]: X(t)=k\}\wedge T_{\text{c}}\wedge t^*$ be the time of the first visit to state $k$ if $\Delta_k = 1$, and otherwise $T_k = T_{\text{c}}\wedge t^*$.
\begin{assumption}[Consistency] \label{ass:con}
$T_{\text{c}}^A = T_{\text{c}}$ and $\overline{X}^{A}(T_{\text{c}}^A\wedge t^*) = \overline{X}(T_{\text{c}}\wedge t^*)$.
\end{assumption}
Consistency links the potential outcomes with the observed data. Specifically, if an individual has not been censored at time $t$, then the potential transition trajectory under the observed treatment assignment $A$ coincides with the observed trajectory.
Assumptions \ref{ass:ign}--\ref{ass:con} ensure that the potential transition hazards can be identified as observed hazards,
\[
\mathrm{d} \Lambda_{kg}^a (t \mid w, \overline{x}(t)) = \mathrm{P} \{t < T_g < t+\mathrm{d} t, \Delta_g = 1 \mid A=a, W=w, \overline{X}(t) = \overline{x}(t), X(t) = k\}.
\]
Consequently, the identifiability of the counterfactual cumulative incidence $F_j^{\overline{a}}(t)$ of state $j\in\mathcal{G}$ under the hypothetical treatment sequence $\overline{a}$ is presented in the following theorem.
\begin{theorem} \label{thm1}
Under Assumptions \ref{ass:ign}--\ref{ass:con}, the counterfactual cumulative incidence $F_j^{\overline{a}}(t)$ of any event $j \in \mathcal{G}$ is identified as
\begin{align}
F_j^{\overline{a}}(t) = \sum_{q\in\mathcal{Q}_{j}}\int_{\mathcal{W}}\int_0^t\int_0^{t_{l(q)}}\cdots\int_0^{t_2} \prod_{r=1}^{l(q)} \mathrm{d} \mathrm{P}_r(t_r \mid a_{q(r-1), \, \boldsymbol{\cdot} \,},w,q[\overline{t}_r]) \mathrm{d} \mathrm{P}(w), \label{eq:id}
\end{align}
where $q[\overline{t}_r]$ is the history of the path $q = (q(0), q(1), \ldots, q(r-1))$ before time $t_r$ with transition times $\overline{t}_r = \{t_s: s=1,\ldots,r-1\}$ and $x(t_r) = q(r-1)$, and $\mathrm{P}(w)$ is the distribution function of covariates. Here, $\mathrm{d} \mathrm{P}_r(t_r \mid a_{q(r-1), \, \boldsymbol{\cdot} \,},w,q[\overline{t}_r])$ is the transition-specific intensity of moving from state $q(r-1)$ to $q(r)$ at time $t_r$,
\begin{equation} \label{eq:dens}
\begin{aligned}
\mathrm{d} \mathrm{P}_r(t_r \mid a_{q(r-1), \, \boldsymbol{\cdot} \,},w,q[\overline{t}_r]) &= \exp\left\{-\sum_{ k \in \underline{\mathcal{G}}_{q(r-1)}} \Lambda_{q(r-1),k}^{a_{q(r-1),k}}(t_r \mid w,q[\overline{t}_r])\right\} \\
&\quad \times \mathrm{d} \Lambda_{q(r-1),q(r)}^{a_{q(r-1),q(r)}}(t_r \mid w,q[\overline{t}_r]).
\end{aligned}
\end{equation}
\end{theorem}
The proof of this identifiability result is provided in Supplementary Material E.1. If all units are given the treatment $a$, then the potential cumulative incidences $F_j^a(t) = \mathrm{P}(j \in \overline{X}^a(t))$ of each event $j \in \mathcal{G}$ is identifiable by integrating event-specific hazards into cumulative incidences using the Kolmogorov forward equation. Based on the identification formula, we observe that the counterfactual cumulative incidence $F_j^{\overline{a}}(t)$ coincides with the potential cumulative incidence $F_j^a(t)$ if all components of the hypothetical treatment vector are identical, i.e., $a_{kg} = a$ for every $k,g\in\mathcal{E}$. That is, $F_j^{\overline{1}}(t) = F_j^1(t)$ and $F_j^{\overline{0}}(t) = F_j^0(t)$, where $\overline{1} = (1, \ldots, 1)$ and $\overline{0} = (0, \ldots, 0)$. Therefore, we can decompose the total effect $\tau_j(t) = F_j^1(t) - F_j^0(t)$ on event $j\in\mathcal{G}$ into a summation of individual effects via transitions.
\begin{remark}
To illustrate how to decompose the total effect into interventional effects, we consider an illness--death (I--D) model. The hypothetical treatment vector $\overline{a} = (a_{\text{OI}},a_{\text{OD}},a_{\text{ID}})$ has three components, corresponding to transitions from the initial state to illness, from the initial state to death, and from illness to death, respectively. The interventional direct effect on death is $F_{\text{D}}^{(0,1,1)}(t) - F_{\text{D}}^{\overline{0}}(t)$, representing the difference in counterfactual cumulative incidences of death resulting from manipulating the hazard of death while keeping the hazard of illness controlled. The interventional indirect effect on death is $F_{\text{D}}^{\overline{1}}(t) - F_{\text{D}}^{(0,1,1)}(t)$, representing the difference in counterfactual cumulative incidences of death due to altering the hazard of disease while controlling the hazard of illness. These two effects sum to the total effect $F_{\text{D}}^1(t) - F_{\text{D}}^0(t)$, which coincides with the natural effects under sequential ignorability and the separable effects under the dismissible components condition in semi-competing risks when there is no treatment-induced confounding \citep{deng2024direct, breum2024estimation}. Furthermore, the direct effect on death can be decomposed into an effect via the transition from the initial state to death $F_{\text{D}}^{(0,1,0)}(t) - F_{\text{D}}^{\overline{0}}(t)$, and an effect via the transition from illness to death $F_{\text{D}}^{(0,1,1)}(t) - F_{\text{D}}^{(0,1,0)}(t)$.
The path-specific effect along the path ${\text{O}}\to{\text{I}}\to{\text{D}}$ is $F_{\text{D}}^{(1,0,1)}(t) - F_{\text{D}}^{\overline{0}}(t)$, and the path-specific effect along the path ${\text{O}}\to{\text{D}}$ is $F_{\text{D}}^{(0,1,0)}(t) - F_{\text{D}}^{\overline{0}}(t)$. Compared with natural effects, identifying interventional effects does not require the cross-world independence assumption required for natural effects; compared with separable effects, identifying interventional effects does not require partial isolation that excludes interacting effect pathways of treatment-induced confounding across different events.
\end{remark}
\section{Estimation and Inference} \label{sec:estim}
\subsection{Efficient Influence Function}
In the real-world study, we observe $n$ independent and identically distributed (i.i.d.) copies of $O = (A, W, \overline{X}(t^*))$, denoted by $\mathcal{O} = \{ O_i = (A_i, W_i, \overline{X}_i(t^*)): i=1,\ldots,n \}$. Equivalently, the trajectory $\overline{X}(t^*)$ can be written as $(T_k,\Delta_k: k\in\mathcal{G})$. Unlike competing risks models, estimating hazard functions nonparametrically is difficult in multi-state models because the origins of transitions between arbitrary states are not aligned. As a compromise, parametric or semiparametric models are usually employed, restricting the dependence of hazards on the history. The set of unordered states in the transition history reflects the current status of an individual; therefore, it is reasonable to expect that the transition hazard obeys the extended Markov property. We assume that the hazard $\mathrm{d}\Lambda_{kg}^{a}(t \mid w, \overline{x}(t))$ at time $t$ depends only on $w$ and the unordered states in the path $\{q(s): s = 0, \dots, l(q), t_s<t\}$ rather than the past transition times $\{t_s: s=1,\ldots,l(q)\}$. Conditional on covariates at a fixed time (see Figure \ref{fig:model}), the Markovness condition implies that the counterfactual hazard of death may take one of four distinct values corresponding to the paths $q_1$, $q_2$, $q_3$, and $\{q_4,q_5\}$. We formalize the Markovness condition in the following statement.
\begin{assumption}[Extended Markovness] \label{cond1}
For every possible state history $\overline{x}(t)$, there exists a set of history $h(\overline{x}(t))$ (unordered states) such that $\mathrm{P} (\overline{X}(t) \in h(\overline{x}(t)) \mid W) > \epsilon$ for some constant $\epsilon>0$ and
\begin{align*}
\mathrm{d}\Lambda_{kg}^{a}(t \mid w, \overline{x}(t)) &= \mathrm{d}\Lambda_{kg}^{a}(t \mid w, h(\overline{x}(t))), \ (k,g)\in\mathcal{E}.
\end{align*}
\end{assumption}
Extended Markovness automatically holds in competing risks data because all transitions share the same at-risk set. The at-risk probability declines from 1 since the time origin. Extended Markovness ensures that there is sufficient data to flexibly estimate the hazards. It excludes the case where the second event in a path occurs just after the baseline. This assumption is testable \citep{grambsch1994proportional}, and can be relaxed by modeling sojourn times. Let $\mathbb{P}_n(\cdot)$ denote the empirical average over the sample.
By replacing the unknown hazards with their corresponding semiparametric estimators, we obtain the \textit{plug-in} (or \textit{regression}) estimator for the conditional counterfactual cumulative incidence $\widehat{F}_j^{\overline{a}}(t \mid w)$ according to the identification result in Theorem \ref{thm1}. Then, the population-level counterfactual cumulative incidence is estimated by $\widehat{F}_j^{\overline{a}}(t) = \mathbb{P}_n \{\widehat{F}_j^{\overline{a}}(t \mid W)\}$.
Plug-in estimators suffer from model misspecification. To increase estimation efficiency and robustness, we derive the efficient influence function (EIF) of the counterfactual cumulative incidence $F_j^{\overline{a}}(t)$, which is in the tangent space restricted by extended Markovness.
\begin{lemma} \label{thm2}
Under Assumptions \ref{ass:ign}--\ref{cond1}, the efficient influence function of the counterfactual cumulative incidence $F_j^{\overline{a}}(t)$ in the semiparametric model space restricted by extended Markovness is given by
\begin{equation} \label{eif1}
\begin{aligned}
\operatorname{EIF}\{F_j^{\overline{a}}(t)\} &= \sum_{q\in\mathcal{Q}_j} \int_0^t\int_0^{t_{l(q)}}\cdots\int_0^{t_2} \prod_{r= 1}^{l(q)} \mathrm{d} \mathrm{P}_{r}(t_{r} \mid a_{q(r-1),\, \boldsymbol{\cdot} \,},W,q[\overline{t}_{r}]) \\
&\quad \times \sum_{r = 1}^{l(q)} \Bigg[ \frac{\mathrm{d} \eta_{q(r-1),q(r)}^{a_{q(r-1),q(r)}}(t_r \mid W,q[\overline{t}_r])}{\mathrm{d} \Lambda_{q(r-1),q(r)}^{a_{q(r-1),q(r)}}(t_r \mid W,q[\overline{t}_r]) } - \sum_{k\in\underline{\mathcal{G}}_{q(r-1)}} \eta_{q(r-1),k}^{a_{q(r-1),k}}(t_r \mid W,q[\overline{t}_r]) \Bigg] \\
&\quad + F_j^{\overline{a}}(t \mid W) - F_j^{\overline{a}}(t),
\end{aligned}
\end{equation}
where
\begin{equation} \label{eq:mart}
\begin{aligned}
\eta_{q(r-1),k}^{a_{q(r-1),k}}(t \mid w,\overline{x}(t)) & = \frac{\mathds{1}(A=a_{q(r-1),k}, W=w)}{\mathrm{P} (A=a_{q(r-1),k}\mid W=w)} \\
&\quad \times \int_0^t \frac{\mathrm{d} M_{q(r-1),k}(s;a_{q(r-1),k}, w, \overline{x}(s))}{\mathrm{P} \{T_{k}\wedge T_{\text{c}}\geq s, \overline{X}(s) \in h(\overline{x}(s)) \mid A=a_{q(r-1),k}, W=w\}},
\end{aligned}
\end{equation}
and $M_{q(r-1),k}(t;a_{q(r-1),k}, w, \overline{x}(t))$ is the martingale process associated with the transition from $q(r-1)$ to $k$,
\[
M_{q(r-1),k}(t;a, w, \overline{x}(t)) = \int_0^t \mathds{1}\{\overline{X}(s)\in h(\overline{x}(s))\} \left\{\mathrm{d}N_{k}(s) - Y_k(s) \mathrm{d}\Lambda_{q(r-1),k}^{a_{q(r-1),k}}(s\mid w,\overline{x}(s)) \right\},
\]
with the event-counting process $N_{k}(t) = \Delta_k \mathds{1} \{T_k \leq t\}$ and at-risk process $Y_k(t) = \mathds{1} \{T_k \wedge T_{\text{c}} \geq t\}$.
\end{lemma}
The proof of Lemma \ref{thm2} is given in Supplementary Material E.2. The EIF is a function of the observed data $O$, characterizing the influence of a single observation on the target estimand. In the context of competing risks, the function $\eta_{0,k}^{a}(t \mid w)$ is exactly the EIF of the potential cause-specific hazard of event $k$ under treatment $a$. Compared with the identification formula based on $F_j^{\overline{a}}(t|W)$, an additional term in the EIF
\begin{align*}
\varphi_j^{\overline{a}}(t; O) &= \sum_{q\in\mathcal{Q}_j} \int_0^t\int_0^{t_{l(q)}}\cdots\int_0^{t_2} \prod_{r= 1}^{l(q)} \mathrm{d} \mathrm{P}_{r}(t_{r} \mid a_{q(r-1),\, \boldsymbol{\cdot} \,},W,q[\overline{t}_{r}]) \\
&\quad \times \sum_{r = 1}^{l(q)} \Bigg[ \frac{\mathrm{d} \eta_{q(r-1),q(r)}^{a_{q(r-1),q(r)}}(t_r \mid W,q[\overline{t}_r])}{\mathrm{d} \Lambda_{q(r-1),q(r)}^{a_{q(r-1),q(r)}}(t_r \mid W,q[\overline{t}_r]) } - \sum_{k\in\underline{\mathcal{G}}_{q(r-1)}} \eta_{q(r-1),k}^{a_{q(r-1),k}}(t_r \mid W,q[\overline{t}_r]) \Bigg]
\end{align*}
whose expectation is zero, incorporates information in every single observation to reduce the bias of the regression estimator.
\begin{remark}
The $\mathrm{d} \Lambda_{q(r-1),q(r)}^{a_{q(r-1),q(r)}}(t_r \mid W,q[\overline{t}_r])$ that appears in the denominator of $\operatorname{EIF}\{F_j^{\overline{a}}(t)\}$ cancels out with $\mathrm{P}_{r}(t_{r} \mid a_{q(r-1),\, \boldsymbol{\cdot} \,},W,q[\overline{t}_{r}])$. This is important for fitting the EIF if we employ semiparametric models for the hazard functions. For example, in the Cox model, the Breslow estimator consistently estimates the cumulative hazard but not the pointwise hazard \citep{breslow1972disussion}. Consistent estimation of cumulative hazards is sufficient to fit the EIF. If the model space is not restricted by extended Markovness, the EIF will involve pointwise hazards (see Supplementary E.5), and consistently estimating pointwise hazards is required in this case, which is much more inefficient.
\end{remark}
\subsection{Estimation Based on the Efficient Influence Function}
The efficient influence function involves transition hazards, censoring hazards, and the propensity score. We fit the transition hazards using semiparametric models. For example, we may use the Cox proportional hazards model with the occurrence statuses of other events as time-varying covariates. In particular, the transition hazard from $k$ to $g$ can be modeled as
\[
\mathrm{d} \Lambda_{kg}^a (t \mid w, \overline{x}(t)) = \mathrm{d} \Lambda_{kg}^a(t) \exp\left\{\beta_{kg}^{{ \mathrm{\scriptscriptstyle T} }}w + \sum_{h\in\mathcal{G}\backslash\{\text{o},k,g\}}\gamma_{k,hg}N_h(t)\right\},
\]
where $\Lambda_{kg}^a(t)$ is the unspecified baseline hazard, $\beta_{kg}$ is the coefficients for baseline covariates, $\gamma_{k, hg}$ is the coefficients for event history, and $N_h(t)$ is the status (counting process) indicating whether event $h$ has occurred at time $t$. If $x(t) \neq k$, we set $\mathrm{d} \Lambda_{kg}^a (t \mid w, \overline{x}(t)) = 0$. If an event $h\in\mathcal{G}$ does not appear in any possible trajectory including $k\to g$, the coefficient $\gamma_{hg}$ is set to zero. Since treatment-induced confounding is integrated out by the definition of a randomized intervention, we do not need to model its intensity. Model parameters can be estimated by fitting Cox models separately for each event $g\in\mathcal{G}\backslash\{\text{o}\}$ using nonparametric maximum likelihood \citep{breslow1972disussion, breslow1974covariance, zeng2007maximum}, where $W$ is treated as time-invariant covariates and $\{N_h(\cdot): h\in\mathcal{G}\backslash\{\text{o},k,g\}\}$ are treated as time-varying covariates. The censoring hazards $\Lambda_{k\text{c}}^a(t\mid w,\overline{x}(t))$ are fitted using similar models. The propensity score can be fitted using generalized linear models such as logistic regression,
\[
\mathrm{P}(A=1 \mid W=w) = 1/\{1 + \exp(-\alpha^{{ \mathrm{\scriptscriptstyle T} }}w)\},
\]
where the parameters are estimated using maximum likelihood.
Using the fitted transition hazards $\widehat\Lambda_{kg}^a(t\mid w,\overline{x}(t))$, we fit the transition density $\mathrm{d} \widehat{\mathrm{P}}_r(t_r \mid a_{q(r-1),\cdot},W,q[\overline{t}_r])$ along each path at every time point for each unit.
In particular, for each path $q = (q(0),q(1),\ldots,q(r-1),q(r))$ with transition times $\overline{t}_r = (t_1,\ldots,t_{r-1},t_r)$,
\begin{align*}
\mathrm{d} \widehat{\mathrm{P}}_r(t_r \mid a_{q(r-1), \, \boldsymbol{\cdot} \,},W,q[\overline{t}_r]) &= \exp\left\{-\sum_{ k \in \underline{\mathcal{G}}_{q(r-1)}} \int_{t_{r-1}}^{t_r} \mathrm{d} \widehat\Lambda_{q(r-1),k}^{a_{q(r-1),k}}(t_r \mid W,q[\overline{t}_r])\right\} \\
&\quad \times \mathrm{d} \widehat\Lambda_{q(r-1),q(r)}^{a_{q(r-1),q(r)}}(t_r \mid W,q[\overline{t}_r]).
\end{align*}
By plugging the estimated densities into Equation \eqref{eq:dens}, we estimate the counterfactual cumulative incidence $\widehat{F}_j^{\overline{a}}(t \mid W;q)$ along each path $q\in\mathcal{Q}_j$. The plug-in estimator for $F_j^{\overline{a}}(t \mid W)$ is obtained by $\widehat{F}_j^{\overline{a}}(t \mid W) = \sum_{q\in\mathcal{Q}_j}\widehat{F}_j^{\overline{a}}(t \mid W;q)$, and the plug-in estimator for $F_j^{\overline{a}}(t)$ is $\widehat{F}_j^{\overline{a}}(t) = \mathbb{P}_n\{\widehat{F}_j^{\overline{a}}(t \mid W)\}$.
Next, we augment the multi-state model by regarding censoring $\text{c}$ as an absorbing state, so that $\mathcal{G}^* = \mathcal{G}\cup\{\text{c}\}$ and $\mathcal{Q}_g^* = \mathcal{Q}_g \cup \{(k,\text{c}): k\in\mathcal{G}\}$. Using the fitted transition hazards $\widehat\Lambda_{kg}^a(t\mid w,\overline{x}(t))$ and censoring hazards $\widehat\Lambda_{k\text{c}}^a(t\mid w,\overline{x}(t))$, we can estimate the counterfactual cumulative incidence $\widehat{F}_{j}^{*a}(t\mid w;q)$ for each $j\in\mathcal{G}$ when all treatment components are set at $a\in\{0,1\}$. Using the fitted propensity score $\widehat{\mathrm{P}}(A=a\mid W=w)$ we fit the weighted integral martingale $\widehat\eta_{q(r-1),k}^{a_{q(r-1),k}}(t_r \mid W,q[\overline{t}_r])$ based on Equation \eqref{eq:mart},
\begin{align*}
&\widehat\eta_{q(r-1),k}^{a_{q(r-1),k}}(t_r \mid W,q[\overline{t}_r]) = \frac{1}{\widehat{\mathrm{P}} (A=a_{q(r-1),k}\mid W)} \\
&\qquad \times \int_{t_{r-1}}^{t_r} \frac{\mathrm{d} N_{k}(s;a_{q(r-1),k},W,q[\overline{t}_r]) - Y_k(s;a_{q(r-1),k},W,q[\overline{t}_r]) \mathrm{d}\widehat\Lambda_{q(r-1),k}^{a_{q(r-1),k}}(s\mid W,q[\overline{t}_r])}{\widehat{\mathrm{P}}\{T_k \wedge T_{\text{c}} \geq s, \overline{X}(s) \in h(q[\overline{t}_r]) \mid A=a_{q(r-1),k}, W\}},
\end{align*}
where
\begin{align*}
&\quad ~ \widehat{\mathrm{P}}\{T_k \wedge T_{\text{c}} \geq s, \overline{X}(s) \in h(q[\overline{t}_r]) \mid A=a_{q(r-1),k}, W\} \\
&=
\sum_{q': h(q')=h(q)} \bigg\{ \widehat{F}_{q(r-1)}^{*a_{q(r-1),k}}(s\mid W;q') - \sum_{g\in\underline{\mathcal{G}}^*_{q(r-1)}} \widehat{F}_g^{*a_{q(r-1),k}}(s\mid W;(q',g)) \bigg\}
\end{align*}
is the estimated probability of remaining in state $q(r-1)$ at time $s$ with a history compatible with $q[\overline{t}_r]$ under treatment $A=a_{q(r-1),k}$ in the presence of censoring.
Let $\widehat\varphi_j^{\overline{a}}(t;O)$ denote the fitted debiasing term obtained by substituting the unknown components with their estimated counterparts.
We update the plug-in estimator $\mathbb{P}_n \{\widehat{F}_j^{\overline{a}}(t \mid W)\}$ by adding the empirical average of the debiasing term,
\begin{equation}\label{cond_2_est_eq}
\begin{aligned}
\widetilde{F}_j^{\overline{a}}(t) = \mathbb{P}_n \left[ \widehat{F}_j^{\overline{a}}(t \mid W) + \widehat\varphi_j^{\overline{a}}(t; O) \right].
\end{aligned}
\end{equation}
The estimator $\widetilde{F}_j^{\overline{a}}(t)$ is known as a \textit{one-step} estimator because it gives a one-step update for the initial plug-in estimator. This EIF-based estimator corrects the first-order bias of the plug-in estimator. This EIF-based estimator solves the estimating equation of the fitted EIF,
\[
\mathbb{P}_n \left[\widehat{\operatorname{EIF}}\{F_j^{\overline{a}}(t)\}\right] := \mathbb{P}_n \left[ \widehat{F}_j^{\overline{a}}(t \mid W) + \widehat\varphi_j^{\overline{a}}(t; O) - \widetilde{F}_j^{\overline{a}}(t)\right] = 0.
\]
The standard error of the EIF-based estimator is given by
\[
\mathrm{s.e.}\{\widetilde{F}_j^{\overline{a}}(t)\} = \left(n^{-1} \mathbb{P}_n \left[\widehat{\operatorname{EIF}}\{F_j^{\overline{a}}(t)\}^2\right]\right)^{1/2}.
\]
\subsection{Asymptotic Properties} \label{sec:asy_pro}
The one-step estimator typically exhibits higher efficiency and robustness to model misspecification, making it more attractive than plug-in methods.
A key advantage of our EIF-based estimator is \textit{multiple robustness}.
\begin{theorem}[Multiple robustness] \label{thm_multi_robust}
Suppose that the assumptions in Lemma \ref{thm2} hold. The EIF-based estimator $\widetilde{F}_j^{\overline{a}}(t)$ is consistent if, for each $q\in\mathcal{Q}_j$, at least one of the following conditions holds:
\begin{enumerate}[(1)]
\item All the transition hazards from the states on the path $q$ to their adjacent following states $\{\Lambda_{q(r - 1), k}^a(\cdot): r = 1, \ldots, l(q), k \in \underline{\mathcal{G}}_{q(r - 1)}\}$ are correctly specified.
\item The propensity score and censoring hazard are correctly specified, while all but one transition hazard in $\{\Lambda_{q(r - 1), k}^a(\cdot): r = 1, \ldots, l(q), k \in \underline{\mathcal{G}}_{q(r - 1)}\}$ are correctly specified.
\end{enumerate}
\end{theorem}
Another appealing property of the EIF-based estimator is \textit{semiparametric efficiency}.
\begin{theorem}[Semiparametric efficiency] \label{thm_asy_normal}
Suppose that the assumptions in Lemma \ref{thm1} and the additional regularity conditions listed in Supplementary Material A hold. Then
\begin{align*}
\sqrt{n} \{\widetilde{F}_j^{\overline{a}}(t) - F_j^{\overline{a}}(t)\} \rightsquigarrow \mathcal{N}\left(0, \ \mathrm{E}[\operatorname{EIF}\{F_j^{\overline{a}}(t)\}]^2\right)
\end{align*}
for $t \in [0,t^*]$. The asymptotic variance attains the semiparametric efficiency bound.
\end{theorem}
The regularity conditions essentially require that the fitted models are not too complex (Donsker condition). They should be correctly specified, and the convergence rate is not too slow. Note that the EIF is Lipschitz continuous with respect to the hazard model and propensity score model. The Cox model with a bounded baseline hazard and time-varying covariates of bounded variation is Donsker \citep{zeng2007maximum}, and the logistic regression model for the propensity score is Donsker. Therefore, the fitted EIF is Donsker. If these working models are correctly specified, then the fitted models converge at the $O_p(n^{-1/2})$ rate. As a result, all regularity conditions for Theorem \ref{thm_asy_normal} hold. Theorem \ref{thm_asy_normal} indicates that the EIF-based estimator $\widetilde{F}_j^{\overline{a}}(t)$ is regular and asymptotically linear (RAL), with an influence function $\operatorname{EIF}\{F_j^{\overline{a}}(t)\}$. The pointwise confidence interval (CI) of $F_j^{\overline{a}}(t)$ can be constructed by the normal approximation.
The proofs of Theorems \ref{thm_multi_robust} and \ref{thm_asy_normal} are given in Supplementary Material E.3 and E.4.
\subsection{Treatment Effects Comparing Sequential Interventions}
The primary goal of causal inference for multi-state data is to answer the following three questions: identifying transition-specific effects to determine the transitions delivering treatment effects, identifying event-specific effects to determine the target of treatment, and identifying path-specific effects to determine the causal pathway. Figure \ref{fig:tes} illustrates three possible ways for envisioning hypothetical interventions. In (a), we let the one-step transition from O to E be treated while other transitions are controlled, and then we can find the treatment effect on the transition (O, E). In (b), we let all one-step transitions to M be treated while other transitions are controlled, and then we can find the treatment effect on M. In (c), we treat the transitions along the path O$\to$E$\to$D while keeping the other transitions controlled, and then we can estimate the treatment effect delivered through this path.
\begin{figure}[!tb]
\centering
\includegraphics[width=0.85\textwidth]{Figure3.pdf}
\caption{Three versions of interventions. Solid arrows represent interventions with hazards under treatment, and dashed lines represent interventions with hazards under control. (a) The treatment effect on the transition O$\to$E. (b) The treatment effect on the event M. (c) The treatment effect on the path O$\to$E$\to$D.} \label{fig:tes}
\end{figure}
Generally, for a target state $j \in \mathcal{G}$, the treatment effect comparing two hypothetical treatment vectors $\overline{a}$ and $\overline{a}^*$ is defined as
\begin{equation}
\tau_j(t;\overline{a},\overline{a}^*) := F_j^{\overline{a}}(t) - F_j^{\overline{a}^*}(t).
\end{equation}
Here, $\overline{a}$ represents the treatment sequence of interest, and $\overline{a}^*$ serves as the baseline treatment \citep{robins2004optimal}. This estimand can be used to quantify whether a treatment effect exists on a specific path.
Specially, we may set all components in $\overline{a}^*$ at 0, then $\tau_j(t;\overline{a},\overline{a}^*)$ represents the effect of partially active treatment. We may also set all components in $\overline{a}^*$ to 1, then $\tau_j(t;\overline{a},\overline{a}^*)$ represents the effect of abstaining from treatment.
Let $\widetilde\tau_j(t;\overline{a},\overline{a}^*) = \widetilde{F}_j^{\overline{a}}(t) - \widetilde{F}_j^{\overline{a}^*}(t)$ be the EIF-based one-step estimator for the treatment effect.
Inference for the treatment effect is straightforward as a result of Theorem \ref{thm_asy_normal} and the linearity of efficient influence functions.
\begin{corollary}\label{cor_inference_path}
Suppose that the conditions in Theorem \ref{thm_asy_normal} hold, then for $t \in [0,t^*]$,
\[
\sqrt{n} \{\widetilde{\tau}_j(t;\overline{a},\overline{a}^*) - \tau_j(t;\overline{a},\overline{a}^*)\} \rightsquigarrow \mathcal{N}\left(0, \ \mathrm{E}[\operatorname{EIF}\{F_j^{\overline{a}}(t)\}-\operatorname{EIF}\{F_j^{\overline{a}^*}(t)\}]^2\right).
\]
\end{corollary}
For hypothesis testing of the treatment effect, we construct a test statistic based on the event-specific restricted mean survival time lost (RMSTL) by the end of the study as
\begin{equation} \label{teststat}
T_j(\overline{a},\overline{a}^*) = \int_0^{t^*} \{\widetilde{F}_j^{\overline{a}}(s) - \widetilde{F}_j^{\overline{a}^*}(s)\} \mathrm{d}s.
\end{equation}
Under the null hypothesis that there is no effect, the test statistic $T_j(\overline{a},\overline{a}^*)$ should converge to 0 in probability. Its asymptotic variance can be derived from the EIF.
\begin{corollary}\label{cor_inference_test}
Suppose that the conditions in Theorem \ref{thm_asy_normal} hold.
Under the null hypothesis $H_0: \tau_j(t;\overline{a},\overline{a}^*) = 0$ for all $t\in[0,t^*]$,
\[
\sqrt{n} T_j(\overline{a},\overline{a}^*) \rightsquigarrow \mathcal{N}\left(0, \ \sigma_j^2(\overline{a},\overline{a}^*)\right),
\]
where
\begin{align*}
\sigma_j^2(\overline{a},\overline{a}^*) &= \operatorname{var}\left(\int_0^{t^*} [\operatorname{EIF}\{F_j^{\overline{a}}(s)\} - \operatorname{EIF}\{F_j^{\overline{a}^*}(s)\}] \mathrm{d}s\right).
\end{align*}
\end{corollary}
In Supplementary Material C, we conduct simulation studies to assess the finite-sample performance of the EIF-based estimator in various settings, confirming its consistency and multiple robustness. The coverage rate of confidence intervals based on the estimated EIF is close to the nominal level. We also find that using an estimated propensity score can yield a slightly lower finite-sample variance than using the true propensity score.
\section{Application to the LEADER Trial} \label{sec:real_data}
\subsection{Data Description}
Vascular events are leading causes of mortality among patients with type 2 diabetes \citep{stratton2000association, advance2008intensive}. In clinical studies, these events are typically classified as expanded major adverse cardiovascular events (EMACE) and microvascular events (MVE), while other types of vascular events are relatively uncommon \citep{emerging2010diabetes, forbes2013mechanisms}. EMACE includes myocardial infarction, stroke, unstable angina pectoris, heart failure, coronary revascularization, and cardiovascular death, while MVE includes nephropathy and retinopathy. These two categories of non-fatal vascular events can occur alternatively, with one potentially leading to or influencing the development of the other before ultimately resulting in death.
The LEADER (Liraglutide Effect and Action in Diabetes: Evaluation of Cardiovascular Outcome Results) Trial was a large, multi-center, double-blind, randomized controlled trial designed to investigate the long-term effects of liraglutide on cardiovascular outcomes in patients with type 2 diabetes at high cardiovascular risk \citep{marso2013design, marso2016liraglutide}. Conducted across 410 clinical research sites in 32 countries as part of a global phase 3a program, the trial enrolled 9,340 participants, who were randomized in equal proportions to receive either liraglutide (a glucagon-like peptide-1 receptor agonist) or placebo. The LEADER trial offers a valuable opportunity to investigate the impact of liraglutide on the progression of vascular events. Standard Kaplan--Meier regression (associated with log-rank test) and Cox regression (associated with Wald or score test) demonstrated that liraglutide significantly reduced the risks of EMACE and MVE. Analyses of individual events further revealed significant reductions in cardiovascular death and nephropathy.
However, these results do not distinguish between the direct effect on a specific event and the indirect effects mediated through other events. Taking death as an example, an individual can either experience death without any vascular events or experience both vascular events and death. To gain a clearer understanding of how treatment improves survival, it is necessary to examine all potential pathways from treatment initiation through subsequent clinical events. Elucidating these causal mechanisms can enhance our understanding of liraglutide’s efficacy and safety. Furthermore, identifying the specific pathways through which liraglutide exerts its effects may foster precision medicine by enabling targeted treatment for individuals at high risk of specific disease pathways. Figure \ref{fig:model} illustrates the five potential paths from baseline to death, where the first occurrences of EMACE and MVE are regarded as intermediate states. The interventional effects framework is well-suited to studying transition-, event-, and path-specific effects in this example.
In the LEADER Trial, 4,668 individuals were assigned to liraglutide and 4,672 to placebo. We adjust for seven baseline covariates: age (continuous), sex (binary: male or female), BMI (binary: normal or high), HbA1c (binary: normal or high), diabetes duration (continuous), cardiovascular risk (binary: high or medium), and insulin naive (binary: yes or no). Although the trial was completely randomized and baseline covariates were balanced, we estimate the propensity score using logistic regression that adjusts for these covariates, which may improve finite-sample efficiency \citep{hirano2003efficient}. There were 4 missing values for age, 9 for BMI, and 19 for diabetes duration. The missingness rate is low, and the hazard estimation is insensitive to imputation methods in pilot analyses. In the formal analysis below, we impute missing values using 5-nearest-neighbor imputation. The LEADER Trial is completely randomized, so ignorability holds. Censoring is administrative, so random censoring holds. The choice of $t^*=60$ ensures positivity. Since there is only one liraglutide formulation, consistency holds.
\subsection{Effect of Liraglutide on Vascular Events and Death}
We consider EMACE (E) and MVE (M) as intermediate events and death (D) as the terminal event. We observe individuals along every path to death in each treatment group. Note that EMACE includes both direct cardiovascular death (CVD) and non-fatal events. To distinguish causes of death, we model the hazards of CVD and other causes of death separately. Since cardiovascular death is part of EMACE by definition, the intervention in the transition to direct CVD is always identical to the intervention in the transition (O,E). Cause-specific hazards for non-fatal EMACE, MVE, direct CVD, and other-cause death were estimated separately in each treatment group using Cox models, adjusting for baseline covariates, with occurrences of other non-fatal events treated as time-varying covariates. The transition hazard from direct CVD to D is set at infinity, meaning that an individual with direct CVD in E transits to D immediately. The estimated coefficients are listed in Supplementary Material D.1.
To determine the direct effect on each event, we estimate interventional effects by switching the treatment on the transitions into the event of interest. For E, we treat transitions into E ($a_{jk}=1$ for $k=\text{E}$) while controlling all other transitions ($a_{jk}=0$ for $k\neq\text{E}$). The difference in counterfactual cumulative incidences of EMACE (including non-fatal EMACE and direct CVD) between this intervention and the all-placebo scenario represents the direct effect on EMACE. Similarly, the difference in counterfactual cumulative incidences of MVE between the world with interventions in the transitions to MVE and the all-placebo scenario represents the direct effect on MVE. The difference in counterfactual cumulative incidences of death (including direct CVD and other-cause death) between the world with interventions in the transitions to other-cause death and the all-placebo scenario represents the direct effect on other-cause death. Figure \ref{fig:ate1} displays the direct effects on these three categories of events. The direct effects on EMACE ($P=0.0013$) and MVE ($P=0.0040$) are significant, while the direct effect on other-cause death ($P=0.8638$) is insignificant.
\begin{figure}[!tb]
\centering
\includegraphics[width=0.96\textwidth]{Figure4.pdf}
\caption{The estimated direct treatment effect on EMACE, MVE, and other-cause death by intervening in the hazards of the event of interest while controlling for other hazards at placebo, with 95\% confidence intervals.} \label{fig:ate1}
\end{figure}
\subsection{Transition-Specific and Path-Specific Effects}
Now we examine the treatment effect on death through each one-step transition; this helps us understand how liraglutide reduces the risks of EMACE and MVE. We first assume all transitions are initially under placebo. We intervene in a single transition and compare the resulting counterfactual cumulative incidence with that in the all-placebo world. Liraglutide has a significant effect on the transition from O to E ($P=0.0152$), whereas the effects on all other transitions are not significant. The reduced risk of EMACE through the transition from baseline to EMACE contributes to a lower risk of EMACE-induced death.
Figure \ref{fig:trans1} shows the restricted mean survival time gained with liraglutide via one-step transitions. Solid lines represent a reduction in risk, while dashed lines represent an increase in risk.
\begin{figure}[!tb]
\centering
\includegraphics[width=0.85\textwidth]{Figure5.pdf}
\caption{The restricted mean survival time gained by liraglutide through one-step transitions, (a) taking all-placebo as the reference, and (b) taking all-liraglutide as the reference. Solid lines represent liraglutide increases survival, and dashed lines represent liraglutide reduces survival.} \label{fig:trans1}
\end{figure}
By intervening in each path from baseline to D, we find that the restricted mean survival time gained by liraglutide is $-$0.0751 (s.e. 0.1166) months along the path O$\to$D, 0.3962 (s.e. 0.1576) months along the path O$\to$E$\to$D, $-$0.0202 (s.e. 0.0416) months along the path O$\to$M$\to$D, 0.0787 (s.e. 0.0773) months along the path O$\to$M$\to$E$\to$D, and 0.3198 (s.e. 0.1468) months along the path O$\to$E$\to$M$\to$D. The estimated treatment effect curves under each path-specific intervention are presented in Supplementary Material D.1. The effects along the paths O$\to$E and O$\to$E$\to$M$\to$D count for the largest proportions, and they are significant at the 0.05 level. The result suggests that treatment in individuals with a high risk of EMACE is crucial to reducing mortality.
In Supplementary Material D.2, we present the total effects on each event. In Supplementary Material D.3, we infer dynamic treatment policies to reduce event risks under sequential ignorability. In Supplementary Material D.4, we conduct a sensitivity analysis by allowing for a parametrically distributed unmeasured confounder using frailty modeling, such that randomized intervention is achieved additionally conditional on this unmeasured confounder \citep{rotolo2016incorporation, li2020shared, gu2024maximum}. We assume the randomized intervention is performed conditional on both observed covariates and frailty; technical details are provided in Supplementary Material B. The results are similar. The small size of frailty indicates that the influence of unmeasured confounding is weak, even if it exists. In Supplementary Material D.5, we conduct a secondary analysis by considering the first occurrence of nine mutually exclusive individual events (myocardial infarction, stroke, unstable angina pectoris, heart failure, coronary revascularization, nephropathy, retinopathy, cardiovascular death, non-cardiovascular death). We find that liraglutide has significant direct effects on cardiovascular death and nephropathy.
The analyses above provide new knowledge on how liraglutide reduces mortality over existing studies. There are three sources of death: direct cardiovascular death (fatal EMACE at first occurrence), non-cardiovascular death, and death following non-fatal EMACE or MVE. Firstly, liraglutide reduces the risk of direct cardiovascular death. Secondly, liraglutide does not affect non-cardiovascular death. Thirdly, liraglutide reduces the risk of non-fatal EMACE and MVE, which in turn lowers the risk of death. Collectively, these findings suggest that liraglutide reduces the risk of all-cause death. These findings are supported by a sensitivity analysis using frailty modeling, which accounts for potential unmeasured confounding.
\section{Discussion} \label{sec:discus}
The key assumption for envisioning the counterfactual cumulative incidence is the randomized intervention in event-counting processes according to hazards. This intervention is counterfactually achievable regardless of unmeasured confounding, allowing inference on path-specific effects. Such a treatment falls within the hypothetical strategy of the International Conference on Harmonization guideline E9 (R1) addendum. If sequential ignorability holds, the randomized interventional effects can be interpreted as natural effects. Under sequential ignorability, the counterfactual cumulative incidences will be identical to the real cumulative incidences in future experiments where the corresponding treatment sequence is applied. Therefore, our framework can potentially inform the development of optimal dynamic treatment policies. Inferring dynamic treatment policies relies on sequential ignorability, that is, the absence of treatment-induced confounding \citep{robins2004optimal}. In future experimental or observational studies, different treatments could be assigned at each stage, providing opportunities to test for unmeasured confounding and evaluate the performance of the inferred treatment policy. Our estimation framework can be extended to accommodate sequential treatments.
Estimating the cumulative incidence functions for multi-state data is much more complex than for competing risks data. First, estimating transition hazards is complicated because they depend on the entire transition history, which comprises multiple ordered events. For competing risks data, all individuals share the same at-risk process; however, the time origins of being at risk in intermediate states for multi-state data are not aligned, so working assumptions on transition hazards should be imposed to ensure robust estimation. Second, censoring renders the complete transition paths unobserved. There may be a shared transition between two paths from the baseline to the terminal state, but the intensities along this shared transition cannot be separated. If an individual is in an intermediate state at the time of censoring, we cannot predict how the individual will transition after censoring. Transition hazards can only be identified conditionally on the observed history up to the censoring and event time, leaving subsequent transitions unconditioned on.
Several promising avenues for future research remain. First, alternative interventions may be considered when time-varying confounders are present. An intuitive modification is intervening in event-counting processes conditional on time-varying covariates. Since such interventions may alter the covariate distribution, it becomes necessary to model the intensity of the time-varying covariates to derive counterfactual cumulative incidences.
Second, as an alternative to the randomized interventional effects framework, the separable effects framework assumes that the initial treatment can be decomposed into components, each exerting a direct effect on a specific event or a pathway \citep{stensrud2021generalized, stensrud2022separable, robins2022interventionist, chen2025definition}. Separable effects place the targets of treatment components on events rather than transitions. Identifying clinically meaningful treatment components with isolated effects is not always feasible, limiting the practical interpretation of separable effects.
Third, interventional effects are not sharp \citep{miles2022on}. Even if a null effect is obtained, it does not necessarily imply that there is no effect, because time-varying confounding may offset the genuine effect. However, in the presence of time-varying confounding, other mediation frameworks suffer from their own limitations, either well-definedness or identifiability, so interventional effects can at least recover useful knowledge that motivates future investigation on biological mechanisms.
Lastly, alternative working models are possible. Nonparametric working models, such as kernel smoothing, can estimate the hazards conditional on prior event times without Assumption \ref{cond1}. The counterfactual cumulative incidence may also be estimated using other EIF-based methods, such as the targeted minimum-loss-based estimation \citep{wang2025targeted}. However, the efficient estimator of CIF may not be solved recursively, as the adjusted model depends on the specific time point in the target CIF.
\section{Data Availability Statement}\label{data-availability-statement}
The data that support the findings of this study are available from Novo Nordisk A/S. Restrictions apply to the availability of these data, which were used under license for this study.
\phantomsection\label{supplementary-material}
\bigskip
\begin{center}
{\large\bf SUPPLEMENTARY MATERIAL}
\end{center}
\begin{description}
\item[Supplementary material]
The online Supplementary Material includes (A) regularity conditions, (B) details of sensitivity analysis, (C) simulation studies, (D) additional data analysis results, and (E) proofs of theoretical results. (pdf)
\item[R code for simulation]
Code for simulation. (.R)
\end{description}
\begin{thebibliography}{xx}
\harvarditem{{ADVANCE Collaborative Group}}{2008}{advance2008intensive}
{ADVANCE Collaborative Group} \harvardyearleft 2008\harvardyearright ,
`Intensive blood glucose control and vascular outcomes in patients with type
2 diabetes', {\em New England Journal of Medicine} {\bf 358}(24),~2560--2572.
\harvarditem[Andersen et~al.]{Andersen, Geskus, de~Witte \harvardand\
Putter}{2012}{andersen2012competing}
Andersen, P.~K., Geskus, R.~B., de~Witte, T. \harvardand\ Putter, H.
\harvardyearleft 2012\harvardyearright , `Competing risks in epidemiology:
possibilities and pitfalls', {\em International Journal of Epidemiology} {\bf
41}(3),~861--870.
\harvarditem{Andersen \harvardand\ Keiding}{2002}{andersen2002multi}
Andersen, P.~K. \harvardand\ Keiding, N. \harvardyearleft
2002\harvardyearright , `Multi-state models for event history analysis', {\em
Statistical Methods in Medical Research} {\bf 11}(2),~91--115.
\harvarditem[Bickel et~al.]{Bickel, Klaassen, Bickel, Ritov, Klaassen, Wellner
\harvardand\ Ritov}{1993}{bickel1993efficient}
Bickel, P.~J., Klaassen, C.~A., Bickel, P.~J., Ritov, Y., Klaassen, J.,
Wellner, J.~A. \harvardand\ Ritov, Y. \harvardyearleft 1993\harvardyearright
, {\em Efficient and adaptive estimation for semiparametric models}, Vol.~4,
Springer.
\harvarditem{Breslow}{1972}{breslow1972disussion}
Breslow, N. \harvardyearleft 1972\harvardyearright , `Disussion of regression
models and life-tables by {D. R. Cox}', {\em Journal of the Royal Statistical
Society Series B: Statistical Methodology} {\bf 34},~216--217.
\harvarditem{Breslow}{1974}{breslow1974covariance}
Breslow, N. \harvardyearleft 1974\harvardyearright , `Covariance analysis of
censored survival data', {\em Biometrics} {\bf 30}(1),~89--99.
\harvarditem{Breslow}{1975}{breslow1975analysis}
Breslow, N.~E. \harvardyearleft 1975\harvardyearright , `Analysis of survival
data under the proportional hazards model', {\em International Statistical
Review/Revue Internationale de Statistique} pp.~45--57.
\harvarditem[Breum et~al.]{Breum, Munch, Gerds \harvardand\
Martinussen}{2024}{breum2024estimation}
Breum, M.~S., Munch, A., Gerds, T.~A. \harvardand\ Martinussen, T.
\harvardyearleft 2024\harvardyearright , `Estimation of separable direct and
indirect effects in a continuous-time illness--death model', {\em Lifetime
Data Analysis} {\bf 30}(1),~143--180.
\harvarditem{Chen}{2012}{chen2012maximum}
Chen, Y.-H. \harvardyearleft 2012\harvardyearright , `Maximum likelihood
analysis of semicompeting risks data with semiparametric regression models',
{\em Lifetime Data Analysis} {\bf 18},~36--57.
\harvarditem{Chen \harvardand\ Lin}{2025}{chen2025definition}
Chen, Y.-L. \harvardand\ Lin, S.-H. \harvardyearleft 2025\harvardyearright ,
`Definition and interpretation of separable path-specific effects with
multiple ordered mediators', {\em Epidemiology} {\bf 36}(5),~677--685.
\harvarditem[Deng et~al.]{Deng, Wang, Zhang \harvardand\
Zhan}{2026}{deng2026randomized}
Deng, Y., Wang, R., Zhang, T. \harvardand\ Zhan, X. \harvardyearleft
2026\harvardyearright , `Randomized interventional effects in semicompeting
risks, with application to a hematopoietic cell transplantation study', {\em
Statistics in Medicine} {\bf 45}(13-14),~e70628.
\harvarditem[Deng et~al.]{Deng, Wang \harvardand\ Zhou}{2024}{deng2024direct}
Deng, Y., Wang, Y. \harvardand\ Zhou, X.-H. \harvardyearleft
2024\harvardyearright , `Direct and indirect treatment effects in the
presence of semicompeting risks', {\em Biometrics} {\bf 80}(2),~ujae032.
\harvarditem{D{\'\i}az \harvardand\ Hejazi}{2020}{diaz2020causal}
D{\'\i}az, I. \harvardand\ Hejazi, N.~S. \harvardyearleft
2020\harvardyearright , `Causal mediation analysis for stochastic
interventions', {\em Journal of the Royal Statistical Society Series B:
Statistical Methodology} {\bf 82}(3),~661--683.
\harvarditem[D'iaz et~al.]{D'iaz, Hoffman \harvardand\
Hejazi}{2022}{Diaz2022Causal}
D'iaz, I., Hoffman, K. \harvardand\ Hejazi, N. \harvardyearleft
2022\harvardyearright , `Causal survival analysis under competing risks using
longitudinal modified treatment policies', {\em Lifetime Data Analysis} {\bf
30},~213--236.
\harvarditem{{Emerging Risk Factors Collaboration}
et~al.}{2010}{emerging2010diabetes}
{Emerging Risk Factors Collaboration} et~al. \harvardyearleft
2010\harvardyearright , `Diabetes mellitus, fasting blood glucose
concentration, and risk of vascular disease: a collaborative meta-analysis of
102 prospective studies', {\em The lancet} {\bf 375}(9733),~2215--2222.
\harvarditem{Fine \harvardand\ Gray}{1999}{fine1999proportional}
Fine, J.~P. \harvardand\ Gray, R.~J. \harvardyearleft 1999\harvardyearright ,
`A proportional hazards model for the subdistribution of a competing risk',
{\em Journal of the American Statistical Association} {\bf
94}(446),~496--509.
\harvarditem[Fine et~al.]{Fine, Jiang \harvardand\
Chappell}{2001}{fine2001semi}
Fine, J.~P., Jiang, H. \harvardand\ Chappell, R. \harvardyearleft
2001\harvardyearright , `On semi-competing risks data', {\em Biometrika} {\bf
88}(4),~907--919.
\harvarditem{Forbes \harvardand\ Cooper}{2013}{forbes2013mechanisms}
Forbes, J.~M. \harvardand\ Cooper, M.~E. \harvardyearleft
2013\harvardyearright , `Mechanisms of diabetic complications', {\em
Physiological Reviews} {\bf 93}(1),~137--188.
\harvarditem{Grambsch \harvardand\ Therneau}{1994}{grambsch1994proportional}
Grambsch, P.~M. \harvardand\ Therneau, T.~M. \harvardyearleft
1994\harvardyearright , `Proportional hazards tests and diagnostics based on
weighted residuals', {\em Biometrika} {\bf 81}(3),~515--526.
\harvarditem[Gu et~al.]{Gu, Zeng, Heiss \harvardand\ Lin}{2024}{gu2024maximum}
Gu, Y., Zeng, D., Heiss, G. \harvardand\ Lin, D. \harvardyearleft
2024\harvardyearright , `Maximum likelihood estimation for semiparametric
regression models with interval-censored multistate data', {\em Biometrika}
{\bf 111}(3),~971--988.
\harvarditem[Ha et~al.]{Ha, Xiang, Peng, Jeong \harvardand\
Lee}{2020}{ha2020frailty}
Ha, I.~D., Xiang, L., Peng, M., Jeong, J.-H. \harvardand\ Lee, Y.
\harvardyearleft 2020\harvardyearright , `Frailty modelling approaches for
semi-competing risks data', {\em Lifetime Data Analysis} {\bf 26},~109--133.
\harvarditem[Hejazi et~al.]{Hejazi, Rudolph, Van Der~Laan \harvardand\
Díaz}{2022}{hejazi2022nonparametric}
Hejazi, N.~S., Rudolph, K.~E., Van Der~Laan, M.~J. \harvardand\ Díaz, I.
\harvardyearleft 2022\harvardyearright , `Nonparametric causal mediation
analysis for stochastic interventional (in)direct effects', {\em
Biostatistics} {\bf 24}(3),~686--707.
\harvarditem[Hines et~al.]{Hines, Dukes, Diaz-Ordaz \harvardand\
Vansteelandt}{2022}{hines2022demystifying}
Hines, O., Dukes, O., Diaz-Ordaz, K. \harvardand\ Vansteelandt, S.
\harvardyearleft 2022\harvardyearright , `Demystifying statistical learning
based on efficient influence functions', {\em The American Statistician} {\bf
76}(3),~292--304.
\harvarditem[Hirano et~al.]{Hirano, Imbens \harvardand\
Ridder}{2003}{hirano2003efficient}
Hirano, K., Imbens, G.~W. \harvardand\ Ridder, G. \harvardyearleft
2003\harvardyearright , `Efficient estimation of average treatment effects
using the estimated propensity score', {\em Econometrica} {\bf
71}(4),~1161--1189.
\harvarditem{Hougaard}{1999}{hougaard1999multi}
Hougaard, P. \harvardyearleft 1999\harvardyearright , `Multi-state models: a
review', {\em Lifetime Data Analysis} {\bf 5},~239--264.
\harvarditem{Huang}{2021}{huang2021causal}
Huang, Y.-T. \harvardyearleft 2021\harvardyearright , `Causal mediation of
semicompeting risks', {\em Biometrics} {\bf 77}(4),~1143--1154.
\harvarditem{Kennedy}{2024}{kennedy2024semiparametric}
Kennedy, E.~H. \harvardyearleft 2024\harvardyearright , `Semiparametric doubly
robust targeted double machine learning: a review', {\em Handbook of
Statistical Methods for Precision Medicine} pp.~207--236.
\harvarditem{Kodell \harvardand\ Nelson}{1980}{kodell1980illness}
Kodell, R. \harvardand\ Nelson, C. \harvardyearleft 1980\harvardyearright ,
`An illness--death model for the study of the carcinogenic process using
survival/sacrifice data.', {\em Biometrics} {\bf 36}(2),~267--277.
\harvarditem[Kormaksson et~al.]{Kormaksson, Lange, Demanse, Strohmaier, Duan,
Xie, Carbini, Bossen, Guettner \harvardand\
Maniero}{2024}{Kormaksson2024Dynamic}
Kormaksson, M., Lange, M.~R., Demanse, D., Strohmaier, S., Duan, J., Xie, Q.,
Carbini, M., Bossen, C., Guettner, A. \harvardand\ Maniero, A.
\harvardyearleft 2024\harvardyearright , `Dynamic path analysis for exploring
treatment effect mediation processes in clinical trials with
time‐to‐event endpoints', {\em Statistics in Medicine} {\bf 43},~4614 --
4634.
\harvarditem[Li et~al.]{Li, Zhang, Bakoyannis \harvardand\
Gao}{2020}{li2020shared}
Li, J., Zhang, Y., Bakoyannis, G. \harvardand\ Gao, S. \harvardyearleft
2020\harvardyearright , `On shared gamma-frailty conditional {Markov} model
for semicompeting risks data', {\em Statistics in Medicine} {\bf
39}(23),~3042--3058.
\harvarditem{Lin \harvardand\ VanderWeele}{2017}{lin2017interventional}
Lin, S.-H. \harvardand\ VanderWeele, T. \harvardyearleft 2017\harvardyearright
, `Interventional approach for path-specific effects', {\em Journal of Causal
Inference} {\bf 5}(1),~20150027.
\harvarditem[Lin et~al.]{Lin, Young, Logan \harvardand\
VanderWeele}{2017}{lin2017mediation}
Lin, S.-H., Young, J.~G., Logan, R. \harvardand\ VanderWeele, T.~J.
\harvardyearleft 2017\harvardyearright , `Mediation analysis for a survival
outcome with time-varying exposures, mediators, and confounders', {\em
Statistics in Medicine} {\bf 36}(26),~4153--4166.
\harvarditem{Lok}{2015}{lok2015defining}
Lok, J.~J. \harvardyearleft 2015\harvardyearright , `Defining and estimating
causal direct and indirect effects when setting the mediator to specific
values is not feasible', {\em Statistics in Medicine} {\bf
35}(22),~4008--4020.
\harvarditem[Marso et~al.]{Marso, Daniels, Brown-Frandsen, Kristensen, Mann,
Nauck, Nissen, Pocock, Poulter, Ravn et~al.}{2016}{marso2016liraglutide}
Marso, S.~P., Daniels, G.~H., Brown-Frandsen, K., Kristensen, P., Mann, J.~F.,
Nauck, M.~A., Nissen, S.~E., Pocock, S., Poulter, N.~R., Ravn, L.~S. et~al.
\harvardyearleft 2016\harvardyearright , `Liraglutide and cardiovascular
outcomes in type 2 diabetes', {\em New England Journal of Medicine} {\bf
375}(4),~311--322.
\harvarditem[Marso et~al.]{Marso, Poulter, Nissen, Nauck, Zinman, Daniels,
Pocock, Steinberg, Bergenstal, Mann et~al.}{2013}{marso2013design}
Marso, S.~P., Poulter, N.~R., Nissen, S.~E., Nauck, M.~A., Zinman, B., Daniels,
G.~H., Pocock, S., Steinberg, W.~M., Bergenstal, R.~M., Mann, J.~F. et~al.
\harvardyearleft 2013\harvardyearright , `Design of the liraglutide effect
and action in diabetes: evaluation of cardiovascular outcome results
{(LEADER)} trial', {\em American Heart Journal} {\bf 166}(5),~823--830.
\harvarditem{Martinussen \harvardand\
Stensrud}{2023}{martinussen2023estimation}
Martinussen, T. \harvardand\ Stensrud, M.~J. \harvardyearleft
2023\harvardyearright , `Estimation of separable direct and indirect effects
in continuous time', {\em Biometrics} {\bf 79}(1),~127--139.
\harvarditem{McLachlan \harvardand\ Krishnan}{2008}{mclachlan2008algorithm}
McLachlan, G.~J. \harvardand\ Krishnan, T. \harvardyearleft
2008\harvardyearright , {\em The EM algorithm and extensions}, John Wiley \&
Sons.
\harvarditem{Miles}{2023}{miles2022on}
Miles, C.~H. \harvardyearleft 2023\harvardyearright , `On the causal
interpretation of randomised interventional indirect effects', {\em Journal
of the Royal Statistical Society Series B: Statistical Methodology} {\bf
85}(4),~1154--1172.
\harvarditem[Miles et~al.]{Miles, Shpitser, Kanki, Meloni \harvardand\
{Tchetgen Tchetgen}}{2020}{miles2020semiparametric}
Miles, C.~H., Shpitser, I., Kanki, P., Meloni, S. \harvardand\ {Tchetgen
Tchetgen}, E.~J. \harvardyearleft 2020\harvardyearright , `On semiparametric
estimation of a path-specific effect in the presence of mediator-outcome
confounding', {\em Biometrika} {\bf 107}(1),~159--172.
\harvarditem[Prentice et~al.]{Prentice, Kalbfleisch, Peterson~Jr, Flournoy,
Farewell \harvardand\ Breslow}{1978}{prentice1978analysis}
Prentice, R.~L., Kalbfleisch, J.~D., Peterson~Jr, A.~V., Flournoy, N.,
Farewell, V.~T. \harvardand\ Breslow, N.~E. \harvardyearleft
1978\harvardyearright , `The analysis of failure times in the presence of
competing risks', {\em Biometrics} pp.~541--554.
\harvarditem[Putter et~al.]{Putter, Fiocco \harvardand\
Geskus}{2007}{putter2007tutorial}
Putter, H., Fiocco, M. \harvardand\ Geskus, R.~B. \harvardyearleft
2007\harvardyearright , `Tutorial in biostatistics: competing risks and
multi-state models', {\em Statistics in Medicine} {\bf 26}(11),~2389--2430.
\harvarditem{Ripatti \harvardand\ Palmgren}{2004}{ripatti2004estimation}
Ripatti, S. \harvardand\ Palmgren, J. \harvardyearleft 2004\harvardyearright ,
`Estimation of multivariate frailty models using penalized partial
likelihood', {\em Biometrics} {\bf 56}(4),~1016--1022.
\harvarditem{Robins}{2004}{robins2004optimal}
Robins, J.~M. \harvardyearleft 2004\harvardyearright , Optimal structural
nested models for optimal sequential decisions, {\em in} `Proceedings of the
Second Seattle Symposium in Biostatistics: Analysis of correlated data',
Springer, pp.~189--326.
\harvarditem[Robins et~al.]{Robins, Richardson \harvardand\
Shpitser}{2022}{robins2022interventionist}
Robins, J.~M., Richardson, T.~S. \harvardand\ Shpitser, I. \harvardyearleft
2022\harvardyearright , An interventionist approach to mediation analysis,
{\em in} `Probabilistic and causal inference: the works of Judea Pearl', ACM,
pp.~713--764.
\harvarditem[Rotolo et~al.]{Rotolo, Rondeau \harvardand\
Legrand}{2016}{rotolo2016incorporation}
Rotolo, F., Rondeau, V. \harvardand\ Legrand, C. \harvardyearleft
2016\harvardyearright , `Incorporation of nested frailties into
semiparametric multi-state models', {\em Statistics in Medicine} {\bf
35}(4),~609--621.
\harvarditem{Shpitser \harvardand\ Tchetgen~Tchetgen}{2016}{shpitser2016causal}
Shpitser, I. \harvardand\ Tchetgen~Tchetgen, E. \harvardyearleft
2016\harvardyearright , `Causal inference with a graphical hierarchy of
interventions', {\em Annals of Statistics} {\bf 44}(6),~2433.
\harvarditem[Stensrud et~al.]{Stensrud, Hern{\'a}n, Tchetgen~Tchetgen, Robins,
Didelez \harvardand\ Young}{2021}{stensrud2021generalized}
Stensrud, M.~J., Hern{\'a}n, M.~A., Tchetgen~Tchetgen, E.~J., Robins, J.~M.,
Didelez, V. \harvardand\ Young, J.~G. \harvardyearleft 2021\harvardyearright
, `A generalized theory of separable effects in competing event settings',
{\em Lifetime Data Analysis} {\bf 27}(4),~588--631.
\harvarditem[Stensrud et~al.]{Stensrud, Young, Didelez, Robins \harvardand\
Hern{\'a}n}{2022}{stensrud2022separable}
Stensrud, M.~J., Young, J.~G., Didelez, V., Robins, J.~M. \harvardand\
Hern{\'a}n, M.~A. \harvardyearleft 2022\harvardyearright , `Separable
effects for causal inference in the presence of competing events', {\em
Journal of the American Statistical Association} {\bf 117}(537),~175--183.
\harvarditem[Stratton et~al.]{Stratton, Adler, Neil, Matthews, Manley, Cull,
Hadden, Turner \harvardand\ Holman}{2000}{stratton2000association}
Stratton, I.~M., Adler, A.~I., Neil, H. A.~W., Matthews, D.~R., Manley, S.~E.,
Cull, C.~A., Hadden, D., Turner, R.~C. \harvardand\ Holman, R.~R.
\harvardyearleft 2000\harvardyearright , `Association of glycaemia with
macrovascular and microvascular complications of type 2 diabetes ({UKPDS
35}): prospective observational study', {\em BMJ} {\bf 321}(7258),~405--412.
\harvarditem[Tai et~al.]{Tai, Lin, Chu, Yu, Puhan \harvardand\
VanderWeele}{2023}{tai2023causal}
Tai, A.-S., Lin, S.-H., Chu, Y.-C., Yu, T., Puhan, M.~A. \harvardand\
VanderWeele, T. \harvardyearleft 2023\harvardyearright , `Causal mediation
analysis with multiple time-varying mediators', {\em Epidemiology} {\bf
34}(1),~8--19.
\harvarditem[Valeri et~al.]{Valeri, Proust-Lima, Fan, Chen \harvardand\
Jacqmin-Gadda}{2023}{valeri2023multistate}
Valeri, L., Proust-Lima, C., Fan, W., Chen, J.~T. \harvardand\ Jacqmin-Gadda,
H. \harvardyearleft 2023\harvardyearright , `A multistate approach for the
study of interventions on an intermediate time-to-event in health disparities
research', {\em Statistical Methods in Medical Research} {\bf
32}(8),~1445--1460.
\harvarditem{Van~der Vaart}{2000}{van2000asymptotic}
Van~der Vaart, A.~W. \harvardyearleft 2000\harvardyearright , {\em Asymptotic
statistics}, Vol.~3, Cambridge university press.
\harvarditem{van~der vaart \harvardand\ Wellner}{2013}{van2013weak}
van~der vaart, A. \harvardand\ Wellner, J. \harvardyearleft
2013\harvardyearright , {\em Weak Convergence and Empirical Processes: With
Applications to Statistics}, Springer Series in Statistics, Springer New
York.
\harvarditem{VanderWeele \harvardand\
Tchetgen~Tchetgen}{2017}{vanderweele2017mediation}
VanderWeele, T.~J. \harvardand\ Tchetgen~Tchetgen, E.~J. \harvardyearleft
2017\harvardyearright , `Mediation analysis with time varying exposures and
mediators', {\em Journal of the Royal Statistical Society Series B:
Statistical Methodology} {\bf 79}(3),~917--938.
\harvarditem{Vansteelandt \harvardand\
Daniel}{2017}{vansteelandt2017interventional}
Vansteelandt, S. \harvardand\ Daniel, R.~M. \harvardyearleft
2017\harvardyearright , `Interventional effects for mediation analysis with
multiple mediators', {\em Epidemiology} {\bf 28}(2),~258--265.
\harvarditem[Vansteelandt et~al.]{Vansteelandt, Linder, Vandenberghe, Steen
\harvardand\ Madsen}{2019}{vansteelandt2019Mediation}
Vansteelandt, S., Linder, M., Vandenberghe, S., Steen, J. \harvardand\ Madsen,
J. \harvardyearleft 2019\harvardyearright , `Mediation analysis of
time‐to‐event endpoints accounting for repeatedly measured mediators
subject to time‐varying confounding', {\em Statistics in Medicine} {\bf
38}(24),~4828 -- 4840.
\harvarditem[Vo et~al.]{Vo, Williams, Liu, Rudolph \harvardand\
D{\'\i}az}{2026}{vo2026recanting}
Vo, T.-T., Williams, N., Liu, R., Rudolph, K.~E. \harvardand\ D{\'\i}az, I.
\harvardyearleft 2026\harvardyearright , `Recanting twins: Addressing
intermediate confounding in mediation analysis', {\em Statistics in Medicine}
{\bf 45}(3-5),~e70432.
\harvarditem[Wang et~al.]{Wang, Laan, Petersen, Gerds, Kvist \harvardand\
Laan}{2025}{wang2025targeted}
Wang, Z., Laan, L. v.~d., Petersen, M., Gerds, T., Kvist, K. \harvardand\ Laan,
M. v.~d. \harvardyearleft 2025\harvardyearright , `Targeted maximum
likelihood based estimation for longitudinal mediation analysis', {\em
Journal of Causal Inference} {\bf 13}(1),~20230013.
\harvarditem[Weir et~al.]{Weir, Rider \harvardand\
Trinquart}{2022}{weir2022counterfactual}
Weir, I.~R., Rider, J.~R. \harvardand\ Trinquart, L. \harvardyearleft
2022\harvardyearright , `Counterfactual mediation analysis in the multistate
model framework for surrogate and clinical time-to-event outcomes in
randomized controlled trials', {\em Pharmaceutical Statistics} {\bf
21}(1),~163--175.
\harvarditem[Xu et~al.]{Xu, Kalbfleisch \harvardand\
Tai}{2010}{xu2010statistical}
Xu, J., Kalbfleisch, J.~D. \harvardand\ Tai, B. \harvardyearleft
2010\harvardyearright , `Statistical analysis of illness--death processes and
semicompeting risks data', {\em Biometrics} {\bf 66}(3),~716--725.
\harvarditem{Zeng \harvardand\ Lin}{2007}{zeng2007maximum}
Zeng, D. \harvardand\ Lin, D. \harvardyearleft 2007\harvardyearright ,
`Maximum likelihood estimation in semiparametric regression models with
censored data', {\em Journal of the Royal Statistical Society Series B:
Statistical Methodology} {\bf 69}(4),~507--564.
\end{thebibliography}
\newpage