# Intervention treatment distributions that depend on the observed treatment
process and model double robustness in causal survival analysis

Lan Wen${}^{\text{a}}$, Julia L. Marcus${}^{\text{b}}$ and Jessica G.
Young${}^{\text{b}}$
${}^{\text{a}}$Department of Statistics and Actuarial Science, University of
Waterloo
200 University Ave W, Waterloo, ON N2L 3G1
${}^{\text{b}}$Department of Population Medicine, Harvard Medical School,
Boston, Massachusetts, 02215, U.S.A.

###### Abstract.

The generalized g-formula can be used to estimate the probability of survival
under a sustained treatment strategy. When treatment strategies are
deterministic, estimators derived from the so-called efficient influence
function (EIF) for the g-formula will be doubly robust to model
misspecification. In recent years, several practical applications have
motivated estimation of the g-formula under non-deterministic treatment
strategies where treatment assignment at each time point depends on the
observed treatment process. In this case, EIF-based estimators may or may not
be doubly robust. In this paper, we provide sufficient conditions to ensure
existence of doubly robust estimators for intervention treatment distributions
that depend on the observed treatment process for point treatment
interventions, and give a class of intervention treatment distributions
dependent on the observed treatment process that guarantee model doubly and
multiply robust estimators in longitudinal settings. Motivated by an
application to pre-exposure prophylaxis (PrEP) initiation studies, we propose
a new treatment intervention dependent on the observed treatment process. We
show there exist 1) estimators that are doubly and multiply robust to model
misspecification, and 2) estimators that when used with machine learning
algorithms can attain fast convergence rates for our proposed intervention.
Theoretical results are confirmed via simulation studies.
keywords: Causal inference, Double robustness, Estimating equations,
Observational study, Stochastic treatment strategies

##  1\. Introduction

The goal of many observational analyses is to estimate the causal effect on
survival of different time-fixed or time-varying treatment strategies,
interventions or rules in a study population. These causal effects can be
formally defined by a contrast (e.g., difference or ratio) in the
distributions of counterfactual outcomes had interventions been implemented to
ensure those strategies are followed in that population. Robins (1986) [1]
showed that, under assumptions that allow complex longitudinal data structures
such that measured time-varying confounders may themselves be affected by past
treatment, the g-formula indexed by a particular treatment strategy identifies
the average counterfactual outcome under that strategy. Therefore, estimators
of the g-formula and associated contrasts indexed by different strategies may
be used to estimate causal effects.

In practice, the g-formula typically depends on high-dimensional nuisance
parameters. In this case, many estimators of the g-formula and associated
contrasts have been proposed including the density-based parametric g-formula
[1], iterated conditional expectation (ICE) estimators [2, 3], inverse
probability weighted (IPW) estimators [4, 5], and estimators derived from the
efficient influence function (EIF) [6, 7]. EIF based estimators (i.e.,
estimators constructed to evaluate the EIF from an empirical sample) have
several theoretical advantages over the other approaches including they may be
$\sqrt{n}$-consistent if the nuisance functions are estimated at slower rates
through flexible nonparametric or machine learning methods [8, 9, 10].

EIF based estimators may also have a model double-robustness property in that,
when nuisance functions are estimated via parametric models, these estimators
may remain consistent and asymptotically normal if models for only one of two
(sets of) nuisance functions are correctly specified, not necessarily both.
This model double robustness property always holds for EIF estimators when the
g-formula is indexed by a deterministic treatment strategy at most dependent
on past treatment and confounders measured in the observational study [6, 11,
12]. However, the identification results of Robins (1986)[1] were not limited
to such deterministic strategies but generalized to allow identification of
stochastic treatment strategies at most dependent on this measured past. The
latter identifying functional or generalized g-formula depends on the
intervention treatment distribution, that is, the distribution of treatment
under an intervention that ensures the strategy of interest is followed
conditional only on the measured past in the observational study. The
generalized g-formula coincides with the more familiar g-formula indexed by a
deterministic strategy when the intervention treatment distribution is chosen
as degenerate conditional on any level of the measured past.

Recently, several practical applications have motivated estimation of the
generalized g-formula indexed by intervention treatment distributions that
depend on the observed treatment process, that is, the observed treatment
distribution conditional on the measured past [13, 14, 15, 16, 17]. The
generalized g-formula indexed by an intervention treatment distribution
dependent on the observed treatment process has the particular advantage of
relying on relatively weak positivity conditions [15, 18] even, for example,
in observational studies where the propensity score is equal or close to zero
for certain measured confounder histories [16]. When the (degenerate or non-
degenerate) intervention treatment distribution does not depend on the
observed treatment process, EIF derived estimators of the generalized
g-formula will be model doubly robust. However, when the intervention
treatment distribution does depend on the observed treatment process, such
estimators may or may not be doubly robust [14, 15, 16, 19].

In this paper, we exploit particular representations of the generalized
g-formula to give sufficient conditions for the existence of doubly robust
estimators for point treatment interventions when the chosen intervention
treatment distribution depends on the observed treatment process with examples
from the recent literature. We also provide a general form of EIFs for a class
of intervention treatment distributions that may depend on the observed
treatment process in longitudinal settings that guarantee model doubly (and
multiply) robust estimators. Motivated by observational studies of the effects
of realistic HIV pre-exposure prophylaxis (PrEP) initiation interventions, we
consider a new class of intervention treatment distributions dependent on the
observed treatment process that is a variation on the incremental propensity
score interventions proposed by Kennedy (2019)[16]. We show that estimators
based on the EIF for our proposed intervention treatment distribution are
model doubly/multiply robust, and can attain fast convergence rates even when
used in combination with machine learning algorithms, where modelling
assumptions are relaxed. We illustrate both EIF-based, as well as simpler
singly robust, estimators of the g-formula indexed by this class of
intervention treatment distribution in simulated data and in an illustrative
data application.

##  2\. Observed data structure

Consider a longitudinal study with $j=0,1,2,\ldots,J$ denoting a follow-up
interval (e.g., week, month) where $J$ is the end of the follow-up of
interest. Assume the following random variables are measured in this study on
each of $n$ individuals meeting some eligibility criteria at baseline. For
each $j=0,1,2,\ldots,J-1$, let $A_{j}$ denote binary or discrete treatment
variable during interval $j$, $\operatorname*{\textbf{L}}_{j}$ a vector of
additional time-varying covariates measured in interval $j$, and $Y_{j+1}$ an
indicator of survival by interval $j+1$. For notational simplicity, we will
assume throughout that all covariates are discrete in that they have
distributions that are absolutely continuous with respect to a counting
measure but arguments naturally extend to settings with continuous covariates
and Lebesgue measures. By definition, ${Y}_{0}=1$ (all individuals are at risk
of failure baseline) and by convention we define
$\bar{\operatorname*{\textbf{L}}}_{-1}=\bar{A}_{-1}=\emptyset$. For a random
variable $X$, we let $\bar{X}_{j}=(X_{0},\ldots,X_{j})$ denote history through
time $j$. We assume the ordering
$O=(\operatorname*{\textbf{L}}_{0},A_{0},Y_{1},\ldots,\operatorname*{\textbf{L}}_{J-1},A_{J-1},Y_{J})$.
Without loss of generality, we will assume no individual is lost to follow-up
until Section 7.3.

##  3\. Intervention treatment distribution

Let $g$ denote a treatment rule that specifies how treatment should be
assigned at each $j=0,1,2,\ldots,J$. Following Richardson and Robins (2013)
[20], denote $(\operatorname*{\textbf{L}}^{g}_{j},Y^{g}_{j})$ and $A_{j}^{+g}$
as the natural values of covariates and survival status and the intervention
value of treatment at $j$ under $g$, respectively. In turn, the distribution
of $A^{+g}_{j}$ evaluated at some treatment level $a_{j}$ conditional on the
“measured past” under $g$
$(Y^{g}_{j}=1,\overline{\operatorname*{\textbf{L}}}^{g}_{j}=\overline{\operatorname*{\textit{{l}}}}_{j},\overline{A}^{+g}_{j-1}=\overline{a}_{j-1})$
is specified by $q^{g}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})\equiv\Pr(A^{+g}_{j}=a_{j}\mid
Y^{g}_{j}=1,\overline{\operatorname*{\textbf{L}}}^{g}_{j}=\overline{\operatorname*{\textit{{l}}}}_{j},\overline{A}^{+g}_{j-1}=\overline{a}_{j-1})$
which we refer to as the intervention treatment distribution at $j$ associated
with $g$.

When treatment assignment at any time under a selected rule $g$
deterministically depends on the measured past, there is only one value
$a^{+}_{j}\in\text{supp}(A^{+g}_{j})$ given any history
$(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})\in\text{supp}(\overline{\operatorname*{\textbf{L}}}^{g}_{j},\overline{A}^{+g}_{j-1})$
for those with $Y^{g}_{j}=1$. In this case, $q^{g}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=1$ when
$a_{j}=a^{+}_{j}$ and 0 otherwise. Examples include static deterministic rules
that assign the same level of treatment to all surviving individuals at all
follow-up times and dynamic deterministic rules that assign treatment based on
the measured past.

By contrast, when a selected rule $g$ assigns treatment stochastically at some
$j$ (as a random draw from a distribution), at most dependent on the measured
past, then there will be multiple values $a^{+}_{j}\in\text{supp}(A^{+g}_{j})$
such that we may have $0<q^{g}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})<1$ when
$a_{j}=a^{+}_{j}$. We focus here on the problem of estimating
$\mbox{E}[Y_{J}^{g}]=\Pr[Y_{J}^{g}=1]$, the cumulative survival probability by
end of follow-up under a choice of $g$, when the intervention treatment
distribution associated with $g$ has this non-degenerate property of a
stochastic rule, in particular, through its dependence on the observed
treatment distribution conditional on the measured past as specified by
$f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})\equiv\Pr(A_{j}=a_{j}|Y_{j}=1,\overline{\operatorname*{\textbf{L}}}_{j}=\overline{\operatorname*{\textit{{l}}}}_{j},\overline{A}_{j-1}=\overline{a}_{j-1}).$
We refer to $f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})$ as the observed
treatment process conditional on the measured past. The observed treatment
process evaluated at $a_{j}=1$ coincides with the so-called propensity score
[21] at $j$ when treatment $A_{j}$ is binary. Next, we consider several
examples of such intervention distributions.

##  4\. Examples of intervention treatment distributions that depend on the
observed treatment process

In this section, we review examples of intervention distributions considered
in the literature that have depended on the observed treatment process:

  * •

Dynamic treatment initiation strategies with grace period:

Motivated by questions about the effects of CD4-based treatment initiation
strategies, previous authors have considered strategies of the form “If a
condition for treatment initiation is met by interval $j$ then start treatment
by $m+j$ for a selected grace period $m$, with no intervention in intervals
$j$ through $j+m-1$. Otherwise, do not start at $j$”, $\forall j$ [4, 22]. For
$A_{j}$ an indicator of treatment initiation by $j$ and
$L_{j}^{\ast}\in\operatorname*{\textbf{L}}_{j}$ an indicator that the
condition for initiating treatment has been met by $j$, the intervention
treatment distribution at each $j$ is

