
\section{Introduction}

In survival analysis, a competing event is one whose realization precludes the
occurrence of the event of interest. As an example, consider a study on the
effect of mechanical ventilation of hospitalized COVID-19 patients on the onset
of acute kidney injury, where the death of a patient by extraneous causes may be
considered a competing event. Commonly used methods for analyzing time-to-event
outcomes under competing events include estimating the cause-specific hazard
function~\citep{prentice1978analysis}, the sub-distribution function~\citep[also
known as cause-specific cumulative incidence;][]{fine1999proportional}, or the
marginal cumulative incidence that treats the competing event as a censoring
event. Early techniques for the estimation of these parameters were developed in
the context of the Cox proportional hazards models~\citep{cox1972regression,
andersen1982cox}.

While these methods have allowed much progress in applied research, they have
two important limitations: (i) they lack a formal framework that clarifies the
conditions required for a causal interpretation of the estimates, and (ii) they
typically rely upon parametric modeling assumptions that can induce
non-negligible bias in the effect estimates, even when the assumptions required
for causal identification hold. Recent developments in causal inference and
semiparametric estimation have yielded important progress in addressing these
significant limitations. For example, when the exposure is
binary,~\citet{young2020causal} present a study on the definition and
identification of several causal effect parameters in the presence of competing
risks, and~\citet{benkeser2018improved} and~\citet{rytgaard2021one} describe
semiparametric--efficient estimators that leverage the flexibility of
data-adaptive regression procedures in order to mitigate model misspecification
bias while retaining asymptotic properties critical to reliable statistical
inference (e.g., $\sqrt{n}$-consistency). Among other
insights,~\cite{young2020causal} clarify that identification of the marginal
cumulative incidence entails considering interventions that eliminate the
competing event, and therefore require no unmeasured confounding of the
competing event--outcome relation. This assumption is not required to identify
the cause-specific cumulative incidence, which is the focus of this work.

Defining causal effects requires considering counterfactual outcomes under
hypothetical interventions on the exposure of interest. In the case of binary
exposures, causal effects are often defined as differences between the expected
outcomes under hypothetical interventions that assign the exposure to all units
or to none --- that is, the average treatment effect. While useful, causal
effects defined through such deterministic interventions are limited in many
applications, including problems where the exposure of interest is numerical
rather than categorical, where the exposure is multivariate, where the exposure
is a time-to-event-variable, or when exposure assignment is deterministic
conditional on potential confounders (cf.~positivity/overlap assumption).

As a solution to these problems, recent work has focused on stochastic
interventions and modified treatment
policies~\citep{stock1989nonparametric, robins2004effects, Diaz12,
  Haneuse2013, young2014identification}.  Stochastic interventions
inquire what would have happened in a hypothetical world where the
post-intervention exposure is a random draw from a user-specified
distribution that may itself depend on the observed data
distribution. Examples of stochastic interventions include
interventional mediation effects~\citep{vanderweele2014effect,
  diaz2021nonparametric, hejazi2022nonparametric}, as well as
incremental propensity score interventions that address violations of
the positivity assumption~\citep{kennedy2019nonparametric,
  wen2021intervention}. Modified treatment policies (MTPs), on the
other hand, are defined as interventions that depend on the natural
value of treatment~\citep{young2014identification}. MTPs allow the
definition and estimation of parameters defined in terms of natural
research questions, such as inquiring what the incidence of wheezing
among children with asthma would have been had the concentration of
PM$_{2.5}$ particles been 5\% lower, or what the mortality rate due to
opioid overdose would have been had naloxone access laws been
implemented earlier~\citep{rudolph2021effects}.

Identifying assumptions for causal effects of longitudinal modified
treatment policies were first articulated
by~\citet{richardson2013single}
and~\citet{young2014identification}. Recently,~\citet{diaz2021lmtp}
developed sequential regression and sequentially doubly robust
estimators for longitudinal modified treatment policies (LMTP). These
estimators are generalizations of estimators developed
by~\citet{diaz2018stochastic} for MTPs in the single time-point case,
and of estimators developed by~\citet{luedtke2017sequential}
and~\citet{rotnitzky2017multiply} for static interventions in the
multiple time-point case. Here, we extend the LMTP methodology for the
estimation of the cause-specific cumulative incidence under competing
events, illustrating how LMTPs can be used to solve problems whose
solution has remained elusive in causal inference, such as the
definition and estimation of effects for continuous, multi-valued,
time-valued, or multivariate exposures. In our motivating example we
illustrate estimation of the effect of time-to-intubation on time to
acute kidney injury, with death by other causes acting as a competing
risk.

Our estimation strategy is rooted in recent developments allowing the
estimation of causal effects using flexible regression techniques from
the machine learning literature~\citep{vanderLaanRose11,
  vanderLaanRose18, chernozhukov2017double}. Central to this theory is
the notion of von-Mises expansion~\citep[e.g.,][]{mises1947asymptotic,
  vanderVaart98, robins2009quadratic}, which together with
sample-splitting techniques~\citep{bickel1982adaptive,
  klaassen1987consistent, zheng2011cross} allows the use of flexible
regression methods for the estimation of nuisance parameters while
still allowing the construction of efficient and $\sqrt{n}$-consistent
estimators.

\section{Defining LMTP effects under competing risks}\label{sec:notation}

Assume that at each time point $t=1,\ldots\tau$, we measure
$(D_t, Y_t, L_t, A_t, C_t)$ on each of $i=1,\ldots, n$ study
participants, where $L_t$ denotes time-varying covariates, $A_t$
denotes a vector of exposure variables at time $t$, $C_t$ denotes an
indicator of loss-to-follow-up (equal to one if the unit remains
on study at time $t+1$ and equal to zero otherwise), $D_t$ is an
indicator of the competing event occurring on or before time $t$, and
$Y_t$ denotes an indicator of the event of interest occurring on or
before time $t$. At the end of the study, we measure the outcome of
interest $Y_{\tau+1}$. We let $D_1=Y_1=0$, meaning that all
participants have not experienced the event of interest or the
competing event at the beginning of the study. We let
$R_t=\mathds{1}\{D_t=0, Y_t=0\}$ denote an indicator that neither the event
of interest nor the competing event have occurred by $t$.  Assume
monotone loss-to-follow-up so that $C_t=0$ implies $C_k=0$ for all
$k>t$, in which case all data become degenerate for $k > t$. In this
work we are interested in a scenario where occurrence of the competing
event precludes occurrence of the event of interest. For example, in a
study looking to assess the effect of hypertension medications on the
likelihood of a stroke, death for reasons other than a stroke is a
competing event. In this notation, we have $Y_t=1$ if $Y_{t-1}=1$,
$D_t=1$ if $D_{t-1}=1$, and $Y_t=0$ if $D_t=1$ and $Y_{t-1}=0$, since
the competing event precludes occurrence of the event of interest. Let
$X_i$ denote the observed data for subject $i$. We let
$\P f = \int f(x)\dd \P(x)$ for a given function $f(x)$. We use $\mathsf{P}_n$
to denote the empirical distribution of $X_1,\ldots\,X_n$, and assume
$\P$ is an element of the nonparametric statistical model defined as
all continuous densities on $X$ with respect to a dominating measure
$\nu$. We let $\mathsf{E}$ denote the expectation with respect to $\P$, i.e.,
$\mathsf{E}\{f(X)\} = \int f(x)\dd\P(x)$. We also let $||f||^2$ denote the
$L_2(\P)$ norm $\int f^2(x)\dd\P(x)$. We use
$\bar W_t = (W_1,\ldots, W_t)$ to denote the past history of a
variable $W$, use $\ubar W_t = (W_t,\ldots, W_\tau)$ to denote the
future of a variable, and use $H_t = (\bar L_t, \bar A_{t-1})$ to
denote the history of the time-varying covariates and the treatment
process up until just before time $t$.  $\bar W$ denotes the entire
history $(W_1,\ldots,W_\tau)$ of a variable. By convention, variables
with an index $t\leq 0$ are defined as the null set, expectations
conditioning on a null set are marginal, products of the type
$\prod_{t=k}^{k-1}b_t$ and $\prod_{t=0}^0b_t$ are equal to one, and
sums of the type $\sum_{t=k}^{k-1}b_t$ and $\sum_{t=0}^0b_t$ are equal
to zero. For vectors $u$ and $v$, $u\le v$ denotes point-wise
inequalities.

We formalize the definition of the causal effects using a
non-parametric structural equation model
\citep{Pearl00}. Specifically, for each time point $t$, we assume the
existence of deterministic functions $f_{D,t}$, $f_{Y,t}$, $f_{L,t}$,
$f_{A,t}$, $f_{C,t}$, such that
\begin{align*}
  D_t&=f_{D,t}(C_{t-1}, A_{t-1}, H_{t-1}, U_{D,t})\\
  Y_t&=f_{Y,t}(D_t, C_{t-1}, A_{t-1}, H_{t-1}, U_{Y,t})\\
  L_t&=f_{L,t}(Y_t, D_t, C_{t-1}, A_{t-1}, H_{t-1}, U_{Y,t})\\
  A_t&=f_{A,t}(H_t, U_{A,t})\\
  C_t&=f_{C,t}(A_t, H_t, U_{C,t})