$q^{g}(1\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})\left.\begin{cases}1,&\text{if
}l^{\ast}_{j-m}=1\\\ 0,&\text{if }l^{\ast}_{j}=0\\\ f(1\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1}),&\text{otherwise}\end{cases}\right.$

or $q^{g}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=(1-l_{j}^{\ast})(1-a_{j})+l_{j-m}^{\ast}a_{j}+(1-l_{j-m}^{\ast})l_{j}^{\ast}f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})$.

  * •

Representative interventions:

Motivated by observational studies to understand the long-term effects of
lifestyle interventions (e.g., interventions that increase daily minutes of
physical activity), previous authors have considered representative
interventions that assign the value of a multi-level treatment to an
individual at each $j=0,\ldots,J$ as a random draw from a particular
distribution: specifically, the observed distribution of treatment in interval
$j$ among those who, in the observational study, (i) had the same measured
confounder and treatment history prior to $j$ as that individual and (ii) had
treatment at $j$ at or above a cutoff $\delta$ (or more generally, within a
pre-specified range), e.g., “at least 30 minutes of daily physical activity”
[23, 17]. In this case,

$q^{g}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1},a_{j}\geq\delta)$

This intervention distribution notably only depends on the observed treatment
process at $j$ among those with treatment in the pre-specified range at $j$.

  * •

Deterministic interventions that depend on the natural value of treatment:

Alternative interventions that maintain a multi-level treatment within a pre-
specified range have been posed that assign treatment at each $j$ as a
function of the natural treatment value at $j$ [24, 13], e.g., “If the natural
value of treatment at $j$ is below $\delta$ than intervene and set treatment
at $j$ to $\delta$. Otherwise, do not intervene at $j$”. The resulting
intervention distribution at each $j$ (conditional only on the measured past
and marginal with respect to the natural value of treatment at $j$) is

$q^{g}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=F_{A_{j}}(\delta\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})I(a_{j}=\delta)+I(a_{j}\geq\delta)f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})$

where $F_{A_{j}}(\delta\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=\sum_{a_{j}<\delta}f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})$.

  * •

Incremental propensity score interventions:  Kennedy (2019)[16] posed
incremental propensity score interventions that at each $j$ assign a binary
treatment according to a strategy $g$ that results in an intervention
treatment distribution defined by a shifted version (on the odds scale) of the
propensity score. Specifically, for a particular $\delta\in(0,\infty)$

$q^{g}(1\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=\dfrac{\delta f(1\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})}{\delta f(1\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})+f(0\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})}$ (1)

or $q^{g}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=\\{{a_{j}\delta
f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})+(1-a_{j})f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})}\\}\\{\delta
f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})+f(1-a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})\\}^{-1}$. Here we
consider a modification of (1) motivated by the PrEP context. Specifically,
for $\delta\in[0,1]$ and $L_{j}^{\ast}\in\operatorname*{\textbf{L}}_{j}$, a
measured marker of risk for HIV acquisition (e.g., receiving a bacterial
sexually transmitted infections or STI test at $j$ and no prior HIV
diagnosis), we consider interventions indexed by the alternative intervention
treatment distribution

$q^{g}(0\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=\left.\begin{cases}f(0\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1}),&\text{if
}l_{j}^{\ast}=0\\\ f(0\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})\delta,&\text{if
}l^{\ast}_{j}=1\end{cases}\right.$ (2)

or $q^{g}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=(1-\delta)l_{j}^{\ast}a_{j}+(l_{j}^{\ast}\delta+1-l_{j}^{\ast})f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})$ after some algebra.

In words, the probability of initiating treatment conditional on the measured
past under $g$ at each $j$ will be larger than the propensity score at $j$
(that under no intervention) by decreasing its complement by a factor of
$\delta$ for those with an indication ($L_{j}^{\ast}=1$). Specifically,
$\delta$ is the risk ratio for not initiating treatment under $g$ vs. no
intervention conditional on the measured past. Choosing $\delta=0$ corresponds
to “always treat” those with $L_{j}^{\ast}=1$, and $\delta=1$ to no
intervention. We will refer to interventions indexed by either (1) or (2) as
incremental propensity score interventions, distinguishing them by the
classifier odds shift or multiplicative shift, respectively.

##  5\. Identification by the generalized g-formula

Consider a treatment assignment rule $g$ at most dependent on the measured
past. Further, let $\mathcal{D}_{g}$ denote the set of all deterministic
strategies at most dependent on this past that individuals could be observed
to follow under the selected rule $g$, with $d$ any element of
$\mathcal{D}_{g}$. In the special case when $g$ is initially selected to be a
deterministic rule then the only element of $\mathcal{D}_{g}$ is $g$.
Otherwise, $\mathcal{D}_{g}$ may contain many elements. Let
$Y_{j}^{d},L_{j}^{d}$ and $A_{j}^{+_{d}}$ denote the natural values of
survival status and covariates and the intervention value of treatment at $j$,
respectively, under a deterministic $d\in\mathcal{D}_{g}$ though
($j=0,\ldots,J)$ and consider the following assumptions:

  1. (1)

Exchangeability: $(Y_{j+1}^{d},\ldots,Y_{J}^{d})\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}A_{j}\mid\bar{\operatorname*{\textbf{L}}}_{j}=\bar{\operatorname*{\textit{{l}}}}_{j},\bar{A}_{j-1}=\bar{a}_{j-1}^{+},Y_{j}=1$.

  2. (2)

Consistency: If $\bar{A}_{j}=\bar{A}_{j}^{+_{d}}$ then
$\bar{Y}_{j+1}=\bar{Y}_{j+1}^{d}$ and
$\bar{\operatorname*{\textbf{L}}}_{j}=\bar{\operatorname*{\textbf{L}}}_{j}^{d}$

  3. (3)

Positivity:
$f_{\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j-1},Y_{j}}(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1}^{+},1)>0\Longrightarrow
f_{A_{j}\mid
Y_{j},\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j-1}}(a_{j}^{+}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1}^{+})>0$

Robins (1986) [1] showed that given these exchangeability, consistency and
positivity conditions hold for all deterministic $d\in\mathcal{D}_{g}$ then
the following function of only the observed data identifies
$\mbox{E}[Y_{J}^{g}]$:

$\displaystyle\psi^{g}=$
$\displaystyle\sum_{\forall\bar{a}_{J-1}}\sum_{\forall\bar{\operatorname*{\textit{{l}}}}_{J-1}}\operatorname*{P}(Y_{J}=1\mid{Y}_{J-1}=1,\bar{\operatorname*{\textbf{L}}}_{J-1}=\bar{\operatorname*{\textit{{l}}}}_{J-1},\bar{A}_{J-1}=\bar{a}_{J-1})\times$
(3)
$\displaystyle\prod_{s=0}^{J-1}\operatorname*{P}(Y_{s}=1\mid{Y}_{s-1}=1,\bar{\operatorname*{\textbf{L}}}_{s-1}=\bar{\operatorname*{\textit{{l}}}}_{s-1},\bar{A}_{s-1}=\bar{a}_{s-1})f({\operatorname*{\textit{{l}}}}_{s}\mid{Y}_{s}=1,\bar{\operatorname*{\textit{{l}}}}_{s-1},\bar{a}_{s-1})q^{\text{\text{g}}}(a_{s}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{s},\bar{a}_{s-1})$

The function $\psi^{g}$ is referred to as the generalized g-formula indexed by
the intervention treatment distribution $q^{\text{\text{g}}}(a_{s}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{s},\bar{a}_{s-1})$. Note that, under
stronger identifying conditions, the generalized g-formula may identify the
outcome mean under a treatment rule $g$ that depends on more than the measured
past [20, 18]. Also see Web Appendix B.

###  5.1. Generalized positivity

Note the assumption that the positivity condition above holds for all
deterministic $d\in\mathcal{D}_{g}$ can be equivalently stated as follows:

$q^{g}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})>0\Longrightarrow
f(a_{j}\mid 1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})>0$ (4)

for all
$\overline{\operatorname*{\textit{{l}}}}_{j},\overline{a}_{j}\in\text{supp}(\overline{\operatorname*{\textbf{L}}}^{g}_{j},\overline{A}^{+g}_{j})$.
The positivity condition (4) generalizes the more familiar definition of
positivity often relied on in the literature that there may be treated and
untreated individuals within any level of the measured past; i.e., the
assumption that the propensity score and its complement are positive for all
possible measured histories and all $j$. It is straightforward to see that the
more general condition (4) reduces to this typical definition of positivity
only for the special case of a static deterministic intervention $g$ on a
binary treatment. By contrast, the more general condition (4) only requires
that, for any level of the past possible in the observational study and also
plausible under $g$, if an intervention level of treatment can occur under $g$
it must also possibly occur in the observational study. Depending on the
choice of $g$, this condition may hold when traditional definitions requiring
positive propensity scores fail. Intervention treatment distributions that
depend on the observed treatment process may help avoid positivity violations
by this more general definition and, in some instances, may guarantee that
positivity violations cannot occur regardless of the observed treatment
process. We discuss this further in the next section.

Similar to arguments given in Kennedy (2019)[16], the odds shift (1) has the
particular advantage that, by construction, the generalized positivity
condition (4) is guaranteed to hold, no matter the nature of the observed
treatment process. By contrast, the multiplicative shift (2) only enjoys this
guarantee for measured pasts consistent with $L_{j}^{\ast}=0$. However,
compared to (1) which is indexed by a shift $\delta$ with no upper bound that
quantifies an odds ratio, (2) may be easier to communicate to subject matter
collaborators as it constrains the choice of $\delta\in[0,1]$ and quantifies a
risk ratio. Notably, the performance of weighted estimators of $\psi^{g}$
indexed by both (1) and (2) are relatively resilient to so-called “near
positivity violations” – such that (4) holds but $f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})$ is still close to
zero for some $(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})$ –
particularly when $\delta$ is chosen to coincide with relatively small
increases in treatment uptake under $g$ (see Section 8).

##  6\. Model double robustness when the intervention treatment distribution
depends on the observed treatment process

Suppose that the observed data $O$ defined in Section 2 follows a law $P$
which is known to belong to
$\mathcal{M}=\\{P_{\theta}\mathrel{\mathop{\mathchar
58\relax}}\theta\in\Theta\\}$ where $\Theta$ is the parameter space. The
efficient influence function (EIF) $U_{\psi^{g}}(O)$ for the causal parameter
$\psi^{g}\equiv\psi^{g}(\theta)$ in a non-parametric model that imposes no
restrictions on the law of $O$ is given by
${d\psi^{g}(\theta_{t})}/{dt}|_{t=0}=E\\{U_{\psi^{g}}(O)S(O)\\}$, where
${d\psi^{g}(\theta_{t})}/{dt}\mid_{t=0}$ is known as the pathwise derivative
of the parameter $\psi^{g}$ along a parametric submodel of the observed data
distribution indexed by $t$, and $S(O)$ is the score function of the
parametric submodel evaluated at $t=0$ [25, 26]. In this section, we provide
results that aid the intuition on the existence of doubly robust estimators of
$\psi^{g}$ when the intervention treatment distribution depends on the
observed treatment process through understanding properties of the EIF for the
parameter $\psi^{g}$.

###  6.1. Point treatment

We begin with the special case of a point treatment where $J=1$ and
$O=(\operatorname*{\textbf{L}}_{0},A_{0},Y_{1})\equiv(\operatorname*{\textbf{L}},A,Y)$.
In this case, (3) reduces to
$\psi^{g}=\sum_{\forall\bar{a}}\sum_{\forall\bar{\operatorname*{\textit{{l}}}}}E(Y|A=a,\operatorname*{\textbf{L}}=\operatorname*{\textit{{l}}})q^{g}(a|\operatorname*{\textit{{l}}})f(\operatorname*{\textit{{l}}})$.

######  Theorem 1.

Suppose $\psi^{g}$ can be written as a linear combination of the form:

$\psi^{g}=c_{1}\underbrace{E\\{h_{1}(O)\\}}_{\nu_{1}}+c_{2}\underbrace{E[E\\{h_{2}(O)\mid
A=a^{\ast},L\\}]}_{\nu_{2}}$ (5)

where $a^{\ast}$, $c_{1}$ and $c_{2}$ are constants, and $h_{1}(O)$ and
$h_{2}(O)$ are known measurable functions of $O$ (i.e., they do not depend on
$\theta$). Then the EIF for $\psi^{g}$ is given by

$U_{\psi^{g}}(O)=c_{1}h_{1}(O)+c_{2}\left[\frac{I(A=a^{\ast})}{f(A\mid\operatorname*{\textbf{L}})}\left[h_{2}(O)-E\\{h_{2}(O)\mid
A,\operatorname*{\textbf{L}}\\}\right]+E\\{h_{2}(O)\mid
A=a^{\ast},\operatorname*{\textbf{L}}\\}\right]-\psi^{g}$ (6)

See Web Appendix A for proof. Clearly, $\psi^{g}$ under a static deterministic
strategy that sets treatment to level $a^{\ast}$ for all individuals trivially
meets the conditions of Theorem 1 by selecting $h_{1}(O)=0$, $c_{2}=1$,
$h_{2}(O)=Y$. In this case, the EIF for the g-formula indexed by $g$ or
$E\\{E(Y\mid A=a^{\ast},\operatorname*{\textbf{L}})\\}$ equals:

$U_{\psi^{g}}(O)=\frac{I(A=a^{\ast})}{f(A\mid\operatorname*{\textbf{L}})}\left\\{Y-m(A,\operatorname*{\textbf{L}})\right\\}+m(a^{\ast},\operatorname*{\textbf{L}})-\psi^{a^{\ast}}$
(7)

where $m(A,\operatorname*{\textbf{L}})\equiv E(Y\mid
A,\operatorname*{\textbf{L}})$ and
$m(a^{\ast},\operatorname*{\textbf{L}})\equiv E(Y\mid
A=a^{\ast},\operatorname*{\textbf{L}})$[6, 27, 7]. A heuristic
justistification for Theorem 6 follows from the fact that the EIF of
${\nu_{1}}\equiv\nu_{1}(\theta)$ is simply $h_{1}(O)$, and the EIF of
${\nu_{2}}\equiv\nu_{2}(\theta)$ can realized by replacing $Y$ with $h_{2}(O)$
in Expression (7), as the function $h_{2}(\cdot)$ does not depend on $\theta$
and therefore its pathwise derivative is zero. Furthermore, it is established
that an estimator derived from the influence function (7) (e.g., an estimator
solving $\sum_{i=1}^{n}U_{\psi^{g}}(O_{i})=0$ for $\psi^{g}$) is model doubly
robust in that it remains consistent if estimated under correctly specified
parametric models for either one of two (sets of) nuisance functions,
specifically $E(Y\mid A,\operatorname*{\textbf{L}})$ or
$f(A\mid\operatorname*{\textbf{L}})$. The following Corollary gives a
sufficient condition for the existence of doubly robust estimators of
$\psi^{g}$ when the intervention treatment distribution depends on the
observed treatment process, provided the conditions of Theorem 1 hold.

######  Corollary 1.1.

Suppose the conditions of Theorem 6 hold. If
$h_{2}(O)=Y\tilde{h}_{2}(A,\operatorname*{\textbf{L}})$, where
$\tilde{h}_{2}(A,\operatorname*{\textbf{L}})$ is a known measurable function
of $(A,\operatorname*{\textbf{L}})$, then an estimator of $\psi^{g}$ derived
from an EIF of the form (6) is model doubly robust.

A proof of Corollary 1.1 is given in Web Appendix A. A similar heuristic
reasoning for Corollary 1.1 is that the estimator of the EIF of a mean outcome
does not rely on any models, and doubly robust estimators exist for $\nu_{2}$
because we have simply replaced $Y$ with $h_{2}(O)$ in Equation (7) which does
not depend on $\theta$. We now consider some applications of Theorem 1 and
Corollary 1.1 to examples where the intervention distribution indexing
$\psi^{g}$ depends on the observed treatment process.

######  Example 1.

Consider a variation of the grace period treatment initiation strategies
defined in Section 4 for $J=1,\leavevmode\nobreak\ m=0$, such that, rather
than withholding treatment when $L^{\ast}=0$, no intervention is made. The
intervention treatment distribution is then given by
$q^{g}(a\mid\operatorname*{\textit{{l}}})=(1-l^{\ast})f(a\mid\operatorname*{\textit{{l}}})+l^{\ast}a$
.

For this choice of intervention distribution we have:

$\displaystyle\psi^{g}$
$\displaystyle=E_{\operatorname*{\textbf{L}}}\left\\{\sum\nolimits_{a=0}^{1}E(Y\mid
a,\operatorname*{\textbf{L}})q^{g}(a\mid\operatorname*{\textbf{L}})\right\\}$
$\displaystyle=E_{\operatorname*{\textbf{L}}}\left[\sum\nolimits_{a=0}^{1}E(Y\mid
a,\operatorname*{\textbf{L}})\left\\{(1-L^{\ast})f(a\mid\operatorname*{\textbf{L}})+I_{\\{1\\}}(a)L^{\ast}\right\\}\right]$
$\displaystyle={E_{\operatorname*{\textbf{L}},A}\left[E\left\\{Y(1-L^{\ast})\mid
A,\operatorname*{\textbf{L}}\right\\}\right]}+{E_{\operatorname*{\textbf{L}}}\left[E\left\\{YL^{\ast}\mid
A=1,\operatorname*{\textbf{L}}\right\\}\right]}$

Selecting $a^{\ast}=1$, $c_{1}=c_{2}=1$, $h_{1}(O)=Y(1-L^{\ast})$, and
$h_{2}(O)=YL^{\ast}$, we have

$\displaystyle\psi^{g}$
$\displaystyle=c_{1}{E\\{h_{1}(O)\\}}+c_{2}{E[E\\{h_{2}(O)\mid
A=a,\operatorname*{\textbf{L}}\\}]}$
$\displaystyle=\underbrace{E\\{Y(1-L^{\ast})\\}}_{\nu_{1}}+\underbrace{E[E\\{YL^{\ast}\mid
A=1,L\\}]}_{\nu_{2}}$

by Theorem 6 and further, the EIF for $\psi^{g}$ is given by

$\displaystyle
U_{\psi^{g}}(O)=Y(1-L^{\ast})+\frac{AL^{\ast}}{f(A\mid\operatorname*{\textbf{L}})}\left\\{Y-m(A,\operatorname*{\textbf{L}})\right\\}+m(1,\operatorname*{\textbf{L}})L^{\ast}-\psi^{g}.$

This can be re-expressed as

$U_{\psi^{g}}(O)=\frac{q^{g}(A\mid\textbf{L})}{f(A\mid\textbf{L})}\left\\{Y-m(A,\textbf{L})\right\\}+m(A,\textbf{L})(1-L^{\ast})+m(1,\textbf{L})L^{\ast}-\psi^{g}$

which is a useful representation for deriving doubly robust estimators. By
Corollary 1.1, estimators based on the EIF for $\psi^{g}$ in this case will be
doubly robust because $\nu_{1}$ can be estimated non-parametrically, and EIF-
based estimators for $\nu_{2}$ are doubly robust. In particular, these will be
consistent if either $f(A\mid\operatorname*{\textbf{L}})$ or
$m(A,\operatorname*{\textbf{L}})$ is consistently estimated, not necessarily
both.

######  Example 2.

Representative interventions for $J=1$. The intervention treatment
distribution is given by
$q^{\text{g}}(a\mid{\operatorname*{\textit{{l}}}})=f(a\mid{\operatorname*{\textit{{l}}}},R=1)$
where $R=I(A\geq\delta)$.

In this case we have

$\displaystyle\psi^{g}$
$\displaystyle=E_{{\operatorname*{\textbf{L}}}}\left\\{\sum\nolimits_{a=0}^{1}E(Y\mid
a,{\operatorname*{\textbf{L}}})q^{g}(a\mid{\operatorname*{\textbf{L}}})\right\\}=E_{{\operatorname*{\textbf{L}}}}\left[E_{A}\left\\{E(Y\mid
A,{\operatorname*{\textbf{L}}},R=1)\mid{\operatorname*{\textbf{L}}},R=1\right\\}\right]=E_{{\operatorname*{\textbf{L}}}}\left\\{E(Y\mid{{\operatorname*{\textbf{L}}}},R=1)\mid\right\\}$

Following previous results [17], we can see that, for this choice of
intervention treatment distribution and $J=1$, $\psi^{g}$ is only a function
of $R$, a coarsening of $A$, and takes the same form as the g-formula indexed
by the static deterministic strategy that sets treatment to 1 but with $R$
playing the role of treatment. We can apply Theorem 6 and Corollary 1.1
replacing the treatment $A$ with $R$ and $a^{\ast}$ with $r^{\ast}$.
Specifically, selecting $r^{\ast}=1$, $c_{2}=1$, $h_{1}(O)=0$, $h_{2}(O)=Y$,
by Theorem 6 we have

$\displaystyle\psi^{g}$ $\displaystyle=c_{2}{E[E\\{h_{2}(O)\mid
R=r^{\ast},L\\}]}=\underbrace{E\left\\{E(Y\mid
R=1,\operatorname*{\textbf{L}})\right\\}}_{\nu_{2}}$

and the EIF for $\psi^{g}$ is