\end{align*}
Here $U=(U_{L,t}, U_{A,t}, U_{C,t}, U_{D,t},U_{Y,t}:t\in \{1,\ldots,\tau+1\})$
is a vector of exogenous variables with unrestricted joint distribution. Without
loss of generality, we assume that variables which are undefined are equal to
zero (e.g., $L_t=0$ if $C_k=1$ for any $k<t$).



We are concerned with the definition and estimation of the causal effect of an
intervention on the treatment process $\bar A$ on the event indicator
$Y_{\tau+1}$ in a hypothetical world where there is no loss to follow-up (i.e.,
$\bar C=0$). The interventions on the treatment process are defined in terms of
longitudinal modified treatment policies \citep{Diaz12,Haneuse2013}, which are
hypothetical interventions where the exposure is assigned as a new random
variable $A_t^\d$ instead of being assigned according to the structural
equation model. An intervention that sets the exposures up to time $t-1$ to $\bar
A_{t-1}^\d$ generates a counterfactual variable $A_t(\bar A_{t-1}^\d)$, which is
referred to as the \textit{natural value of treatment}
\citep{richardson2013single, young2014identification}, and represents the value
of treatment that would have been observed at time $t$ under an intervention
carried out up until time $t-1$ but discontinued thereafter. An intervention on
all the treatment and censoring variables up to $t=\tau$ generates
a counterfactual outcome $Y_{\tau+1}(\bar A^\d)$. Causal effects are defined as
contrasts between the counterfactual probability $\P[Y_{\tau+1}(\bar
A^\d)=1]$ of a failure event by the end
of the study  under interventions implied by
different interventions $\d$, as defined below.

The causal effects we define are characterized by a user-given function
$\d(a_t, h_t,\varepsilon_t)$ that maps a given treatment value $a_t$ (i.e., the
``natural value'') and a history $h_t$ into a new exposure value. The function
$\d$ is also allowed to depend on a randomizer $\varepsilon_t$ drawn
independently across units and independently of $U$, and whose distribution does
not depend on $\P$. For fixed values $\bar a_t$, and $\bar l_t$, we recursively
define $a_t^\d=\d(a_t, h^\d_t, \varepsilon_t)$, where $h^\d_t=(\bar a_{t-1}^\d,
\bar l_t)$. The intervention is undefined if $h^\d_t$ is such that $d_s^\d=1$ or
$y_s^\d=1$ for any $s\leq t$. The LMTPs that we study are thus defined as
$A^\d_t = \d(A_t(\bar A_{t-1}^\d), H_t(\bar A_{t-1}^\d), \varepsilon_t)$ for
a user-given function $\d$. In what follows, we let $\mathsf{g}_{A,t}(a_t \mid h_t)$
denote the density or probability mass function of $A_t$ conditional on
$C_{t-1}=R_t=1$ and $H_t=h_t$, and we let $\mathsf{g}_{C,t}(a_t, h_t)$ denote
$\P(C_t=1\mid C_{t-1}=R_t=1, A_t=a_t, H_t=h_t)$. Furthermore, we use
$\mathsf{g}_t^\d(a_t \mid h_t)$ to denote the density or probability mass function of
$\d(A_t,H_t,\varepsilon_t)$ conditional on $C_{t-1}=R_t=1$ and $H_t=h_t$.

To ground the ideas and illustrate the usefulness of this framework, we now
consider several examples of regimes $\d$ that yield interesting and
scientifically informative causal contrasts. More examples appear in
\citet{robins2004effects,young2014identification,diaz2021lmtp}.

\begin{example}[Additive shift LMTP]
  Let $A_t$ denote a vector of numerical exposures, such as a drug dose or
  energy expenditure through physical activity. Assume that $A_t$ is supported
  as $\P(A_t \leq u_t(h_t)\mid H_t = h_t)=1$ for some $u_t$. Then, for
  a user-given vector $\delta$, we let
  \begin{equation}\label{eq:defdshift}
    \d(a_t,h_t)=
    \begin{cases}
      a_t + \delta & \text{if } a_t+\delta \leq u_t(h_t) \\
      a_t & \text{if } a_t +\delta > u_t(h_t).
    \end{cases}
  \end{equation} We define
  $A_t^\d=\d(A_t(\bar A_{t-1}^\d), H_t(\bar A_{t-1}^\d))$. This
  intervention has been discussed by~\citet{Diaz12}, \citet{diaz2018stochastic},
  and~\citet{Haneuse2013} in the context of single time-point exposures;
  by~\citet{diaz2020causal} and~\citet{hejazi2022nonparametric} in the
  context of causal mediation analysis; and by~\citet{hejazi2020efficient} in
  the context of studies with two-phase sampling designs. This intervention
  considers hypothetical worlds in which the natural exposure at time $t$ is
  increased by a user-given value $\delta$, whenever such an increase is
  feasible for a unit with history $H_t(A_{t-1}^\d)$.
\end{example}

\begin{example}[Multiplicative shift LMTP]\label{ex:mshift}
  Let $A_t$ denote a vector of numerical exposures such as pollutant
  concentrations for various pollutants. To define this intervention, assume
  that $A_t$ is supported as $\P(A_t \geq l_t(h_t)\mid H_t = h_t)=1$ for some
  $l_t$. Then, for a user-given value $0 < \delta < 1$, we let
  \begin{equation}\label{eq:mshift}
    \d(a_t,h_t)=
    \begin{cases}
      a_t \times \delta & \text{if } a_t \times \delta \geq l_t(h_t)  \\
      a_t & \text{if } a_t \times \delta < l_t(h_t).
    \end{cases}
  \end{equation}
  The intervention implied by this function considers hypothetical
  worlds where the natural exposure at time $t$ is increased by a
  user-given factor $\delta$, whenever such increase is feasible for a
  unit with history $H_t(A_{t-1}^\d)$. This intervention seems useful,
  for example, to solve current problems in environmental epidemiology
  related to the effect of environmental mixtures
  \citep{gibson2019overview}.
\end{example}

\begin{example}[Incremental propensity score interventions based on
  the risk ratio]\label{ex:ipsi}
  \citet{kennedy2019nonparametric} proposed an intervention tailored to binary
  exposures in which the post-intervention exposure is a draw from a user-given
  distribution $\mathsf{g}_\delta(a_t\mid h_t)$, where this distribution is chosen such
  that the odds ratio comparing the likelihood of exposure under intervention
  versus the actual exposure is a user-given value $\delta$. As identifying the
  effect of this intervention does not require the positivity assumption, these
  incremental propensity score interventions are very useful in settings where
  violations of the positivity assumption preclude identification of static
  effects such as the average treatment effect.

  Recently, \cite{wen2021intervention} proposed a similar intervention that
  instead uses the risk ratio to quantify the likelihood of exposure under the
  intervention. This intervention may be expressed as an LMTP as follows
  \begin{equation}\label{eq:defipsi}
    \d(a_t,\varepsilon_t)=
    \begin{cases}
      a_t & \text{if } \varepsilon_t < \delta  \\
      0 & \text{otherwise},
    \end{cases}
  \end{equation}
  where $\varepsilon_t$ is a random draw from the uniform distribution on
  $(0,1)$. As before, we make use of the definition
  $A_t^\d=\d(A_t(\bar A_{t-1}^\d), H_t(\bar A_{t-1}^\d))$. Note that the
  probability of exposure under the intervention is equal to $\mathsf{g}_\delta(1\mid
  h_t) = \delta \times \mathsf{g}(1\mid h_t)$; thus, the risk ratio of exposure under
  the intervention versus the natural exposure value is $\delta$. Importantly,
  a unit with zero probability of exposure conditional on their history also has
  zero probability of exposure under the intervention. This weakens the
  positivity assumptions required for identification and estimation, as we will
  discuss below.
\end{example}

\begin{example}[Dynamic treatment initiation strategies with grace
  period]\label{ex:grace}
  Let $A_t$ denote a binary treatment variable.  Consider an intervention of the
  form ``for a user-given grace period $m$, if a condition for treatment
  initiation is met at or before $t-m$ then start treatment at $t$; otherwise,
  do not intervene.'' Letting $L'_t$ denote an indicator that the condition has
  been met, this intervention can be operationalized as an LMTP as follows:
  \begin{equation}\label{eq:defgrace}
    \d(a_t,h_t)=
    \begin{cases}
      1 & \text{if } l'_{t-m} = 1  \\
      a_t & \text{otherwise},
    \end{cases}
  \end{equation}
  This type of intervention has been considered in the context of treatment
  initiation for HIV~\citep{van2005history,cain2010start}, where the condition
  $L_t'$ is that CD4+ cell count drops below a pre-specified value. In this type
  of application, realistic treatment policies may involve a delay of treatment
  initiation for logistical reasons. Accordingly, allowing for a grace period
  quantifies an effect that is arguably of more practical relevance.
\end{example}