$U_{\psi^{g}}(O)=\frac{I(R=1)}{f(R\mid\operatorname*{\textbf{L}})}\left\\{Y-m(R,\operatorname*{\textbf{L}})\right\\}+m(1,\operatorname*{\textbf{L}})-\psi^{g}$

By Corollary 1.1, estimators based on the EIF will be model doubly robust,
i.e., consistent if either models for $f(R\mid L)$ or $E(Y\mid R,L)$ are
correctly specified.

######  Example 3.

Multiplicative shift incremental propensity score interventions for $J=1$. The
intervention treatment distribution is given by
$q^{g}(a\mid\operatorname*{\textit{{l}}})=(1-\delta)al^{\ast}+f(a\mid\operatorname*{\textit{{l}}})(l^{\ast}\delta+1-l^{\ast})$.
Note the intervention distribution in Example 1 is a special case with
$\delta=0$.

In this case, for a choice of $\delta\in\\{0,1\\}$ we have

$\displaystyle\psi^{g}(\delta)$
$\displaystyle=E_{\operatorname*{\textbf{L}}}\left\\{\sum\nolimits_{a=0}^{1}E(Y\mid
a,\operatorname*{\textbf{L}})q^{g}(a\mid\operatorname*{\textbf{L}})\right\\}$
$\displaystyle=E_{\operatorname*{\textbf{L}}}\left[\sum\nolimits_{a=0}^{1}E(Y\mid
a,\operatorname*{\textbf{L}})\left\\{f(a\mid\operatorname*{\textbf{L}})(L^{\ast}\delta+1-L^{\ast})+L^{\ast}a(1-\delta)\right\\}\right]$
$\displaystyle=E_{\operatorname*{\textbf{L}},A}[E\left\\{Y(L^{\ast}\delta+1-L^{\ast})\right\\}\mid
A,\operatorname*{\textbf{L}}]+E_{\operatorname*{\textbf{L}}}\left[E\\{YL^{\ast}(1-\delta)\mid
A=1,L\\}\right]$

Selecting $a^{\ast}=1$, $c_{1}=1$, $c_{2}=(1-\delta)$,
$h_{1}(O)=Y(L^{\ast}\delta+1-L^{\ast})$, $h_{2}(O)=YL^{\ast}$, we have

$\displaystyle\psi^{g}(\delta)$
$\displaystyle=c_{1}{E\\{h_{1}(O)\\}}+c_{2}{E[E\\{h_{2}(O)\mid A=a,L\\}]}$
$\displaystyle=\underbrace{E\\{Y(L^{\ast}\delta+1-L^{\ast})\\}}_{\nu_{1}}+(1-\delta)\underbrace{E[E\\{YL^{\ast}\mid
A=1,L\\}]}_{\nu_{2}}$

by Theorem 6 and the EIF for $\psi^{g}(\delta)$ is given by

$\displaystyle
U_{\psi^{g}(\delta)}(O)=Y(L^{\ast}\delta+1-L^{\ast})+(1-\delta)\left[\frac{L^{\ast}A}{f(A\mid\operatorname*{\textbf{L}})}\left\\{Y-m(A,\operatorname*{\textbf{L}})\right\\}+m(1,\operatorname*{\textbf{L}})L^{\ast}\right]-\psi^{g}(\delta).$

This can be re-expressed as

$\displaystyle
U_{\psi^{g}(\delta)}(O)=\frac{q^{g}(A\mid\textbf{L})}{f(A\mid\textbf{L})}\left\\{Y-m(A,\textbf{L})\right\\}+m(A,\textbf{L})(L^{\ast}\delta+1-L^{\ast})+m(1,\textbf{L})L^{\ast}(1-\delta)-\psi^{g}(\delta)$

which is useful for deriving doubly robust estimators. By Corollary 1.1, the
estimators based on the EIF for $\psi^{g}(\delta)$ will be model doubly
robust. As Kennedy (2019) [16] noted, the EIF for $\psi^{g}$ indexed by an
odds shift (1) is not model doubly robust and, therefore, does not meet the
conditions of Corollary 1.1.

Note that the conditions of Corollary 1.1 are sufficient for model double
robustness of the EIF but are not necessary conditions. In Web Appendix B, we
consider model double robustness and the EIF for deterministic strategies that
depend on the natural value of treatment discussed in Section 4 [14, 15, 18,
19]. These are examples of $\psi^{g}$ indexed by intervention treatment
distributions that do not meet the conditions of Corollary 1.1 yet model
doubly robust estimators still exist.

###  6.2. Time-varying treatments

Recently, Molina (2017)[28] showed that, in time-varying treatment settings,
estimators derived from the EIF for a $\psi^{g}$ indexed by any intervention
treatment distribution that does not depend on the observed treatment process
[29, 6, 7] confer more protection against model misspecification than model
double robustness. Rather, they showed that these estimators are $J+1$ model
multiply robust, which implies model double robustness. The following Theorem
gives a sufficient condition for the existence of $J+1$ model multiple robust
estimators of $\psi^{g}$ when the intervention treatment distribution may
depend on the observed treatment process, and a simple approach to deriving
the EIFs for a particular class of such intervention treatment distributions.

######  Theorem 2.

Suppose an intervention treatment distribution can be written as the
following:

$\displaystyle q^{g}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=$ $\displaystyle
c_{1}h_{1}(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})I(a_{j}=a_{j}^{\ast})+c_{2}h_{2}(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})+$ (8) $\displaystyle
c_{3}h_{3}(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})p^{*}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})$

where $a_{j}^{\ast}$, $c_{1},\leavevmode\nobreak\ c_{2}$ and $c_{3}$ are
constants; $h_{1}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j-1})$,
$h_{2}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j-1})$ and
$h_{3}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j-1})$ are known
measurable functions of
$(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j-1})$; and $p^{*}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})$ is a non-degenerate
known probability distribution for $A_{j}$. Then the EIF for $\psi^{g}$
indexed by this intervention treatment distribution is

$\displaystyle
U_{\psi^{g}}(O)=\sum_{j=1}^{J}({T_{j}}-Q_{j-1})\prod_{k=0}^{j-1}\frac{q^{g}(A_{k}\mid
Y_{k}=1,\bar{\operatorname*{\textbf{L}}}_{k},\bar{A}_{k-1})}{f(A_{k}\mid
Y_{k}=1,\bar{\operatorname*{\textbf{L}}}_{k},\bar{A}_{k-1})}+T_{0}-\psi^{g}$
(9)

where $Q_{j}\equiv
Q_{j}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},\bar{Y}_{j})$ and
$T_{j}\equiv
T_{j}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},\bar{Y}_{j})$ are
iteratively defined from $j=J-1,\ldots,0$ such that for $T_{J}\equiv Y_{J}$,
we have $Q_{j}\equiv
E(T_{j+1}\mid\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},\bar{Y}_{j})$
and

$\displaystyle T_{j}=$ $\displaystyle
c_{1}Q_{j}^{A_{j}=a_{j}^{\ast}}h_{1}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j-1})+c_{2}Q_{j}h_{2}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j-1})+c_{3}\left\\{\sum_{a_{j}}p^{*}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{A}_{j-1})Q_{j}^{A_{j}=a_{j}}\right\\}h_{3}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j-1})$

with $Q_{j}^{A_{j}=a^{\ast}_{j}}\equiv
Q_{j}(\bar{\operatorname*{\textbf{L}}}_{j},A_{j}=a_{j}^{\ast},\bar{A}_{j-1},\bar{Y}_{j})$.
Estimators based on this EIF are $J+1$ model multiply robust in that they are
consistent if models for $Q_{j}$ are correctly specified for
$j={k,\ldots,J-1}$ and the observed treatment models are correctly specified
from $j=0,\ldots,k-1$ (for $k=0,\ldots,J$), where $j=s,s-1$ is $\emptyset$
$\forall s$.

Theorem 2 makes the derivation of the EIF and the corresponding estimators far
more straightforward and accessible when intervention distributions are in the
form given by (8). In Web Appendix D we prove that Expression (9) is the EIF
under a nonparametric model that imposes no restriction on the observed data
law for $\psi^{g}$ indexed by (8). In Web Appendix E we prove that estimators
based on this EIF are $J+1$ model multiply robust. Note that, by the
monotonicity of the survival indicators, we have $Y_{j+1}=Y_{j+1}Y_{j}$. This
implies that
$Q_{j}=Y_{j}Q_{j}=Y_{j}Q_{j}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},{Y}_{j}=1)$,
where
$Q_{j}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},{Y}_{j}=1)=E(T_{j+1}\mid\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},{Y}_{j}=1)$.
We now illustrate Theorem 2 in an example where the intervention treatment
distribution depends on the observed treatment process.

######  Example 4.

Consider the multiplicative shift incremental propensity score interventions
from Section 4, recalling the intervention distribution is $q^{g}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=(1-\delta)l_{j}^{\ast}a_{j}+(l_{j}^{\ast}\delta+1-l_{j}^{\ast})f(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})$.