\begin{remark}
  In contrast to LMTPs, stochastic interventions proceed by specifying
  an exposure density $\mathsf{g}_t^\d(a_t\mid h_t)$ \textit{a priori} and considering
  a sequence of draws $Q^\d_t$ from $\mathsf{g}_t^\d(a_t\mid h_t)$. Counterfactual
  outcomes under stochastic interventions can then be defined as
  $Y(\bar Q^\d)$. \citet{wen2021intervention} study stochastic
  interventions where the post-intervention density can be
  expressed as a mixture between a known density $\mathbbm{f}_t(a_t\mid h_t)$
  and the density of the observed data $\mathsf{g}_t(a_t\mid h_t)$. Any such
  intervention has an LMTP counterpart defined by
  \begin{equation*}
    \d(a_t,h_t,\varepsilon_t)=\mathds{1}\{\varepsilon_t\leq c_t(h_t)\} q_t +
    \mathds{1}\{\varepsilon_t> c_t(h_t)\} a_t,
  \end{equation*}
  where $\varepsilon_t$ is a random draw from a uniform distribution in
  $[0,1]$, $0<c_t(h_t)<1$ is a mixture constant, and $q_t$ represents a
  draw from $\mathbbm{f}_t(a_t\mid h_t)$. Efficient estimators for this type of
  LMTP have been previously proposed
  by~\cite{diaz2021lmtp}. \citet{wen2021intervention} discuss the
  stochastic interventions implied by Examples~\ref{ex:ipsi}
  and~\ref{ex:grace}, but incorrectly claim that the estimators they
  develop are not particular cases of the methods of~\cite{diaz2021lmtp}.
\end{remark}

\subsection{Identification of the effect of LMTPs under competing risks}

The following assumptions will be sufficient to prove identification:

\begin{assumptioniden}[Supported exposure intervention]\label{ass:support1}
  If $(a_t,h_t)\in \mathop{\mathrm{supp}}\{A_t,H_t\}$ then
  $(\d(a_t,h_t,\varepsilon_t),h_t)\in \mathop{\mathrm{supp}}\{A_t,H_t\}$ for $t\in\{1,\ldots,\tau\}$.
\end{assumptioniden}

\begin{assumptioniden}[Positivity of loss-to-follow-up mechanism]\label{ass:support2}
  $\P\{\mathsf{g}_{C,t}(A_t, H_t)>0\} = 1$ for $t\in\{1,\ldots,\tau\}$.
\end{assumptioniden}