This intervention distribution can be written in the form of Equation (8) by
selecting $a_{j}^{\ast}=1$, $c_{1}=1-\delta,\leavevmode\nobreak\
c_{2}=1,\leavevmode\nobreak\
h_{1}(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=l_{j}^{\ast},\leavevmode\nobreak\
h_{2}(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=l_{j}^{\ast}\delta+1-l_{j}^{\ast},\leavevmode\nobreak\
h_{3}(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=0$. By Theorem 2,
the EIF for this intervention distribution is then given by:

$\displaystyle U_{\psi^{g}(\delta)}(O)=$
$\displaystyle(Y_{J}-Q_{J-1})\prod_{j=0}^{J-1}\frac{q^{g}(A_{j}\mid
Y_{j}=1,\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j-1})}{f(A_{j}\mid
Y_{j}=1,\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j-1})}+$ (10)
$\displaystyle\sum_{j=1}^{J-1}\Big{\\{}\underbrace{(1-\delta)Q_{j}^{A_{j}=1}L_{j}^{\ast}+Q_{j}(L_{j}^{\ast}\delta+1-L_{j}^{\ast})}_{T_{j}}-Q_{j-1}\Big{\\}}\prod_{k=0}^{j-1}\frac{q^{g}(A_{k}\mid
Y_{k}=1,\bar{\operatorname*{\textbf{L}}}_{k},\bar{A}_{k-1})}{f(A_{k}\mid
Y_{k}=1,\bar{\operatorname*{\textbf{L}}}_{k},\bar{A}_{k-1})}+$
$\displaystyle\underbrace{(1-\delta)Q_{0}^{A_{0}=1}L_{0}^{\ast}+Q_{0}(L_{0}^{\ast}\delta+1-L_{0}^{\ast})}_{T_{0}}-\psi^{g}$

It is also straightforward to see that any $g$ corresponding to a
deterministic static treatment rule meets the conditions of Theorem 2 by
selecting
$h_{2}(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=h_{3}(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=0$,
$h_{1}(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=1$ and $c_{1}=1$.
In Web Appendix E, we further illustrate the application of Theorem 2 to
deterministic dynamic treatment rules, as well as other examples of
intervention distributions that depend on the observed treatment process for
the time-varying case including representative interventions and dynamic
treatment initiation strategies with a grace period. Note that, in these
examples and the incremental propensity score intervention example above,
Theorem 2 holds by selecting
$h_{3}(\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})=0$ such that
$p^{*}(a_{j}\mid 1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})$ need
not be specified. More generally, the applicability of Theorem 2 may require
specification of $p^{*}(a_{j}\mid
1,\bar{\operatorname*{\textit{{l}}}}_{j},\bar{a}_{j-1})$. For example, this
applies to an alternative grace period strategy where initiation within the
grace period is assigned such that there is a uniform probability of
initiating at each $j$ [4].

Note that the EIF given in Diaz et al. (2020) [19] cover a comprehensive class
of interventions that also guarantees estimators with model double robustness,
including interventions that depend on the natural value of treatment [14, 15]
also known as modified treatment policies. Following the results of Theorem 2
in Diaz et al. (2020), the EIF of the implied modified treatment policy from
our proposed intervention necessarily involve randomizer terms, but their
derivation of the corresponding EIF assumes that the distributions of the
randomizers are not known, when they certainly will be. Our Theorem 2 provides
the EIF for the g-formula indexed by a class of stochastic interventions that
may depend on the observed treatment process. It can be shown that nearly all
functionals in this class are captured by the g-formula functionals for which
Diaz et al. (2020) provides the EIF. Diaz et al.’s results would capture all
of the functionals in this class, including those indexed by our proposed
multiplicative incremental propensity score interventions, provided they
projected their EIF onto a tangent space corresponding to smaller models
whereby the distributions of some conceptualized randomizers are known. We do
not take this approach and allow one to derive the EIF directly from a
stochastic treatment distribution without requiring one to define an implied
modified treatment policy first111There are multiple ways to define a modified
treatment policy that has an identifying formula to equal the g-formula for
our multiplicative shift intervention. It will involve randomizer terms that
can potentially also depend on $L_{j}^{*}$. This alternative derivation of the
EIF may be more intuitive for treatment distributions that do not depend on
the natural value of treatment.

Finally, in Web Appendix C we use a similar line of reasoning to Theorem 6 and
Corollary 1.1 to derive the EIF for $\psi^{g}$ and to assess the existence of
doubly robust estimators for $\psi^{g}$ indexed by an intervention
distribution that depends on the observed treatment process when $J=2$ with
examples. However, this approach to deriving the EIF is cumbersome for large
$J$, providing no simplification over Theorem 2.

##  7\. Estimators of $\psi^{g}$ indexed by multiplicative shift incremental
propensity score interventions

In this section, we consider various estimators of $\psi^{g}$ under the
multiplicative shift incremental propensity score interventions defined by
(2).

###  7.1. EIF-based estimators

Several EIF-based estimators for $\psi^{g}$ have been proposed for
deterministic treatment interventions including Bang and Robins (2005)’s
estimator [6, 30, 31], weighted ICE estimator [32, 12] and targeted maximum
likelihood estimator (TMLE) [7, 33, 34]. Weighted ICE and TMLE are variations
of Bang and Robins (2005)[6]. Compared with Bang and Robins (2005)[6],
weighted ICE can give better performance [2]. Unlike both Bang and Robins
(2005)[6] and weighted ICE, TMLE can incorporate machine learning algorithms
[2]. In the absence of machine learning algorithms, weighted ICE and TMLE
perform similarly [2], but weighted ICE is easier to implement. In this
section we will consider two estimators: 1) weighted ICE estimator that uses
parametric models to estimate the nuisance functions thereby allowing for
$J+1$ model multiple robustness; and 2) TMLE that also uses sample-splitting
and cross fitting [26, 35, 10] to allow one to incorporate machine learning
algorithms to estimate the the nuisance functions.

####  7.1.1. Weighted ICE estimator

Let $\pi_{j}\equiv f(A_{j}\mid
Y_{j}=1,\bar{A}_{j-1},\bar{\operatorname*{\textbf{L}}}_{j})$ and let
$\pi_{j}(\alpha_{j})=f(A_{j}\mid
Y_{j}=1,\bar{A}_{j-1},\bar{\operatorname*{\textbf{L}}}_{j};\alpha_{j})$ be a
working parametric model for $\pi_{j}$ with
$\alpha=(\alpha_{0},\ldots,\alpha_{J-1})$. Denote estimates
$\hat{\pi}_{j}\equiv\pi_{j}(\hat{\alpha}_{j})$ of $\pi_{j}$ with
$\hat{\alpha}_{j}$ the maximum likelihood estimate (MLE) of ${\alpha}_{j}$
computed from the observed data. Subsequently, let $\hat{q}^{g}_{j}\equiv
q^{g}_{j}(\hat{\pi}_{j})$ be an estimate of $q^{g}(A_{j}\mid
Y_{j}=1,\bar{A}_{j-1},\bar{\operatorname*{\textbf{L}}}_{j})$ as defined in (2)
for a choice of $\delta\in[0,1]$, replacing the observed treatment process
with the estimate $\hat{\pi}_{j}$. Let $\hat{Q}_{j}$ be a working parametric
model for $Q_{j}$ defined in Theorem 2. In the following algorithm, each
$\hat{T}_{j}$ is calculated by replacing $Q_{j}$ in formula (10) with the
estimate $\hat{Q}_{j}$. The weighted ICE algorithm is specifically implemented
as follows:

1:Compute the MLEs $\hat{\alpha}$ of ${\alpha}$ from the observed data. Set
$\hat{T}_{J}=Y_{J}$.

2:Recursively from $j=J-1,\ldots,0$:

  1. (A)

Fit a logistic regression model
$Q_{j}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},Y_{j}=1;\theta_{j})=\operatorname*{expit}\\{\theta_{j}^{T}\phi(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j})\\}$
for
$E(\hat{T}_{j+1}\mid\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},Y_{j}=1)$
with observational weight $\prod_{k=0}^{j}({\hat{q}^{g}_{k}}/{\hat{\pi}_{k}})$
in those who survive by time $j$. Here,
$\phi(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j})$ is a known function
of $\bar{\operatorname*{\textbf{L}}}_{j}$ and $\bar{A}_{j}$. More
specifically, we solve for $\theta_{j}$ in the following estimating equation:

$\mathbb{P}_{n}\left[Y_{j}\prod_{k=0}^{j}\frac{\hat{q}^{g}}{\hat{\pi}_{k}}\phi_{j}({\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j}})\left\\{\hat{T}_{j+1}-Q_{j}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},Y_{j}=1;{\theta}_{j})\right\\}\right]=0$

  2. (B)

Compute $\hat{T}_{j}$ from $\hat{Q}_{j}\equiv
Q_{j}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},\bar{Y}_{j};\hat{\theta}_{j})$
ensuring $\hat{T}_{j}=0$ when $Y_{j}=0$.

3:Estimate $\hat{\psi}^{g}(\delta)_{WICE}=\mathbb{P}_{n}(\hat{T}_{0})$

Algorithm 1 Algorithm for Weighted ICE

where $\mathbb{P}_{n}\\{f(X)\\}=n^{-1}\sum_{i=1}^{n}f(X_{i})$. Following
arguments in Section 6.2, this estimator is $J+1$ model multiply robust.

####  7.1.2. TMLE with sample-splitting and cross-fitting

This algorithm utilizes sample-splitting and cross-fitting to allow flexible
machine learning algorithms for estimating nuisance functions while
circumventing Donsker class conditions [36, 26]. In Web Appendix G, we prove
the asymptotic normality of this estimator under the condition that the
nuisance functions are estimated consistently at rates faster than $n^{-1/4}$
when $\psi^{g}$ is indexed by the interventions (2).

Suppose that a sample of size $n$ is split into $M$ disjoint subsets. Let
$S_{m}$ denote the subset of individuals in split $m=1,\ldots,M$ and let
$S_{-m}$ denote individuals not in split $m$ (i.e., $S_{-m}=\\{i\notin
S_{m}\\}$). Moreover, let $\hat{\pi}_{j}^{(-m)}$, $\hat{q}_{j}^{(-m)}$ and
$\hat{Q}_{j}^{(-m)}$ denote estimates of $\pi_{j},\leavevmode\nobreak\
q_{j}^{g}$ and $Q_{j}$ obtained from machine learning algorithms to
individuals in $S_{(-m)}$.

1:For each $m=1,\ldots,M$:

  1. (A)

For individuals in $S_{-m}$: compute $\hat{\pi}_{j}^{(-m)}$, $\forall j$. Set
$\hat{T}_{J}=Y_{J}$.

  2. (B)

Recursively from $j=J-1,\ldots,0$ for individuals in $S_{-m}$:

    1. (a)

Compute
$\hat{Q}_{j}^{(-m)}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},Y_{j}=1)$
by regressing $\hat{T}_{j+1}$ on
$(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j})$ in those alive at time
$j$

    2. (b)

Compute $\hat{T}_{j}$ from
$\hat{Q}_{j}^{(-m)}\equiv\hat{Q}_{j}^{(-m)}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},\bar{Y}_{j})$
by formula (10), setting $\hat{T}_{j}=0$ if ${Y}_{j}=0$

  3. (C)

For individuals in $S_{m}$, set $\hat{T}_{J}=Y_{J}$. Then recursively from
$j=J-1,\ldots,0$:

    1. (a)

Solve for $\gamma_{j}$ in the following set of estimating equations:

$\mathbb{P}_{n}^{m}\left(Y_{j}\prod_{k=0}^{j}\frac{\hat{q}_{k}^{{g}^{(-m)}}}{\hat{\pi}_{k}^{(-m)}}\left[\hat{T}_{j+1}-\operatorname*{expit}\left\\{\operatorname*{logit}\left(\hat{Q}_{j}^{(-m)}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},Y_{j}=1)\right)+\gamma_{j}\right\\}\right]\right)=0$

    2. (b)

Compute $\hat{T}_{j}$ from
$\hat{Q}_{j}^{\Delta}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},Y_{j}=1)\equiv\operatorname*{expit}\left\\{\operatorname*{logit}\left(\hat{Q}_{j}^{(-m)}(\bar{\operatorname*{\textbf{L}}}_{j},\bar{A}_{j},Y_{j}=1)\right)+\hat{\gamma}_{j}\right\\}$
if $Y_{j}=1$, otherwise set $\hat{T}_{j}=0$ if ${Y}_{j}=0$

2:Calculate
$\hat{\psi}^{g}(\delta)_{TMLE}=\frac{1}{M}\sum_{m=1}^{M}\mathbb{P}_{n}^{m}(\hat{T}_{0})$

Algorithm 2 Algorithm for TMLE with sample-splitting and cross-fitting

Here $\mathbb{P}_{n}^{m}\\{f(X)\\}=\frac{1}{\mid S_{m}\mid}\sum_{i\in
S_{m}}f(X_{i})$ where ${\mid S_{m}\mid}=n/M$ is the cardinality of $S_{m}$.

###  7.2. Singly robust estimators

We also consider less optimal but computationally simple singly robust
estimators of $\psi^{g}$ indexed by (2). An IPW estimator
$\hat{\psi^{g}}_{IPW}(\delta)$ can be obtained by the product
$\hat{\psi}^{g}_{IPW}(\delta)=\prod_{j=0}^{J-1}\hat{\Upsilon}^{g}_{IPW,j}(\delta)$
where $\hat{\Upsilon}^{g}_{IPW,j}(\delta)$ can be interpreted as an estimate
of the discrete hazard at $j$ under a stochastic strategy $g$ where treatment
assignment is a draw from (2) given the identifying conditions of Section 5.
Each $\hat{\Upsilon}^{g}_{IPW,j}(\delta)$ can be obtained by solving for
$\Upsilon^{g}_{IPW,j}(\delta)$ in the following estimating equations:

$\mathbb{P}_{n}\left[Y_{j}\prod_{k=0}^{j}\frac{\hat{q}^{g}_{k}}{\hat{\pi}_{k}}\\{Y_{j+1}-\Upsilon^{g}_{IPW,j}(\delta)\\}\right]=0$

Alternatively, the singly-robust ICE estimator, which we will denote
$\hat{\psi}^{g}_{ICE}(\delta)$, can be obtained as a special case of the
algorithm for weighted ICE above where the observational weights are set to 1.

###  7.3. Censoring

Straightforward extensions of the identification arguments in Section 5 in
studies with censoring follow by implicitly including in $g$ a hypothetical
intervention that eliminates censoring throughout follow-up [22] with
straightforward extensions of the g-formula $\psi^{g}$, properties of its EIF
and associated estimation procedures. Briefly, denote $C_{j}$ as the indicator
of censoring by time $j$ and adopt the order $(L_{j},A_{j},C_{j+1},Y_{j+1})$.
Extensions to accommodate censoring for singly robust weighted estimators and
the various EIF-based estimators considered, require, in addition to
estimating $\alpha_{j}$ in $f(A_{j}\mid
Y_{j}=1,\bar{A}_{j-1},\bar{\operatorname*{\textbf{L}}}_{j},C_{j}=0;\alpha_{j})$,
also estimating $\alpha^{c}_{j}$ in $P(C_{j+1}=1\mid
C_{j}=0,\bar{A}_{j},\bar{\operatorname*{\textbf{L}}}_{j},Y_{j}=1;\alpha^{c}_{j})$
for $j=0,\ldots,J-1$ with $\alpha^{c}=(\alpha^{c}_{1},\ldots,\alpha^{c}_{J})$.
Further details of modifications to the weighted ICE and TMLE to accommodate
censoring are provided in Web Appendix F.

##  8\. Simulation studies

We conducted two different simulation studies. The first simulation study aims
to compare the performance of the weighted ICE, IPW and ICE estimators when
the nuisance functions are estimated through parametric models under various
model misspecification scenarios. The second simulation study aims to compare
the performance of TMLE with sample splitting and cross fitting, IPW and ICE
when the nuisance functions are estimated through machine learning algorithms.

###  8.1. Simulation study 1: using parametric models

In this simulation study we compare the performance of the weighted ICE
estimator with the singly robust estimators (IPW and ICE estimators) for
$\psi^{g}$ indexed by the intervention distribution (2) which, under
identifying conditions discussed in Section 5, equals the cumulative
probability of survival at $J$ under an intervention that increases the
probability of treatment initiation in those with $L_{j}^{\ast}=1$ as a
function $\delta$. Recall that this increase is defined such that decreasing
values of $\delta$ correspond to an increasing probability of treatment
initiation (with $\delta=1$ coinciding with no treatment intervention).

We simulated 1000 samples of $n=(500,\leavevmode\nobreak\
1000,\leavevmode\nobreak\ 2500)$ individuals selecting $J=5$ and
$\delta=(0.75,\leavevmode\nobreak\ 0.50,\leavevmode\nobreak\ 0.25)$. We
simulated the following variables:
$(\operatorname*{\textbf{L}}_{0},A_{0},C_{1},Y_{1},\operatorname*{\textbf{L}}_{1},A_{1},...,C_{5},Y_{5})$,
where $\operatorname*{\textbf{L}}_{j}=(L_{j}^{\ast},L_{1j},L_{2j})$ is the
vector of measured confounders. Specifically, we generated $L_{1}^{\ast}$ and
$L_{11}\sim\text{Ber}\\{\text{expit}(-1)\\}$, and
$L_{2}\sim\text{Ber}\\{\text{expit}(1+L_{1}^{\ast})\\}$. The censoring
indicator at each time $j$ ($j=1,\ldots,5$) was simulated from
$C_{j}\sim\text{Ber}\\{\text{expit}(-2+L_{1j}-L_{2j})\\}$ if $C_{j-1}=0$ and
$Y_{j}=1$. The outcome at each time $j$ ($j=1,\ldots,5$) is simulated from
$Y_{j}\sim\text{Ber}\\{\text{expit}(1+3A_{j-1}-2L_{j-1}^{\ast}+L_{1,j-1}-L_{2,j-1})\\}$
if $Y_{j-1}=1$ and $C_{j}=0$. The time-varying confounders at time $j$
($j=0,\ldots,4$) are simulated from
$L_{j}^{\ast}\sim\text{Ber}\\{\text{expit}(-1-A_{j-1}+L_{j-1}-L_{1,j-1}+L_{2,j-1})\\}$,
$L_{1j}\sim\text{Ber}\\{\text{expit}(-1+A_{j-1}+L_{1,j-1}-L_{2,j-1})\\}$ and
$L_{2j}\sim\text{Ber}\\{\text{expit}(1+A_{j-1}+L_{j}^{\ast}+L_{2,j-1})\\}$ if
$Y_{j}=1$. Treatment at time $j$ ($j=0,\ldots,4$) is simulated from
$A_{j}\sim\text{Ber}\\{\text{expit}(-1-2L_{j}^{\ast}-L_{1j}+L_{2j}+2A_{j-1})\\}$
if $Y_{j}=1$. In addition
$(Y_{j},\operatorname*{\textbf{L}}_{j},A_{j},\ldots)=(\emptyset,\emptyset,\emptyset,\ldots)$
if $C_{j}=1$, and
$(\operatorname*{\textbf{L}}_{j},A_{j},C_{j+1},\ldots)=(\emptyset,\emptyset,\emptyset,\ldots)$
if $Y_{j}=0$.

The true cumulative probabilities of survival were calculated by using the
true parametric models to generate a Monte Carlo sample of size $10^{7}$ under
all interventions of interest. Our selection of parameters resulted in a
scenario where selecting smaller $\delta$ (that is, interventions with larger
increases in the probability of treatment initiation at each $j$) improves
survival.

We considered three estimation scenarios for each choice of $\delta$ and
sample size such that (1) all models are correctly specified, (2) only the
outcome regression models are correctly specified, and (3) only the treatment
(propensity score) and censoring models are correctly specified. The true
functional forms of the treatment and censoring models are known under our
simulation because treatment and censoring were generated to only depend on
past measured variables. Similarly, the functional form of the outcome
regression model for $Q_{J-1}=E(Y_{J}\mid
Y_{J-1}=1,C_{J}=0,\bar{A}_{J-1},\overline{L}_{J-1})$ is known due to the
absence of unmeasured common causes. However, the true functional forms of the
outcome regressions $Q_{j}$ for $0\leq j<J-1$ are not known under our
simulation. To ensure correctly specified models for $Q_{j}$, $0\leq j<J-1$,
saturated models were fit, i.e., all main terms and interaction terms for
$(A_{j},L_{j}^{\ast},L_{1j},L_{2j})$. In scenarios with misspecified models,
at each time $j$, the misspecified treatment model ignores the censoring
process and excludes $A_{j-1}$ in the model, and the misspecified outcome
regression model excludes any pairwise interactions between the covariates and
treatment.

Figure 1. Results for the Simulation study 1

Figure 1 compares performance of the three estimators of $\psi^{g}$ indexed by
(2) for $\delta=(0.75,\leavevmode\nobreak\ 0.50,\leavevmode\nobreak\ 0.25)$.
Complementary results are given in Tables 3–5 in Web Appendix H. As expected,
all estimators were nearly unbiased under correctly specified models. Under
our model misspecification scenarios, $\hat{\psi}^{g}(\delta)_{WICE}$ is
nearly unbiased, but the IPW estimator is biased when the treatment models are
misspecified, and the ICE estimator is biased when the outcome models are
misspecified. In addition, under correctly specified models
$\hat{\psi}^{g}(\delta)_{ICE}$ is the most efficient, and
$\hat{\psi}^{g}(\delta)_{IPW}$ is the least efficient estimator.
Interestingly, the simulation results show that
$\hat{\psi}^{g}(\delta)_{WICE}$ has smaller MSE than the IPW estimator in all
scenarios.

The simulation results also show that as $\delta$ decreases, the standard
error (and MSE) in all three estimators increases. This is due to an increase
in the effect of near positivity violations as $\delta$ nears zero. In fact,
we would expect all three estimators to have the largest standard errors when
$\delta=0$, which is equivalent to a strategy that treats all individuals with
$L_{j}^{\ast}=1$ at all times. In Web Appendix H, we also show $J+1$ model
robustness of our weighted ICE estimator in a model misspecification scenario
that requires more than model double robustness.

###  8.2. Simulation study 2: using machine learning methods