\begin{assumptioniden}[Strong sequential randomization of exposure mechanism]\label{ass:exch1}
  $U_{A,t}\mbox{$\perp\!\!\!\perp$}} \newcommand{\rs}{R (\underline U_{D, {t+1}}, \underline U_{Y,
    {t+1}},\allowbreak \underline U_{L, {t+1}}, \underline U_{A, t+1})
  \mid (H_t, R_t=C_t=1)$ for all $t\in\{1,\ldots,\tau\}$.
\end{assumptioniden}

\begin{assumptioniden}[Sequential randomization of loss-to-follow-up
  mechanism]\label{ass:exch2}
  $U_{C,t}\mbox{$\perp\!\!\!\perp$}} \newcommand{\rs}{R (\underline U_{D, {t+1}}, \underline U_{Y,
    t+1},\allowbreak \underline U_{L, t+1}) \mid (A_t,H_t, R_t=1)$ for
  all $t\in\{1,\ldots,\tau\}$.
\end{assumptioniden}

Assumptions~\ref{ass:support1} and~\ref{ass:support2} are required so that the
interventions considered are supported in the data. In words,
Assumption~\ref{ass:support1} requires that if there is an individual with
treatment level $a_t$ and covariate values $h_t$ who is event-free and
uncensored at time $t$, there must also be an individual with treatment level
$\d(a_t,h_t,\varepsilon_t)$ and covariate values $h_t$ who is event-free and
uncensored at time $t$. Assumptions~\ref{ass:support2} states that for every
possible value of the covariates $h_t$ that has positive probability of
occurring, there are uncensored individuals.
Assumption~\ref{ass:exch1} is required to identify the effect of an LMTP. It
states that, among event-free uncensored individuals, $H_t$ contains all the
common causes of exposure at time $t$ and future events, exposures, and
covariates. Assumption~\ref{ass:exch2} is required to identify the effect of
a static intervention that prevents loss to follow-up. It states that, among
event-free uncensored individuals, $H_t$ contains all the common causes of
censoring at time $t$ and future events and covariates.

\begin{proposition}[Identification of the effect of LMTPs with
  time-to-event outcomes subject to competing events]\label{theo:iden}
  Set $\mathsf{q}_{\tau+1}= Y_{\tau+1}$ and $R_{\tau+1}=1$. For $t=\tau,\ldots,1$,
  recursively define
  \begin{equation}
    \mathsf{q}_t:
    (a_t, h_t)
    \mapsto \mathsf{E}\left[R_{t+1}\times\mathsf{q}_{t+1}(A_{t+1}^\d,
      H_{t+1})\mid C_t= R_t=1, A_t=a_t,
      H_t=h_t\right],\label{eq:seqreg}
  \end{equation}
  and define $\theta = \mathsf{E}[\mathsf{q}_1(A_1^\d, L_1)]$. Under Assumptions
  \ref{ass:support1}-\ref{ass:exch2}, $\P[Y_{\tau+1}(\bar A^\d)=1]$ is identified as
  $\theta$.
\end{proposition}
The identification result above, proved in the supplementary materials, extends the identification result
of~\citet{young2020causal} to LMTPs, and the identification result
of~\cite{diaz2021lmtp} to the case of competing events.

In the next section we present efficiency theory and efficient estimators for
$\theta$, based on an extension of our prior work~\citep{diaz2021lmtp}. These
estimators allow the use of slow-converging data-adaptive regression techniques
to obtain estimators of $\theta$ that converge at parametric rates. The
general approach involves finding a von-Mises-type
approximation~\citep{mises1947asymptotic, vanderVaart98,
robins2009quadratic, Bickel97} for the parameter $\mathsf{q}_{L, t}$, which can be
intuitively understood to be first-order expansions with second-order error
(remainder) terms. As the errors in the expansion are second-order, the
resulting estimator of $\theta$ will be $\sqrt{n}$-consistent and asymptotically
normal as long as the second-order error terms converge to zero at rate
$\sqrt{n}$. This would be satisfied, for example, if all the regression
functions used for estimation converge at rate $n^{1/4}$ (see
\S\ref{sec:estima}). The reader interested in this general theory is encouraged
to consult the works of \citet{buckley1979linear,rubin2007doubly,
Robins00,Robins&Rotnitzky&Zhao94,vanderLaan2003,Bang05, vdl2006targeted,
vanderLaanRose11, vanderLaanRose18, luedtke2017sequential,
rotnitzky2017multiply}.

\section{Efficient estimation via sequential doubly robust
  regression}\label{sec:estima}

As noted by~\citet{diaz2021lmtp}, $\sqrt{n}$-consistent estimation is not
possible for arbitrary functions $\d$ when the exposure is continuous. Thus, we
impose the following assumption, which is a generalization of an assumption
first used by~\citep{Haneuse2013}. In what follows, we
will assume that the function $\d$ does not depend on the distribution of the
observed data $\P$.

\begin{assumption}[Piecewise smooth invertibility]\label{ass:inv}
  Assume that the cardinality of $A_t$ is $p$. For each $h_t$, assume
  the support of $A_t$ conditional on $H_t=h_t$ may be partitioned
  into $p$-dimensional subintervals
  ${\cal I}_{t,j}(h_t):j = 1, \ldots, J_t(h_t)$ such that
  $\d(a_t, h_t, \varepsilon_t)$ is equal to some
  $\d_j(a_t, h_t, \varepsilon_t)$ in ${\cal I}_{t,j}(h_t)$ and
  $\d_j(\cdot,h_t, \varepsilon_t)$ has inverse function
  $\b_j(\cdot, h_t, \varepsilon_t)$ with Jacobian
  $\b_j'(\cdot, h_t, \varepsilon_t)$ with respect to $a_t$.
\end{assumption}

In what follows, we assume that either (i) (i) $A_t$ is a discrete random
variable for all $t$, or (ii) $A_t$ is a continuous random variable and the modified treatment
policy $\d$ satisfies Assumption~\ref{ass:inv}. Next, for continuous
$A_t$ we define
{\footnotesize\begin{equation}\label{eq:gdelta}
    \mathsf{g}_{A,t}^\d(a_t \mid h_t) =
    \sum_{j=1}^{J_t(h_t)} \int \mathds{1}_{t, j} \{\b_j(a_t, h_t, \varepsilon_t),
    h_t\}\times \mathsf{g}_{A,t}\{\b_j(a_t, h_t, \varepsilon_t)\mid h_t\}
    \times\big|\det\{\b_j'(a_t,h_t, \varepsilon_t)\}\big|\dd\P(\varepsilon_t),
  \end{equation}}%
where $\mathds{1}_{t,j} \{u, h_t\} = 1$ if $u \in {\cal I}_{t, j}(h_t)$ and
$\mathds{1}_{t,j} \{u, h_t\} = 0$ otherwise. Under Assumption~\ref{ass:inv},
it is easy to show that the p.d.f.~of $A_t^\d$ conditional on the
history $h_t$ is $\mathsf{g}_{A,t}^\d(a_t \mid h_t)$. In the case of
Example~\ref{ex:mshift}, the post-intervention p.d.f.~becomes
\begin{equation*}
  \mathsf{g}_{A,t}^\d(a_t\mid h_t) = \mathsf{g}_{A,t}(a_t\mid h_t) \mathds{1}\{a_t \delta <
  l_t(h_t)\} + \delta^{-1}\mathsf{g}_{A,t}(a_t\delta^{-1}\mid h_t)
  \mathds{1}\{a_t\geq l_t(h_t)\},
\end{equation*}
which shows that piecewise smoothness is sufficient to handle interventions such
as~\eqref{eq:mshift} which are not smooth in the full range of the exposure.
Assumption~\ref{ass:inv} and expression~\eqref{eq:gdelta} ensure that we can use
the change of variable formula when computing integrals over $\mathcal A_t$ for
continuous exposures. This is useful for studying properties of the parameter
and estimators we propose.

For discrete exposure variables, we let
\begin{equation}
  \mathsf{g}_{A,t}^\d(a_t\mid h_t) = \sum_{s_t\in \mathcal A_t}\int
  \mathds{1}\{\d(s_t, h_t,\varepsilon_t)=a_t\}\mathsf{g}_{A,t}(s_t\mid
  h_t)\dd\P(\varepsilon_t).\label{eq:disc}
\end{equation}
Evaluated using the definition of the incremental propensity score intervention
based on the risk ratio given in Example~\ref{ex:ipsi}, this expression yields
\[
  \mathsf{g}_{A,t}^\d(a_t\mid h_t) = a_t\delta \mathsf{g}_{A,t}(1\mid h_t) + (1-a_t)[1-\delta
    \mathsf{g}_{A,t}(1\mid h_t)],
\]
which clarifies that the risk ratio for receiving the exposure under the
intervention vs the observed data is equal to $\delta$. In what follows, it will
be useful to define the parameter
\[
  \mathsf{w}_t(c_t, a_t, h_t) = \frac{\mathsf{g}_{A,t}^\d(a_t\mid
    h_t)}{\mathsf{g}_{A,t}(a_t\mid h_t)}\frac{c_t\times r_t}{\mathsf{g}_{C,t}(a_t, h_t)},
\]
and the function
\begin{equation}
  \label{eq:eifgamma}
  \varphi_t : x\mapsto
  \sum_{s=t}^\tau\left(\prod_{k=t}^s \mathsf{w}_k(c_k, a_k, h_k)\right)
  \{r_{s+1}\mathsf{q}_{s+1}(a_{s+1}^\d, h_{s+1}) -
  \mathsf{q}_s(a_s, h_s)\} + \mathsf{q}_t(a_t^\d, h_t)
\end{equation}
for $t=\tau,\ldots,1$. When necessary, we use the notation $\varphi_t(z;\eta)$ or
$\varphi_t(z;\underline\eta_t)$ to highlight the dependence of $\varphi_t$ on
$\underline\eta_t=(\mathsf{w}_t,\mathsf{m}_t,\ldots,\mathsf{w}_\tau, \mathsf{m}_\tau)
. We also use $\eta$ to denote $(\mathsf{w}_1,\mathsf{m}_1,\ldots,\mathsf{w}_\tau, \mathsf{m}_\tau)$, and define
$\varphi_{\tau+1}(Z;\eta) = Y$.

Estimation of $\theta$ will proceed by regression of the transformation
$\varphi_{t+1}$ sequentially from $t=\tau$ to $0=1$ to obtain estimates $\mathsf{m}_t$. We
say that $\varphi_{t+1}$ is a multiply robust unbiased transformation for $\mathsf{m}_t$ due
to the following proposition:

\begin{proposition}[Doubly robust transformation]\label{prop:vonm}
  Let $\eta'$ be such that either $\mathsf{m}'_s=\mathsf{m}_s$ or $\mathsf{w}'_s=\mathsf{w}_s$ for all
  $s>t$. Then, in the event $R_t=1$, we have
  \[\mathsf{E}\big[R_{t+1}\times \varphi_{t+1}(Z;\eta')\mid
    C_t=R_t=1,A_t=a_t,H_t=h_t\big] =\mathsf{q}_t(a_t,h_t).\]
\end{proposition}
This proposition is a straightforward application of the von-Mises expansion
given in Lemma~\ref{supp:lemma:vm} in the supplementary materials. Another
consequence of Lemma~\ref{supp:lemma:vm} is that the efficient influence
function for estimating $\theta = \mathsf{E}[\mathsf{m}_1(A^\d,L_1)]$ in the non-parametric
model is given by $\varphi_1(Z)-\theta$~\citep[see][]{diaz2021lmtp}.
Proposition~\ref{prop:vonm} motivates the construction of the sequential
regression estimator by iteratively regressing an estimate of the data
transformation $\varphi_{t+1}(Z;\eta)$ on $(A_t,H_t)$ among individuals with
$C_t=R_t=1$, starting at $\varphi_{\tau+1}(Z;\eta)=Y$.

In order to avoid imposing entropy conditions on the initial estimators, we use
sample splitting and cross-fitting~\citep{klaassen1987consistent,
zheng2011cross, chernozhukov2018double}. Let ${\cal V}_1, \ldots, {\cal V}_J$
denote a random partition of the index set $\{1, \ldots, n\}$ into $J$
prediction sets of approximately the same size. That is, ${\cal V}_j\subset
\{1, \ldots, n\}$; $\bigcup_{j=1}^J {\cal V}_j = \{1, \ldots, n\}$; and ${\cal
V}_j\cap {\cal V}_{j'} = \emptyset$. In addition, for each $j$, the associated
training sample is given by ${\cal T}_j = \{1, \ldots, n\} \setminus {\cal
V}_j$. We let $\hat \eta_{j}$ denote the estimator of $\eta$ obtained by
training the corresponding prediction algorithm using only data in the sample
${\cal T}_j$. Further, we let $j(i)$ denote the index of the validation set
containing unit $i$. For preliminary cross-fitted estimates
$\hat\mathsf{w}_{1,j(i)},\ldots,\hat \mathsf{w}_{\tau, j(i)}$, the estimator is defined as
follows:
\begin{enumerate}[label=Step \arabic*:, align=left, leftmargin=*]
\item Initialize  $\varphi_{\tau+1}(Z_i;\underline{\check\eta}_{\tau,j(i)})=
  Y_i$ for $i=1,\ldots,n$.
\item For $t=\tau,\ldots,1$:
  \begin{itemize}
  \item Compute the pseudo-outcome
    $\check Y_{t+1,i} = R_{t+1}\times
    \varphi_{t+1}(Z_i;\underline{\check\eta}_{t,j(i)})$ for all
    $i=1,\ldots,n$.
  \item For $j=1,\ldots,J$:
    \begin{itemize}
    \item Regress $\check Y_{t+1,i}$ on $(A_{t,i,}H_{t,i})$ using any
      regression technique and using only data points
      $i\in \mathcal T_{j}$.\label{step:2}
    \item Let $\check \mathsf{q}_{t,j}$ denote the output, update
      $\underline{\check\eta}_{t, j} = (\hat\mathsf{w}_{t,j}, \check
      \mathsf{m}_{t,j},\ldots,\hat\mathsf{w}_{\tau,j}, \check \mathsf{m}_{\tau,j})$, and
      iterate.
    \end{itemize}
  \end{itemize}
\item Define the sequential regression estimator as
  \[\hat\theta=\frac{1}{n}\sum_{i=1}^n\varphi_1(Z_i;\check\eta_{j(i)}).\]
\end{enumerate}
Implementation of the above estimator requires an algorithm for estimation of
the parameters $\mathsf{w}_t$. This parameter may be estimated component-wise as
follows. First, the censoring mechanism $\mathsf{g}_{C,t}$ may be estimated through any
regression or classification procedure by regressing the indicator $C_t$ on data
$(A_t,H_t)$ among individuals with $R_t=C_t=1$. Thus, it only remains to
estimate the density ratio $\mathsf{g}_t^\d(a_t\mid h_t)/\mathsf{g}_t(a_t\mid h_t)$. An option
to estimate this density ratio is to first obtain an estimate of the density
$\mathsf{g}_t(a_t\mid h_t)$, and then to obtain an estimator of $\mathsf{g}_t^\d(a_t\mid h_t)$
by applying formula~\eqref{eq:gdelta} or~\eqref{eq:disc}. This approach is
highly problematic in the case of continuous or multivariate exposures due both
to the curse of dimensionality and to the fact that flexible and computationally
efficient estimation of conditional densities is an under-developed area of
machine learning (relative to regression). This approach can also be problematic
for discrete exposures of high cardinality. Since the conditional densities
themselves are unnecessary for estimation, and direct estimation of the ratio
$\mathsf{g}_t^\d(a_t\mid h_t)/\mathsf{g}_t(a_t\mid h_t)$ suffices, we propose to use
\citet{diaz2021lmtp}'s method for direct estimation of density ratios. Briefly,
this approach starts by constructing an augmented artificial dataset of size
$2n$ in which each observation is duplicated. For each observation, one
duplicate gets assigned the post-intervention exposure value $\d(A_t, H_t)$ and the
other gets the actual observed value $A_t$. Then, the task of density ratio
estimation reduces to the task of classifying which record belongs to the
intervened-upon unit, based on data on the exposure and the history $H_t$.
Specifically, \citet{diaz2021lmtp} show that the odds of this classification
problem is equal to the density ratio.

The above estimation algorithm requires the implementation of several regression
procedures. Many algorithms from the machine learning literature are available
for such regression problems, including those based on ensembles of regression
trees, splines, etc. As stated in the asymptotic normality result below, it is
important that the regression procedures used approximate as well as possible
the relevant components of the true data-generating mechanism. Increasing the
predictive performance of these regression fits can entail estimator selection
from a large list of candidate estimators. If the sample size is large enough,
a principled solution to this problem is the use of Super
Learning~\citep{vanderLaanDudoitvanderVaart06, vanderLaanPolleyHubbard07}, a
form of model stacking~\citep{Wolpert1992,Breiman1996} that constructs an
ensemble model as a convex combination of candidate algorithms from
a user-supplied regression library, where the weights in the combination are
computed based on cross-validation and are such that they minimize the
prediction error of the resulting estimator.

The asymptotic distribution of $\theta$ is Gaussian with variance equal to the
efficiency bound, under conditions. To state the conditions, we first define the
data-dependent parameter
\[
  \check\mathsf{m}^\dagger(a_t, h_t) = \mathsf{E}\big[\varphi_{t+1}(Z;\underline{\check\eta}_t)\mid
    A_t=a_t,H_t=h_t\big],
\]
where the outer expectation is only with respect to the distribution $\P$ of $Z$
(i.e., $\check\eta$ is fixed). The following proposition is a consequence of the
results of~\citet{diaz2021lmtp}:
\begin{proposition}[Consistency and weak convergence of SDR
  estimator]\label{theo:asrem}
  We have
  \begin{enumerate}[label=(\roman*)]
  \item If, for each time $t$, either $||\hat\mathsf{w}_t - \mathsf{w}_t|| = o_\P(1)$
    or $||\check\mathsf{q}_t - \check\mathsf{q}^\dagger_t|| = o_\P(1)$, then
    $\hat\theta =\theta + o_\P(1)$.
  \item If
    $\sum_{t=1}^\tau||\hat\mathsf{w}_t - \mathsf{w}_t||\, ||\check\mathsf{q}_t -
    \check\mathsf{q}_t^\dagger|| = o_\P(n^{-1/2})$ and
    $\P\{\mathsf{w}_t(A_t,H_t) < c\}=\P\{\hat\mathsf{w}_t(A_t,H_t) < c\}=1$ for some
    $c<\infty$, then
    \[\sqrt{n}(\hat\theta - \theta) \rightsquigarrow
      N(0,\sigma^2),\] where
    $\sigma^2=\var\{\varphi_1(Z;\eta)\}$ is the non-parametric efficiency bound.
  \end{enumerate}
\end{proposition}
Note that this means that $\hat\theta$ is \textit{sequentially doubly
  robust} in the sense that, at each time point, it is robust to
misspecification of at most one of the two regression procedures
used. This property is also known in the literature under the name of
$2^\tau$-multiple robustness. We prefer the term \textit{sequentially
  doubly robust} for two reasons. First, contrasted to the double
robustness property which requires one out of two nuisance estimators
to be consistent, the term $2^\tau$-multiple robustness erroneously
conveys the message that consistency would follow from one out of
$2^\tau$ nuisance estimators being consistent. Second, retaining the
term ``double'' in this description conveys a fundamental property of
the estimation procedure: it is based on an expansion (given in
Lemma~\ref{supp:lemma:vm} in the supplementary material) with
second-order remainder terms. In this nomenclature, a triply robust
estimating procedure would be one based on an expansion with
third-order remainder terms, and so forth.\footnote{Some of these
  points were brought to our attention by Edward Kennedy in a recent
  conversation on Twitter (see this thread,
  \url{https://twitter.com/edwardhkennedy/status/1446952129218949128?s=20&t=8o_-20szw0PPfqocYiw7Dg}).}

Under the conditions of the above theorem, a Wald-type procedure, in which the
standard error may be estimated based on the empirical variance of
$\varphi_1(Z;\check\eta_{j(i)})$, yields asymptotically valid confidence intervals
and hypothesis tests. A drawback of the above estimation approach is that there
is no guarantee that the parameter will remain within the bounds of the
parameter space. That is, it is possible, in principle, that the sequential
regression estimator will provide time-to-event probabilities outside of the
closed unit interval $[0,1]$. Furthermore, when the survival function must be
estimated at several points in time, there is no guarantee that the resulting
point estimates will be antitonic (i.e., monotonic decreasing). In the
supplementary materials, we present a targeted minimum loss estimator (TMLE)
capable of producing estimates that always lie within the unit interval.
However, boundedness of the TMLE comes with a robustness tradeoff. In
particular, the robustness properties of the TMLE are inferior to those of the
sequentially doubly robust estimator. Specifically, consistency of the TMLE
requires that either $\lVert \hat\mathsf{w}_t - \mathsf{w}_t \rVert = o_\P(1)$ or $\lVert
\hat\mathsf{q}_t - \hat\mathsf{q}_t \rVert = o_\P(1)$, for an estimator $\hat\mathsf{q}_t$ constructed
based on sequential regression using~\eqref{eq:seqreg}. As a result, consistent
estimation of $\mathsf{q}_t$ in the TMLE procedure requires consistent estimation of
$\mathsf{q}_s: s>t$. By contrast, the sequentially doubly robust estimator provides
consistent estimation of $\mathsf{q}_t$ when either $\mathsf{q}_s$ or $\mathsf{w}_s$ is consistently
estimated for all $s>t$. Both~\citet{luedtke2017sequential} and
\citet{diaz2021lmtp} provide further in-depth discussion on this topic.

To resolve the above issues, we leverage a set of techniques proposed
by~\citet{westling2020correcting}. Under their approach, out-of-bounds estimates
are truncated to remain within bounds, and the possibly non-monotonic estimate
of the survival curve, alongside any confidence region limits, is projected onto
the space of monotonic functions by way of isotonic regression. Importantly,
\citet{westling2020correcting} show that the resulting projected estimator
retains the same asymptotic properties as the original estimator. In our case,
this implies that the projected estimator retains the properties of asymptotic
normality and consistency under the assumptions of Proposition~\ref{theo:asrem}.








\section{Illustrative application}

\subsection{Background on motivating study and data}

To illustrate the proposed methods, we aim to answer the following question:
\textit{what is the effect of invasive mechanical ventilation (IMV) on
acute kidney injury (AKI) among COVID-19 patients?}  This question has been
recently identified by an expert panel on lung--kidney interactions in critically
ill patients~\citep{joannidis2020lung} as an area in need of further research.
AKI is a common condition in the ICU, complicating about 30\% of ICU admissions,
and causing increased risk of in-hospital mortality and long-term morbidity and
mortality~\citep{kes2010acute}.
AKI is often viewed as a tolerable consequence of interventions to support other
failing organs, such as IMV~\citep{husain2016lung}.

To answer this question, we will use data on approximately 3,300 patients
hospitalized with COVID-19 at the New York Presbyterian Cornell, Queens, or
Lower Manhattan hospitals. Our analysis focuses on patients hospitalized between
March 3 and May 15, 2020, who did not have a previous history of Chronic Kidney
Disease (CKD). This choice of time frame is rooted in there being a lack of
treatment guidelines in the initial months of the COVID-19 pandemic, which
resulted in a large variability in provider practice regarding mechanical
ventilation. This clinical heterogeneity in the deployment of IMV makes the
problem particularly well-suited to assessment via the causal effects of
LMTPs. The analytical dataset was built from two distinct databases.
Demographics, comorbidity, intubation, death, and discharge data were abstracted
from electronic health records by trained medical professionals into a secure
RedCAP database~\citep{goyal2020clinical}. These data were supplemented with the
Weill Cornell Critical carE Database for Advanced Research (WC-CEDAR),
a comprehensive data repository housing laboratory, procedure, diagnosis,
medication, microbiology, and flowsheet data documented as part of standard
care~\citep{schenck2021critical}. This dataset is unique in that it contains
a large sample of patients treated in the same hospital system within a short
period of time, while still exhibiting extensive treatment variability.

The most frequent reason for morbidity and mortality among patients with
COVID-19 is pneumonia leading to acute respiratory distress syndrome
(ARDS). Respiratory support via devices such as nasal cannulae, face
masks, or endotracheal tubes (i.e., intubation and subsequent
mechanical ventilation) is a primary treatment for
ARDS~\citep{hasan2020mortality}. For severe ARDS patients, an
important decision must be made regarding the optimal timing of
intubation based on clinical factors such as oxygen saturation,
dyspnea, respiratory rate, chest radiograph, and, occasionally,
external resource considerations such as ventilator or provider
availability~\citep{tobinresp, thomson2021timing}. At the outset of
the COVID-19 pandemic, multiple international guidelines recommended
early intubation in an effort to protect healthcare workers from
infection and to avoid complications stemming from ``crash''
intubations~\citep{papoutsi2021effect}. Yet, as the pandemic
progressed and multiple reports documented high mortality for
mechanically ventilated patients, guidance changed to delaying
intubation with the rationale that IMV increases risk of
secondary infections, ventilator-induced lung injury, damage to other
organs, and death~\citep{bavishi2021timing, tobin2006principles}. The
lungs and kidneys are physiologically closely tied, and it has been
hypothesized that mechanical ventilation may cause AKI in COVID-19
patients via oxygen toxicity and capillary endothelial damage leading
to inflammation, hypotension, and
sepsis~\citep{durdevic2020progressive}.

Despite these hypotheses, a causal link between intubation and AKI in patients
with COVID-19 has not been definitively established. Furthermore, such a link is
difficult to study in a randomized fashion due to the clinically subjective
nature of intubation assignment --- that is, it is impossible to pick an oxygen
saturation breakpoint to intubate across all patients and hospital
settings~\citep{tobin2020caution, perkins2020recovery}.

\subsection{Time-varying variables and modified treatment policy}

For each patient, study time begins on the day of hospitalization. The
exposure of interest is daily maximum respiratory support, which can
take three levels: no supplemental oxygen, supplemental oxygen not
including IMV, and IMV. Our goal is to estimate the overall causal
effect on daily AKI rates of an intervention that would delay intubation for
IMV by a single day among 449 patients who received IMV in the first
14 days of hospitalization, where death by other causes is a competing
risk for AKI development. In notation, this LMTP may be expressed as
follows. Let $A_t\in\{0,1,2\}$, where $0$ indicates no supplemental
oxygen at end of day $t$, $1$ indicates supplemental oxygen not
including IMV, and $2$ indicates IMV. Consider the following
intervention:
\begin{equation}\label{eq:exmtp}
\d_t(A_t,H_t) =
  \begin{cases}
    1 &\text{ if } A_t=2 \text{ and } A_s \leq 1 \text{ for all } s < t,\\
    A_t & \text{ otherwise.}
  \end{cases}
\end{equation}
Under this intervention, a patient who is first intubated at day $t$
would receive non-invasive oxygen support at day $t$. Otherwise, the
patient would receive the oxygen support actually (i.e.,
``naturally'') observed at day $t$. This intervention assesses what
would have happened in a hypothetical world where intubation
procedures were more conservative than they were at the beginning of
the pandemic in New York City's surge conditions, when limited
information and treatments were available for COVID-19.
Figure~\ref{fig:os} shows the observed
supplemental oxygen levels and observed AKI and death rates across days 0-14. 

Baseline confounders include age, sex, race, ethnicity, body mass index (BMI),
hospital location, home oxygen status, and comorbidities (e.g., hypertension,
history of stroke, Diabetes Mellitus, Coronary Artery Disease, active cancer,
cirrhosis, asthma, Chronic Obstructive Pulmonary Disease, Interstitial Lung
Disease, HIV infection, and immunosuppression). Time-dependent confounders
include vital signs (e.g., highest and lowest respiratory rate, oxygen
saturation, temperature, heart rate, and blood pressure), laboratory results
(e.g., C-Reactive Protein, BUN Creatinine Ratio, Activated Partial
Thromboplastin time, Creatinine, lymphocytes, neutrophils, bilirubin, platelets,
D-dimer, glucose, arterial partial pressure of oxygen, and arterial partial
pressure of carbon dioxide), and an indicator for daily corticosteroid
administration greater than 0.5 mg/kg body weight. In cases of missing baseline
confounders, multiple imputation using chained equations
(MICE)~\citep{mice2011jss} is used. For missing data at later time points, the
last observation is carried forward. Patients are censored at their day of
hospital discharge, as AKI and vital status were unknown after this point.

\begin{figure}[!htp]
  \centering
  \includegraphics[scale = 0.6]{figure_alluvial}
  \caption{Alluvial plot showing the movement of patients between exposure and
  outcome statuses across study days. The exposure contains three levels: no
  supplemental oxygen, supplemental oxygen not including IMV, and IMV. The
  outcome of interest is AKI, and death is a competing risk. Patients who are
  discharged are right-censored.}
  \label{fig:os}
\end{figure}

\subsection{Results}

We estimate the effect of the LMTP with the SDR estimator (on account of its
enhanced robustness relative to TMLE), using a version of the \texttt{lmtp}
\texttt{R} package~\citep{williams2022lmtp} modified to support the competing
risks setting. The nuisance functions $\r_t$ and $\mathsf{m}_t$ are estimated using the
Super Learner ensemble modeling algorithm~\citep{vanderLaanPolleyHubbard07} via
the \texttt{sl3} \texttt{R} package~\citep{coyle2022sl3}. The candidate library
considered by the Super Learner included a total of 32 learning algorithms,
inclusive of hyperparameter variations. The algorithms included multivariate
adaptive regression splines~\citep[MARS;][]{friedman1991multivariate};
random forests~\citep{breiman2001random, wright2017ranger}; extreme gradient
boosted trees~\citep{friedman2001greedy, chen2016xgboost}; lasso, ridge, and
elastic net (with equal weights, i.e., $\alpha = 0.5$) regularized
regression~\citep{friedman2009glmnet}; main terms generalized linear
models~\citep{mccullagh1989generalized, enea2009fitting}; generalized linear
models with Bayesian priors on parameter estimates~\citep{gelman2006data}; and
an unadjusted outcome mean.  The estimators of $\r_t$ and $\mathsf{m}_t$ were computed
under a Markov assumption with lag $r=2$ to reflect the collection of laboratory
results at minimum 48-hour intervals.
Figure~\ref{fig:sdr_survest} displays the SDR-based cumulative event curve estimates for
each treatment condition (panel A) after the correction for monotonicity
proposed by~\citet{westling2020correcting}, and the incidence difference between
the two (panel B), together with 95\% simultaneous confidence bands.

\begin{figure}[!htb]
  \centering
  \includegraphics[scale=0.315]{sdr_surv_paneled}
  \caption{Panel A: Estimated AKI incidence across the arms defined
    by the two treatment strategies, 1-day delay-of-intubation and no
    intervention. Panel B: The time-specific difference between
    the two treatment conditions. In both panels, 95\% simultaneous confidence
    bands cover the sets of point estimates.}
  \label{fig:sdr_survest}
\end{figure}
Examination of Figure~\ref{fig:sdr_survest} panel A reveals a slight
but increasingly widening gap (across days) between the SDR-based incidence
estimates for each arm. In agreement with this, inspection of panel B of
Figure~\ref{fig:sdr_survest} makes apparent the protective effect of the
delay-in-intubation policy against AKI. The estimates of the incidence difference
are consistently negative and decrease steadily in time; moreover, the
simultaneous confidence bands do not overlap with the null value of zero except
for on the first day. Hypothesis testing of the incidence difference yields
adjusted p-values $< 0.001$ for all days after the first, indicating
a statistically significant protective effect of the LMTP; the Bonferroni
correction was used to adjust for testing multiplicity. Based on these
estimates, the delay-in-intubation policy would be expected to result in
a decrease in 14-day AKI rates by approximately 3.5\% (simultaneous 95\% CI:
[2.7\%, 4.3\%]) relative to no intervention.



\section{Discussion}

Modified treatment policies represent a recent but powerful advance in modern
causal inference. These intervention regimes have a variety of applications,
including (i) practical solutions to violations of the positivity assumption;
(ii) definition and estimation of the effects of continuous, time-valued, and
multivariate exposures; and (iii) definition of effects that depend on the
natural value of treatment. Here, we have advanced the LMTP methodology,
extending it to the case of survival outcomes with competing events. We have
presented identification assumptions for the causal parameters, as well as
a sequentially doubly robust and efficient estimation algorithm capable of
leveraging the many data-adaptive regression procedures available for nuisance
parameter estimation. We ensure that the estimates of the survival function
remain within the parameter space (antitonic and bounded in the closed unit
interval) by leveraging the monotonic projection framework
of~\citet{westling2020unified} and~\citet{westling2020correcting}. We have
illustrated the scientific utility of our approach through application to the
study of negative side effects of invasive mechanical ventilation on acute
kidney injury in patients hospitalized with COVID-19, for which death by other
causes serves as a competing event.

Our Proposition~\ref{prop:vonm} shows that doubly robust estimation is possible
for any intervention that is defined in terms of an LMTP satisfying the stated
conditions. Namely, any randomized intervention defined by a known function
$\d(a_t,h_t,\varepsilon_t)$ through either an LMTP or a general stochastic
intervention accommodates doubly robust estimation, provided that: (1) either
(i) $A_t$ is a discrete random variable for all $t$, or (ii) $A_t$ is
a continuous random variable and the modified treatment policy $\d$ satisfies
Assumption~\ref{ass:inv}, and (2) the randomizer $\varepsilon_t$ is drawn
independently across units and independently of $U$, and its distribution does
not depend on $\P$. Notably, this rich class contains all interventions
considered by Theorem 2 of~\citet{wen2021intervention}.


\section*{Conflict of Interest}

The authors declare that they have no conflict of interest.


\section{Auxiliary results}
The proofs of results in the paper rely heavily on existing
identification and optimality results for general longitudinal
modified treatment policies~\citep{diaz2021lmtp}. We first
present those results. Consider a data structure
$Z=(W_1, E_1, \ldots,W_{\tau}, E_{\tau}, Y)$. We assume these data are
generated by a non-parametric structural equation model with exogenous
errors $U_{W_t}$ and $U_{E_t}$ as detailed by
\cite{diaz2021lmtp}. Let $\mathsf{m}_{C, \tau+1}=Y$,
$V_t = (\bar W_t, \bar E_{t-1})$, and let $\mathbbm{h}(e_t, v_t)$ denote a
modified treatment policy on $E$, with
$E_t^\mathbbm{h}=\mathbbm{h}(E_t(\bar E_{t-1}^\mathbbm{h}),V_t(\bar E_{t-1}^\mathbbm{h}))$ denoting the
intervened exposure variable. For $t=\tau,\ldots,0$, let
\[\mathsf{m}_t(e_t, v_t)=\mathsf{E}[\mathsf{m}_{t+1}(E_{t+1}^\mathbbm{h}, V_{t+1})\mid E_t= e_t,
  V_t=v_t].\] Furthermore, let $\mathsf{g}_t(e_t\mid v_t)$ denote the
probability density function of $E_t$ conditional on $V_t=v_t$,
$\mathsf{g}_t^\mathbbm{h}(e_t\mid v_t)$ the probability density function of $E_t^\mathbbm{h}$
conditional on $V_t=v_t$, and define
\[\r_t(e_t,v_t) = \frac{\mathsf{g}_t^\mathbbm{h}(e_t\mid v_t)}{\mathsf{g}_t(e_t\mid v_t)}.\]
Define
\[\psi_t(z) = \sum_{s=t}^\tau \left(\prod_{k=t}^s\r_k(e_k, v_k)
  \right)\{\mathsf{m}_{s+1}(e_{s+1}^\mathbbm{h},v_{s+1}) - \mathsf{m}_s(e_s,v_s)\} +
  \mathsf{m}_t(e_t^\mathbbm{h},v_t).\]
Then \cite{diaz2021lmtp} proved the following three
results.
\begin{theorem}\label{theo:gcompid}
  Let $Y(\bar E^\mathbbm{h})$ denote a counterfactual outcome in a hypothetical
  world where $\bar E = \bar E^\mathbbm{h}$. Assume
  $U_{E,t}\mbox{$\perp\!\!\!\perp$}} \newcommand{\rs}{R \ubar U_{W,t+1}\mid V_t$ and that
  $(e_t,v_t)\in\mathop{\mathrm{supp}}\{E_t,V_t\}$ implies $(\mathbbm{h}(e_t,v_t),v_t)\in\mathop{\mathrm{supp}}\{E_t,V_t\}$. Then
  $\lambda=\mathsf{E}[\mathsf{m}_1(E_1^\mathbbm{h},V_1)]$ identifies $\mathsf{E}\{Y(\bar E^\mathbbm{h})\}$.
\end{theorem}

\begin{theorem}\label{theo:eifgcomp}
  The efficient influence function for estimating $\lambda$ in the
  non-parametric model is given by $\psi_1(Z)-\lambda$.
\end{theorem}

\begin{theorem}\label{theo:dr} For any
  fixed sequence
  $(\mathsf{m}_{1}', \r_{1}', \ldots, \mathsf{m}_{\tau}', \r_{\tau}')$,
  define
  \[\mathsf{K}_{t,s}' =
    \prod_{l=t+1}^{s-1}\r_l'(E_l, V_l),\]
  and
  \[\mathsf{U}_t(e_t, v_t) =
    \sum_{s=t+1}^\tau\mathsf{E}\left[\mathsf{K}_{t,s}'\{\r_s'(E_s,V_s) - \r_s(E_s,V_s)\}\{\mathsf{m}_s'(E_s,V_s) - \mathsf{m}_s(E_s,V_s)\}\mid
    E_t = e_t, V_t=v_t\right].\]
  Then we have
  \[\mathsf{m}_t(e_t, v_t) = \mathsf{E}[\psi_{t+1}'(Z)\mid E_t=e_t, V_t=v_t] +
    \mathsf{U}_t(e_t, v_t)\] where $\psi_{t+1}'$ is defined as above with
  $\mathsf{m}_s$ and $\r_s$ replaced by $\mathsf{m}_s'$ and $\r_s'$, respectively.
\end{theorem}

\section{Identification (Proposition \ref{main:theo:iden})}
\begin{proof}
  Let $V_t=(D_t, Y_t, L_t)$ and $E_t=(A_t, C_t)$. Let
  $\mathbbm{h}(E_t, V_t) = (\d(A_t,H_t), 1)$. This theorem follows from
  application of Theorem~\ref{theo:gcompid} above. First, note that
  for $t=\tau$, we have
  \begin{align*}
    \mathsf{m}_t(e_t^\mathbbm{h},v_t)
    &= \mathsf{E}[Y_{\tau+1}\mid D_t=d_t, Y_t=y_t, H_t=h_t,
      A_t=a_t^\d, C_t=1]\\
    &= \mathds{1}\{d_t=1,y_t=1\}\mathsf{E}[Y_{\tau+1}\mid D_t= Y_t=C_t=1, H_t=h_t,
      A_t=a_t^\d]\\
    &= r_t\mathsf{q}_t(a_t^\d, h_t).
  \end{align*}
  For $t<\tau$, we have
  \begin{align*}
    \mathsf{m}_t(e_t^\mathbbm{h},v_t)
    &= \mathsf{E}[\mathsf{m}_{t+1}(E_{t+1}^\mathbbm{h},V_{t+1})\mid D_t=d_t, Y_t=y_t, H_t=h_t,
      A_t=a_t^\d, C_t=1]\\
    &= \mathsf{E}[R_{t+1}\mathsf{q}_{t+1}(A_{t+1}^\d, H_{t+1})\mid  D_t=d_t, Y_t=y_t, H_t=h_t,
      A_t=a_t^\d, C_t=1]\\
    &= r_t\mathsf{E}[R_{t+1}\mathsf{q}_{t+1}(A_{t+1}^\d, H_{t+1})\mid D_t= Y_t=C_t=1, H_t=h_t,
      A_t=a_t^\d]\\
    &= r_t\mathsf{q}_t(a_t^\d, h_t).
  \end{align*}
  Since $R_1=1$, we have $\mathsf{m}_1(E_1^\mathbbm{h},V_1)=\mathsf{q}_1(A_1^\d, H_1)$,
  concluding the proof of the result.
\end{proof}

\section{Efficient influence function}
\begin{proof}
  As before, let $V_t=(D_t, Y_t, L_t)$ and $E_t=(A_t, C_t)$. Let
  $\mathbbm{h}(E_t, V_t) = (\d(A_t,H_t), 1)$. First, notice that
  \begin{align*}
    \r_t(e_t, v_t) & = \frac{\mathsf{g}_t^\mathbbm{h}(e_t\mid v_t)}{\mathsf{g}_t(e_t\mid
                     v_t)}\\
                   & = \frac{\mathsf{g}_{A,t}^\d(a_t\mid h_t)}{\mathsf{g}_{A,t}(a_t\mid h_t)}\frac{c_t}{\mathsf{g}_{C,t}(a_t,h_t)}.
  \end{align*}
  Then we have
  \begin{align}
    \psi_t(z)
    &= \sum_{s=t}^\tau \left(\prod_{k=t}^s\r_k(e_k,
      v_k)\right)\{\mathsf{m}_{s+1}(e_{s+1}^\mathbbm{h},v_{s+1}) -
      \mathsf{m}_s(e_s,v_s)\} + \mathsf{m}_t(e_t^\mathbbm{h},v_t)\notag\\
    &= \sum_{s=t}^\tau \left(\prod_{k=t}^s\r_k(e_k,
      v_k)\right)\{r_{s+1}\mathsf{q}_{s+1}(a_{s+1}^\mathbbm{h},h_{s+1}) -
      r_s\mathsf{q}_s(a_s,h_s)\} + r_t\mathsf{q}_t(a_t^\d,h_t)\notag\\
    &= \sum_{s=t}^\tau \left(\prod_{k=t}^s\r_k(e_k,
      v_k)\right)r_s\{r_{s+1}\mathsf{q}_{s+1}(a_{s+1}^\mathbbm{h},h_{s+1}) -
      \mathsf{q}_s(a_s,h_s)\} + r_t\mathsf{q}_t(a_t^\d,h_t)\notag\\
    &= \sum_{s=t}^\tau \left(\prod_{k=t}^s\r_k(e_k,
      v_k)r_k\right)\{r_{s+1}\mathsf{q}_{s+1}(a_{s+1}^\mathbbm{h},h_{s+1}) -
      \mathsf{q}_s(a_s,h_s)\} + r_t\mathsf{q}_t(a_t^\d,h_t)\notag\\
    &= r_t\left[\sum_{s=t}^\tau \left(\prod_{k=t}^s\mathsf{w}_t(c_k,a_k,
      h_k)\right)\{r_{s+1}\mathsf{q}_{s+1}(a_{s+1}^\mathbbm{h},h_{s+1}) -
      \mathsf{q}_s(a_s,h_s)\} + \mathsf{q}_t(a_t^\d,h_t)\right]\notag\\
    &=r_t\varphi_t(x)\label{eq:eqeif}.
  \end{align}
  An application of Theorem~\ref{theo:eifgcomp} together with the
  above results for $t=1$ yields the desired result after noticing
  that $R_1=1$.
\end{proof}

\section{von-Mises expansion}
\begin{lemma}[von-Mises expansion]\label{lemma:vm}
  Let
  {\small\[\mathsf{R}_t(a_t,h_t)=\sum_{s=t+1}^\tau\mathsf{E}\left[\H_{t,s}'R_s\{\mathsf{w}_s'(A_s,H_s) -
      \mathsf{w}_s(A_s,H_s)\}\{\mathsf{q}_s'(A_s,H_s) - \mathsf{q}_s(A_s,H_s)\}\mid R_t=C_t=1,
      A_t=a_t, H_t=h_t
    \right].\]}%
  Then, in the event $R_t=1$:
  \[\mathsf{q}_t(a_t, h_t)=\mathsf{E}[R_{t+1}\varphi_{t+1}'(X)\mid R_t= C_t=1, A_t=a_t, H_t=h_t] + \mathsf{R}_t(a_t, h_t).\]
\end{lemma}
\begin{proof}
  This results follows from an application of Theorem~\ref{theo:dr} to
  a data structure $V_t=(D_t, Y_t, L_t)$ and $E_t=(A_t, C_t)$ with
  $\mathbbm{h}(E_t, V_t) = (\d(A_t,H_t), 1)$, together with expression
  (\ref{eq:eqeif}). First, note that
  {\scriptsize
    \begin{align*}
      \mathsf{U}_t(e_t, v_t)
      & = \sum_{s=t+1}^\tau\mathsf{E}\left[\mathsf{K}_{t,s}'\{\r_s'(E_s,V_s) - \r_s(E_s,V_s)\}\{\mathsf{m}_s'(E_s,V_s) - \mathsf{m}_s(E_s,V_s)\}\mid
        E_t = e_t, V_t=v_t\right]\\
      & = \sum_{s=t+1}^\tau\mathsf{E}\left[\H_{t,s}'R_s\{\mathsf{w}_s'(A_s,H_s) -
        \mathsf{w}_s(A_s,H_s)\}\{\mathsf{q}_s'(A_s,H_s) - \mathsf{q}_s(A_s,H_s)\}\mid D_t=d_t, Y_t=y_t, C_t=c_t,
        A_t=a_t, H_t=h_t
        \right]\\
      & = \sum_{s=t+1}^\tau\mathsf{E}\left[\H_{t,s}'R_s\{\mathsf{w}_s'(A_s,H_s) -
        \mathsf{w}_s(A_s,H_s)\}\{\mathsf{q}_s'(A_s,H_s) - \mathsf{q}_s(A_s,H_s)\}\mid D_t=d_t, Y_t=y_t, C_t=c_t,
        A_t=a_t, H_t=h_t
        \right]\\
      &= r_t\sum_{s=t+1}^\tau\mathsf{E}\left[\H_{t,s}'R_s\{\mathsf{w}_s'(A_s,H_s) -
        \mathsf{w}_s(A_s,H_s)\}\{\mathsf{q}_s'(A_s,H_s) - \mathsf{q}_s(A_s,H_s)\}\mid R_t=C_t=1,
        A_t=a_t, H_t=h_t
        \right]\\
      &= r_t\mathsf{R}_t(a_t,h_t)\\
  \end{align*}}%
Specifically, we have
  \begin{align*}
    r_t\mathsf{q}_t(a_t, h_t) &= \mathsf{m}_t(e_t, h_t)\\
    &= \mathsf{E}[\psi_{t+1}'(Z)\mid E_t=e_t, V_t=v_t] + \mathsf{U}_t(e_t, v_t)\\
    &= \mathsf{E}[R_{t+1}\varphi_{t+1}'(X)\mid D_t=d_t, Y_t=y_t, C_t=c_t,
      A_t=a_t, H_t=h_t] + \mathsf{U}_t(e_t, v_t)\\
    &= r_t\mathsf{E}[R_{t+1}\varphi_{t+1}'(X)\mid R_t= C_t=1, A_t=a_t, H_t=h_t] + r_t\mathsf{R}_t(a_t, h_t),\\
  \end{align*}
  which concludes the proof.
\end{proof}





\section{Targeted Minimum Loss Based Estimator}

The targeted minimum loss-based estimator is computed as a
substitution estimator that uses an estimate $\tilde\mathsf{q}_{1,j(i)}$
carefully constructed to solve the cross-validated efficient influence
function estimating equation
$\mathsf{P}_n \{\varphi_1(\cdot, \tilde\eta_{j(\cdot)}) - \hat\theta_{\mbox{\scriptsize tmle}}\}= 0$. The
construction of $\tilde\mathsf{q}_{1,j(i)}$ is motivated by the observation
that the efficient influence function of $\theta$ can be expressed as
a sum of terms of the form:
\[\left(\prod_{k=1}^t \mathsf{w}_k(c_k,a_k, h_k)\right)\{r_{t+1}\mathsf{q}_{t+1}(a_{t+1}^\d, h_{t+1}) -
  \mathsf{q}_t(a_t, h_t)\},\] which take the form of score functions
$\lambda(W)\{M - \mathsf{E}(M\mid W)\}$ for appropriately defined variables $M$
and $W$ and some weight function $\lambda$. It is well known that if
$\mathsf{E}(M\mid W)$ is estimated within a weighted generalized linear model
with canonical link that includes an intercept, then the weighted MLE
estimate solves the score equation
$\sum_i\lambda(W_i)\{M_i - \hat\mathsf{E}(M_i\mid W_i)\}=0$. TMLE uses this
observation to iteratively tilt preliminary estimates of $\mathsf{q}_t$
towards a solution of the efficient influence function estimating
equation. The algorithm is defined as follows:

\begin{enumerate}[label=Step \arabic*., align=left, leftmargin=*]
\item Initialize $\tilde\eta =\hat\eta$ and
  $\tilde\mathsf{q}_{\tau+1,j(i)}(A^\d_{\tau+1,i}, H_{\tau+1,i}) = Y_i$.
\item For $s=1,\ldots,\tau$, compute the weights
  \[\lambda_{s,i} = \prod_{k=1}^s \hat\mathsf{w}_{k,j(i)}(C_{k,i},A_{k,i}, H_{k,i})\]  
\item For $t=\tau,\ldots,1$:
  \begin{itemize}
  \item Compute the pseudo-outcome $\check Y_{t+1,i} = R_{t+1}\tilde\mathsf{q}_{t+1,j(i)}(A^\d_{t+1,i}, H_{t+1,i})$
  \item Fit the generalized linear tilting model
    \[\link \tilde\mathsf{q}_t^\epsilon(A_{t,i},H_{t,i}) = \epsilon + \link
      \tilde\mathsf{q}_{t,j(i)}(A_{t,1}, H_{t,i})\] where $\link(\cdot)$ is
    the canonical link. Here, $\link \tilde\mathsf{q}_{t,j(i)}(a_t,h_t)$ is an
    offset variable (i.e., a variable with known parameter value equal
    to one). The parameter $\epsilon$ may be estimated by running a
    generalized linear model of the pseudo-outcome
    $\check Y_{t+1,i}$ with only an
    intercept term, an offset term equal to
    $\link \tilde\mathsf{q}_{t,j(i)}(A_{t,i},H_{t,i})$, and weights
    $\lambda_{t,i}$, using all the data points in the sample such that
    $R_{i,t}=C_{i,t}=1$. An outcome bounded in an interval $[a,b]$ may
    be analyzed with logistic regression (i.e., $\link=\logit$) by
    mapping it to an outcome $(0,1)$ through the transformation
    $(Y - a) / (b - a)(1-2\gamma) + \gamma$ for some small value
    $\gamma>0$. This approach has robustness advantages compared to
    fitting a linear model as it guarantees that the predictions in
    the next step remain within the outcome space \citep[see
    ][]{Gruber2010t}.
  \item Let $\hat\epsilon$ denote the maximum likelihood estimate, and
    update the estimates as
    \begin{align*}
      \link \tilde\mathsf{q}_{t, j(i)}^{\hat\epsilon}(A_{t,i},H_{t,i})&=
                                                                \hat\epsilon + \link\tilde\mathsf{q}_{t, j(i)}(A_{t,i},H_{t,i})\\
      \link \tilde\mathsf{q}_{t, j(i)}^{\hat\epsilon}(A_{t,i}^\d,H_{t,i})&=
                                                                   \hat\epsilon + \link\tilde\mathsf{q}_{t, j(i)}(A_{t,i}^\d,H_{t,i}).
    \end{align*}
    The above procedure with canonical link guarantees the following
    score equation is solved: {\small
      \[\frac{1}{n}\sum_{i=1}^n \left(\prod_{k=1}^t \hat\mathsf{w}_{k,j(i)}(C_{k,i},A_{k,i},
          H_{k,i})\right)\{R_{t+1,i}\tilde \mathsf{q}_{t+1,j(i)}(A^\d_{t+1,i}, H_{t+1,i}) -
        \tilde\mathsf{q}_{t, j(i)}^{\hat\epsilon}(A_{t,i},H_{t,i})\}=0\]
    }
  \item Update $\tilde\mathsf{q}_{t, j(i)}= \tilde\mathsf{q}_{t, j(i)}^{\hat\epsilon}$, $t = t-1$, and iterate. 
  \end{itemize}
\item The TMLE is defined as
  \[\hat\theta=\frac{1}{n}\sum_{i=1}^n\tilde\mathsf{q}_{1,j(i)}(A_{1,i}^\d,L_{1,i}).\]
\end{enumerate}
The iterative procedure and the score equation argument above
guarantee that
\[\frac{1}{n}\sum_{i=1}^n\sum_{t=1}^\tau\left(\prod_{k=1}^t
    \hat\mathsf{w}_{k,j(i)}(C_{k,i},A_{k,i}, H_{k,i})\right)\{R_{t+1,
    i}\tilde\mathsf{q}_{t+1,j(i)}(A_{t+1,i}^\d, H_{t+1,i}) -
  \tilde\mathsf{q}_{t,j(i)}(A_{t,i}, H_{t,i})\}=0,\] and thus that
$\mathsf{P}_n \{\varphi_1(\cdot;\tilde\eta_{j(\cdot)}) - \hat\theta\}= 0$. This fact is
crucial to prove weak convergence and consistency of the TMLE.


 
 
 
 
 
 
 
 
 
 