In the second simulation, we compare the performance of algorithms that use
machine learning to estimate the nuisance functions for $\psi^{g}$ indexed by
(2) with $J=5$. Specifically, we compare TMLE with sample splitting and cross
fitting, IPW and ICE. Given much longer computation times, we limited
consideration to one choice of $\delta=0.5$. Unlike Simulation 1, we add model
complexity to the data generating mechanism by considering continuous
covariates, which might mimic real-life data more closely. We simulated 1000
hypothetical cohorts of $n=(250,\leavevmode\nobreak\ 500,\leavevmode\nobreak\
1000)$ comprising the following variables:
$(\operatorname*{\textbf{L}}_{0},A_{0},C_{1},Y_{1},\operatorname*{\textbf{L}}_{1},A_{1},...,C_{5},Y_{5})$,
where $\operatorname*{\textbf{L}}_{j}=(L_{1j},L_{j}^{\ast})$. In addition,
$\operatorname*{\textbf{L}}_{0}=(L_{0}^{1},L_{0}^{2},L_{10},L_{0}^{\ast})$,
where $L_{0}^{1}$ and $L_{0}^{2}$ are baseline covariates. In particular,
$L_{0}^{1}\sim\text{Ber}(0.5)$, $L_{0}^{2}\sim\text{N}(0,1)$,
$L_{10}\sim\text{N}(2+L_{0}^{1},1)$ and
$L_{0}^{\ast}\sim\text{Ber}\\{\text{expit}(1.5-0.5L_{0}+L_{0}^{1}+0.25L_{0}^{2})\\}$.
For $j\geq 1$,
$L_{1j}\sim\text{N}(2+A_{j-1}-L_{j-1}^{\ast}+0.5L_{1,j-1}+L_{0}^{1},1)$ and
$L_{j}^{\ast}\sim\text{Ber}\\{\text{expit}(1.5-A_{j-1}-0.5L_{1j}+L_{j-1}^{\ast}+L_{0}^{1}+0.25L_{0}^{2})\\}$
if $Y_{j}=1$. Censoring indicator at each time $j$ ($j=1,\ldots,5$) is
simulated from
$C_{j}\sim\text{Ber}[\text{expit}\\{-4-A_{j-1}-L_{j-1}^{\ast}-0.5\sqrt{\mid
L_{1,j-1}L_{0}^{2}\mid}+1.5\mid L_{1,j-1}\mid/(1+\exp(L_{0}^{2}))\\}]$ if
$C_{j-1}=0$ and $Y_{j}=1$. The outcome at each time $j$ ($j=1,\ldots,5$) is
simulated from
$Y_{j}\sim\text{Ber}[\text{expit}\\{-1+2A_{j-1}-2L_{j-1}^{\ast}+0.25L_{j-1}^{\ast}L_{1,j-1}+0.5L_{0}^{1}+0.75\mid
L_{1,j-1}+L_{0}^{2}\mid^{1.5}\\}]$ if $Y_{j-1}=1$ and $C_{j}=0$. Treatment at
time $j$ ($j=0,\ldots,4$) is simulated from
$A_{j}\sim\text{Ber}\\{\text{expit}(-3+L_{j}^{\ast}-0.5L_{1j}+0.25L_{j}^{\ast}L_{j}+0.5L_{0}^{1}+0.25L_{0}^{2}+0.5\mid
L_{0}^{2}\mid)\\}$ if $Y_{j}=1$ and $A_{j-1}=0$, and is set to $1$ if
$Y_{j}=1$ and $A_{j-1}=1$.

Nuisance functions were estimated using the Super Learner ensemble, which uses
cross validation to select the best convex combination of predictions from a
pool of prediction algorithms [37]. The library of potential candidates used
here consisted of: generalized linear models and its variants (SL.glm,
SL.glm.interaction), Bayesian generalized linear models (SL.bayesglm),
generalized additive models with smoothing splines (SL.gam), multivariate
adaptive regression Splines (SL.earth), neural networks (SL.nnet) and random
forest (SL.ranger).

Table 1 compares the performance of the 3 estimators. The ICE and IPW
estimators show bias as they are not expected to converge at $\sqrt{n}$ rates
when machine learning is used for nuisance parameter estimation. TMLE, on the
other hand, show little to no bias in all instances. This agrees with theory
as TMLE allows the nuisance functions to converge at slower nonparametric
rates. Moreover, the estimated coverage probability of the confidence
intervals for TMLE based on the asymptotic variance (see Web Appendix G) is
very close to the nominal 95%: $(94.7,\leavevmode\nobreak\
96.2,\leavevmode\nobreak\ 95.2,\leavevmode\nobreak\ 94.4)$ for
$n=(250,\leavevmode\nobreak\ 500,\leavevmode\nobreak\
1000,\leavevmode\nobreak\ 2500)$, respectively.

Table 1. Simulation study 2 for proposed treatment intervention distribution
and incorporating machine learning algorithms ($M=2$). True probability of
survival at time 5 is $0.629$. All values are multiplied by 100.

| $n=250$ | $n=500$ | $n=1000$ | $n=2500$
---|---|---|---|---
Estimator | BIAS | SE | RMSE | BIAS | SE | RMSE | BIAS | SE | RMSE | BIAS | SE | RMSE
$\hat{\psi}^{g}(\delta)_{ICE}$ | -1.50 | 4.35 | 4.61 | -0.82 | 2.83 | 2.95 | -0.47 | 2.14 | 2.19 | -0.16 | 1.38 | 1.39
$\hat{\psi}^{g}(\delta)_{IPW}$ | -1.50 | 4.91 | 5.13 | -1.50 | 3.32 | 3.64 | -1.35 | 2.60 | 2.93 | -1.10 | 1.71 | 2.03
$\hat{\psi}^{g}(\delta)_{TMLE}$ | -0.19 | 5.79 | 5.79 | -0.09 | 3.61 | 3.62 | -0.07 | 2.59 | 2.59 | 0.03 | 1.65 | 1.65

##  9\. Application

Randomized trials suggest that antiretroviral pre-exposure prophylaxis (PrEP)
is highly effective in preventing HIV infection among men who have sex with
men (MSM) [38, 39, 40]. At the same time, there is widespread concern that
PrEP may decrease condom use, thereby increasing incidence of bacterial
sexually transmitted infections (STIs) among MSM [41, 42]. With this backdrop,
PrEP uptake has been low in practice, particularly in MSM with markers for
higher HIV risk [43]. This is precisely the setting where near or true
positivity violations will occur if the analyst attempts to query
observational data about the effects of deterministic interventions such as
“always treat” versus “never treat” with PrEP because propensity scores will
be close to (or equal to) zero for individuals with certain levels of the
measured confounders. In conjunction with these challenges, such deterministic
treatment effects are not of greatest interest for treatments like PrEP, where
biological benefits are established but population disease burden may be
impacted with even small increases in treatment uptake.

We illustrate the application of the estimators of survival by $J$ discussed
in Section 7 using electronic health record data from the Cambridge Health
Alliance – a large community healthcare system in Eastern Massachusetts – to
estimate effects of increasing PrEP uptake on bacterial STI diagnosis.
Specifically, we consider interventions that, beginning at the time of an HIV
negative test, successfully increase the proportion initiating PrEP ($A_{j}$)
in each follow-up week $j$ only among those receiving an STI test and no prior
diagnosis of HIV at time $j$ ($L_{j}^{\ast}=1$, being tested for STIs suggests
recent condomless sex and PrEP would not be used after an HIV diagnosis). No
intervention is made for the remainder of the population at time $j$
($L_{j}^{\ast}=0$). Increases in treatment uptake under these interventions
are quantified by a specified $\delta\in[0,1]$ as defined in (7), which
quantifies the factor by which the probability of treatment non-initiation is
decreased (relative to no intervention) at $j$. We consider $J=26$ weeks and,
as in simulation study 1, consider $\delta=(0.95,\leavevmode\nobreak\
0.85,\leavevmode\nobreak\ 0.75)$ representing realistic interventions that
result in “low”, “medium” and “high” success in PrEP uptake relative to no
intervention. We use $\delta=1$ (corresponding to no intervention) as the
reference in defining causal effects.

Our analytic data set was restricted to patients who met all the following
inclusion criteria at some point during 2012–2017: 1) Cis male with report of
male gender of sex partner(s); 2) 15 years of age or older; 3) an HIV-negative
test; 4) had no PrEP prescription in the 3 months prior to baseline; and 5)
had no STI diagnosis in the 12 months prior to baseline. Baseline (week $j=0$)
for an individual was defined as the first week that all of these inclusion
criteria are met. For simplicity, we excluded one individual who met these
criteria but died without a bacterial STI diagnosis during the 26 week follow-
up period. Our final analytic data set consisted of $n=1103$ individuals. As
expected, few initiated PrEP over the follow-up (cumulatively 5.1% over the 26
weeks). The cumulative proportion of those receiving an STI test while being
HIV-free over the 26 weeks was 70.7%. Note that no individual was treated as
censored in this analysis, requiring additional assumptions that medical care
was not sought outside of the Cambridge Health Alliance by any individual
included at baseline over the 26 week follow-up.

Baseline covariates $\operatorname*{\textbf{L}}_{0}$ included age and calendar
year at baseline, race/ethnicity, and time-varying covariates
$\operatorname*{\textbf{L}}_{j}$ included indicator of any ambulatory
encounter, indicator of HIV, indicator of any HIV testing and indicator of any
STI testing. We used the Super Learner ensemble (with the same potential
candidates as in the simulation) to estimate all nuisance functions for TMLE
with $M=5$. We compared these results with the IPW, ICE and weighted ICE
estimators described in this paper where the nuisance functions are specified
by parametric models. Confidence intervals for each of the methods are
obtained from 1000 bootstrap samples by taking the 2.5th and 97.5th
percentiles of the resulting estimates.

Our estimate of the probability of not receiving an STI diagnosis under no
intervention by 26 week follow-up ($\delta=1$) was 93.7%. Table 2 shows
results from the four methods for $\delta<1$. In this case point estimates
from all of the methods are similar. The results do not provide sufficient
evidence that increasing PrEP uptake increases risk of STI diagnosis. For
instance, compared with no intervention ($\delta=1$), the relative survival
estimates under low, medium and high increases in PrEP uptake were 0.99 (95%
CI = $(0.96,\leavevmode\nobreak\ 1.01)$), 0.97 (95% CI =
$(0.91,\leavevmode\nobreak\ 1.01)$) and 0.96 (95% CI =
$(0.87,\leavevmode\nobreak\ 1.02)$), respectively under
$\hat{\psi}^{g}(\delta)_{TMLE}$. The relative survival estimates using other
estimators were very similar (see Web Appendix I). We also note that due to an
increase in the presence of near positivity violations as $\delta$ nears zero,
observational weights calculated under smaller $\delta$ were more variable
than larger $\delta$ (see Web Appendix I). We would expect standard errors
from all of the estimators to be the largest for $\delta=0$.

Table 2. Point estimates and 95% confidence intervals from analysis of MSM
from the Cambridge Health Alliance on the effect of incremental PrEP
initiation on incident STI diagnosis. All values are multiplied by 100.

|  $\hat{\psi}^{g}(\delta)_{TMLE}$ (with ML) | $\hat{\psi}^{g}(\delta)_{WICE}$ | $\hat{\psi}^{g}(\delta)_{ICE}$ | $\hat{\psi}^{g}(\delta)_{IPW}$
---|---|---|---|---
$\uparrow$ in PrEP | Est. | 95% C.I. | Est. | 95% C.I. | Est. | 95% C.I. | Est. | 95% C.I.
Low | 92.9 | (89.3, 95.2) | 93.0 | (90.9, 94.9) | 92.9 | (91.1, 94.6) | 92.9 | (90.5, 94.9)
Medium | 91.4 | (84.8, 95.5) | 91.6 | (87.1, 95,1) | 91.4 | (88.5, 94.1) | 91.3 | (85.8, 95.0)
High | 90.9 | (80.9, 95.8) | 90.3 | (83.2, 95.4) | 90.0 | (85.8, 93.8) | 89.9 | (80.9, 95.3)

##  10\. Discussion

Many methods that have been proposed for estimating causal estimands in time-
varying treatment settings for survival analysis, and among these methods are
estimators that offer protection against model misspecification and can also
attain semiparametric efficiency bound. However, most of these doubly robust
estimators have been in the setting of deterministic treatment interventions.
In this paper we provided some sufficient conditions for the existence of
doubly robust estimators when a treatment intervention distribution can depend
on the observed treatment process for point treatment processes. We also
discussed a class of intervention distributions that are always guaranteed to
give doubly/multiply robust estimators and gave a general form of the EIFs
that are associated with these intervention distributions. Among these
intervention distributions is our multiplicative shift incremental propensity
score intervention distribution, which aims to increase treatment uptake in a
group of individuals who are at high risk of the outcome but have low exposure
to treatment. We provided various estimators that can be used for our proposed
treatment intervention for both parametric and machine learning algorithms.

We conducted two simulation studies for our proposed multiplicative shift
intervention distribution. Our first study show that the weighted ICE is more
robust to model misspecification than IPW and ICE when the nuisance functions
are estimated using parametric models. Our second study show that the TMLE
with sample-splitting and cross-fitting is consistent as long as the nuisance
functions are estimated consistently at fast enough rates using machine
learning methods, which may not necessarily be $n^{-1/2}$. We also illustrated
our an application of our estimators to a real world dataset in the PrEP
context.

Note that our proposed intervention treatment distribution (2) is guaranteed
under a stochastic intervention such that treatment initiation status at time
$j$ for an individual with $L^{\ast}_{j}=1$ is a random draw from (2). The
identifying conditions reviewed in Section 5 are sufficient to give our effect
estimates this interpretation. A more policy relevant interpretation might,
for example, be an intervention where individuals with $L^{\ast}_{j}=1$ are
always offered PrEP counseling. In these individuals, the intervention
distribution (2) quantifies the hypothesized “success” of such an intervention
where $g$ in this case really refers to a deterministic strategy relative to
the unmeasured treatment “offered PrEP counseling”. Additional assumptions are
needed to give our effect estimates this interpretation following similar
arguments to those given in Richardson and Robins (2013)[20] and Young et al.
(2014) [18].

Finally, we note that while machine learning algorithms are more robust to
model form misspecification, they are also computationally complex and may be
practically infeasible for very large datasets without powerful computing
systems. Moreover, issues related to data privacy make access to advanced
computational resources impossible in many cases. Therefore, estimators that
offer model double/multiple robustness are useful in practice as they offer
protection against model misspecification and can be easily computed using
standard off-the-shelf regression software in R.

## Acknowledgment

This research was funded by NIAID grant R21AI143386-01A1. The authors thank
Dr. Gerard Coste at Cambridge Health Alliance for assistance with collection
of clinical data.

## References

  * [1] Robins JM.  A new approach to causal inference in mortality studies with a sustained exposure period: application to control of the healthy worker survivor effect.  _Mathematical Modelling_ 1986; 7: 1393–1512.
  * [2] Tran L, Yiannoutsos C, Wools-Kaloustian K et al.  Double robust efficient estimators of longitudinal treatment effects: Comparative performance in simulations and a case study.  _The International Journal of Biostatistics_ 2019; .
  * [3] Wen L, Young JG, Robins JM et al.  Parametric g-formula implementations for causal survival analyses.  _Biometrics_ 2021; .
  * [4] Cain L, Robins JM, Lanoy E et al.  When to start treatment? a systematic approach to the comparison of dynamic regimes using observational data.  _The International Journal of Biostatistics_ 2010; 6.
  * [5] Neugebauer R, Schmittdiel J and Van Der Laan M.  Targeted learning in real-world comparative effectiveness research with time-varying interventions.  _Statistics in Medicine_ 2014; 33: 2480–2520.
  * [6] Bang H and Robins J.  Doubly robust estimation in missing data and causal inference models.  _Biometrics_ 2005; 61: 962–973.
  * [7] Van Der Laan MJ and Rose S.  _Targeted learning: causal inference for observational and experimental data_.  Springer Science & Business Media, 2011.
  * [8] Robins JM, Li L, Tchetgen E et al.  Quadratic semiparametric von mises calculus.  _Metrika_ 2009; 69: 227–247.
  * [9] Robins JM, Li L, Tchetgen E et al.  Asymptotic normality of quadratic estimators.  _Stochastic processes and their applications_ 2016; 126: 3733–3759.
  * [10] Chernozhukov V, Chetverikov D, Demirer M et al.  Double/debiased machine learning for treatment and structural parameters, 2018.
  * [11] Stitelman OM, De G V and Van Der Laan MJ.  A general implementation of tmle for longitudinal data applied to causal inference in survival analysis.  _The International Journal of Biostatistics_ 2012; 8.
  * [12] Rotnitzky A, Robins J and Babino L.  On the multiply robust estimation of the mean of the g-functional.  _arXiv preprint arXiv:170508582_ 2017; .
  * [13] Taubman SL, Robins JM, Mittleman MA et al.  Intervening on risk factors for coronary heart disease: an application of the parametric g-formula.  _International Journal of Epidemiology_ 2009; 38: 1599–1611.
  * [14] Muñoz ID and Van Der Laan M.  Population intervention causal effects based on stochastic interventions.  _Biometrics_ 2012; 68: 541–549.
  * [15] Haneuse S and Rotnitzky A.  Estimation of the effect of interventions that modify the received treatment.  _Statistics in Medicine_ 2013; 32: 5260–5277.
  * [16] Kennedy E.  Nonparametric causal effects based on incremental propensity score interventions.  _Journal of the American Statistical Association_ 2019; 114: 645–656.
  * [17] Young JG, Logan R, Robins JM et al.  Inverse probability weighted estimation of risk under representative interventions in observational studies.  _Journal of the American Statistical Association_ 2019; 114: 938–947.
  * [18] Young JG, Hernán MA and Robins JM.  Identification, estimation and approximation of risk under interventions that depend on the natural value of treatment using observational data.  _Epidemiologic Methods_ 2014; 3: 1–19.
  * [19] Díaz I, Williams N, Hoffman KL et al.  Non-parametric causal effects based on longitudinal modified treatment policies.  _arXiv preprint arXiv:200601366_ 2020; .
  * [20] Richardson TS and Robins JM.  Single world intervention graphs (swigs): A unification of the counterfactual and graphical approaches to causality.  _Center for the Statistics and the Social Sciences, University of Washington Series Working Paper_ 2013; 128: 2013.
  * [21] Rosenbaum PR and Rubin DB.  The central role of the propensity score in observational studies for causal effects.  _Biometrika_ 1983; 70: 41–55.
  * [22] Young JG, Cain LE, Robins JM et al.  Comparative effectiveness of dynamic treatment regimes: an application of the parametric g-formula.  _Statistics in Biosciences_ 2011; 3: 119–143.
  * [23] Picciotto S, Hernán MA, Page JH et al.  Structural nested cumulative failure time models to estimate the effects of interventions.  _Journal of the American Statistical Association_ 2012; 107: 886–900.
  * [24] Robins JM, Hernán MA and Siebert U.  Effects of multiple interventions.  _Comparative quantification of health risks: global and regional burden of disease attributable to selected major risk factors_ 2004; 1: 2191–2230.
  * [25] Newey WK.  The asymptotic variance of semiparametric estimators.  _Econometrica: Journal of the Econometric Society_ 1994; : 1349–1382.
  * [26] Van Der Vaart AW.  _Asymptotic statistics_ , volume 3.  Cambridge university press, 2000.
  * [27] Tsiatis A.  _Semiparametric Theory and Missing Data_.  1 ed. New York, NY: Springer, 2006.
  * [28] Molina J, Rotnitzky A, Sued M et al.  Multiple robustness in factorized likelihood models.  _Biometrika_ 2017; 104: 561–581.
  * [29] Robins JM.  Robust estimation in sequentially ignorable missing data and causal inference models.  In _Proceedings of the American Statistical Association_ , volume 1999\. Indianapolis, IN, pp. 6–10.
  * [30] Scharfstein D, Rotnitzky A and Robins J.  Adjusting for nonignorable drop-out using semiparametric nonresponse models.  _Journal of the American Statistical Association_ 1999; 94: 1096–1120.
  * [31] Robins JM.  Robust estimation in sequentially ignorable missing data and causal inference models.  In _Proceedings of the American Statistical Association_ , volume 1999\. Indianapolis, IN, pp. 6–10.
  * [32] Robins JM, Sued M, Lei-Gomez Q et al.  Comment: Performance of double-robust estimators when “inverse probability” weights are highly variable.  _Statistical Science_ 2007; 22: 544–559.
  * [33] Petersen M, Schwab J, Gruber S et al.  Targeted maximum likelihood estimation for dynamic and static longitudinal marginal structural working models.  _Journal of causal inference_ 2014; 2: 147–185.
  * [34] Lendle S, Schwab J, Petersen ML et al.  ltmle: an r package implementing targeted minimum loss-based estimation for longitudinal data.  _Journal of Statistical Software_ 2017; 81: 1–21.
  * [35] Zheng W and Van Der Laan MJ.  Asymptotic theory for cross-validated targeted maximum likelihood estimation 2010; .
  * [36] Van Der Vaart AW and Wellner JA.  Weak convergence.  In _Weak convergence and empirical processes_. Springer, 1996.  pp. 16–28.
  * [37] Van Der laan MJ, Polley EC and Hubbard AE.  Super learner.  _Statistical applications in genetics and molecular biology_ 2007; 6.
  * [38] Grant RM, Lama JR, Anderson PL et al.  Preexposure chemoprophylaxis for hiv prevention in men who have sex with men.  _New England Journal of Medicine_ 2010; 363: 2587–2599.
  * [39] Molina J, Capitant C, Spire B et al.  On-demand preexposure prophylaxis in men at high risk for hiv-1 infection.  _N Engl J Med_ 2015; 373: 2237–2246.
  * [40] McCormack S, Dunn DT, Desai M et al.  Pre-exposure prophylaxis to prevent the acquisition of hiv-1 infection (proud): effectiveness results from the pilot phase of a pragmatic open-label randomised trial.  _The Lancet_ 2016; 387: 53–60.
  * [41] Tellalian D, Maznavi K, Bredeek U et al.  Pre-exposure prophylaxis (prep) for hiv infection: results of a survey of hiv healthcare providers evaluating their knowledge, attitudes, and prescribing practices.  _AIDS patient care and STDs_ 2013; 27: 553–559.
  * [42] Krakower D, Ware N, Mitty JA et al.  Hiv providers’ perceived barriers and facilitators to implementing pre-exposure prophylaxis in care settings: a qualitative study.  _AIDS and Behavior_ 2014; 18: 1712–1721.
  * [43] Marcus JL, Hurley LB, Hare CB et al.  Disparities in uptake of hiv preexposure prophylaxis in a large integrated healthcare system.  _American journal of public health_ 2016; 106: e2.