# Heterogeneous interventional indirect effects with multiple mediators: non-
parametric and semi-parametric approaches

Max Rubinsteinlabel =<EMAIL_ADDRESS>[ Zach Bransonlabel = e2 [
<EMAIL_ADDRESS>Edward H. Kennedylabel = e3 [<EMAIL_ADDRESS>RAND Corporation  Department of Statistics & Data Science, Carnegie Mellon
University,

, , and

We propose semi- and non-parametric methods to estimate conditional
interventional indirect effects in the setting of two discrete mediators whose
causal ordering is unknown. Average interventional indirect effects have been
shown to decompose an average treatment effect into a direct effect and
interventional indirect effects that quantify effects of hypothetical
interventions on mediator distributions. Yet these effects may be
heterogeneous across the covariate distribution. We therefore consider the
problem of estimating these effects at particular points. We first propose an
influence-function based estimator of the projection of the conditional
effects onto a working model, and show that under some conditions we can
achieve root-n consistent and asymptotically normal estimates of this
parameter. Second, we propose a fully non-parametric approach to estimation
and show the conditions where this approach can achieve oracle rates of
convergence. Finally, we propose a sensitivity analysis for the conditional
effects in the presence of mediator-outcome confounding given a bounded
outcome. We propose estimating bounds on the conditional effects using these
same methods, and show that these results easily extend to allow for
influence-function based estimates of the bounds on the average effects. We
conclude by demonstrating our methods to examine heterogeneous mediated
effects with respect to the effect of COVID-19 vaccinations on depression via
social isolation and worries about health during February 2021.

##  1 Introduction

A goal of causal mediation analysis is to understand the mechanisms through
which interventions work. “Natural effects” most directly pertain to the idea
of mechanism (Miles, 2022) and decompose the individual-level treatment effect
into pathways that work directly or via changes in mediator values. However,
the identifying assumptions required to estimate these effects are
unenforceable even in randomized experiments. These effects are also not
generally identified in common applied settings that involve multiple
mediators unless the mediators are considered jointly. “Interventional
effects” were proposed as alternative causal estimands that are identifiable
under weaker assumptions and in settings with multiple mediators (see, e.g.,
Didelez, Dawid and Geneletti (2006), Vansteelandt and Daniel (2017)). These
effects conceptualize hypothetical interventions on the mediator distributions
defined at specific covariate values. Unlike natural effects, these effects
are identifiable in a sequentially randomized experiment. While the
relationship between interventional effects and mechanisms acting at an
individual-level is unclear (Miles, 2022), interventional effects have gained
popularity in applied research over the past several years.

This popularity is also in part because the same statistical functionals yield
alternative causal interpretations that are often of substantive interest.
Specifically, under weaker assumptions these same methods can quantify
disparity reductions achieved via interventions on some possibly mediating
factor(s). For example, Vansteelandt and Daniel (2017) analyze disparities in
breast-cancer survival among high and low socioeconomic status (SES) women.
They consider a model where SES causes breast-cancer survival via a direct
pathway and via cancer screening and treatment choices, so that SES takes the
role of an exposure. They estimate that if low SES women had the same
(conditional) distribution of cancer screening and treatment choices as high
SES women, the observed disparity in breast-cancer survival between high and
low SES women would be reduced by half. This effect requires conceiving of an
unspecified hypothetical intervention that could shift the covariate-stratrum
specific distributions of cancer screening and treatment choices among low SES
women to match that of high SES women. However, this intervention does not
require conceiving of potential outcomes with respect to SES, or more
generally of the types of controversial counterfactual quantities required to
conceive of natural effects (Vansteelandt and Daniel, 2017).

To date the literature on interventional effects has primarily focused on
estimating average effects. We instead consider estimating conditional effects
across covariates that are possibly continuous. For example, consider the
application from Vansteelandt and Daniel (2017). One natural follow-up
question might be how these disparity reductions change as a function of a
woman’s age. Even if the total disparity in cancer survival were constant
across age, it remains possible that age moderates the interventional effects,
and therefore also the proportion of the disparity that would be eliminated
via such an intervention. These questions pertain to conditional
interventional indirect effects (CIIE). To our knowledge, proposed strategies
to estimate the CIIE have been limited to parametric methods (see, e.g.,
Vansteelandt and Daniel (2017), Loh et al. (2020)), and the validity of the
inferences are tied to assumptions that the models are correctly specified.

Our first contribution is to propose methods that allow for flexible non-
parametric and machine learning methods for estimation. Specifically, we
consider the setting of a binary intervention and two discrete-valued
mediators whose causal ordering is unknown. We propose two estimation
procedures: first, a semi-parametric projection-based approach; second, a
fully non-parametric approach. Both procedures are conceptually simple, and
involve a regression of an estimate of the uncentered influence function for
the average effect onto the covariates. However, the semi-parametric approach
targets a projection of the CIIE onto a parametric model rather than the CIIE
itself. Projection-based estimators have frequently been proposed in the
context of different causal estimands (see, e.g., Kennedy, Lorch and Small
(2019), Cuellar and Kennedy (2020), Kennedy, Balakrishnan and Wasserman
(2021)). Our proposal extends this idea to this setting, and we show that
under some conditions, root-n consistent and asymptotically normal estimates
of the projection parameter are possible. The second proposal considers a
fully non-parametric estimation procedure (a “DR-Learner”) that targets the
CIIE directly, extending results from Kennedy (2020). While directly targeting
the CIIE instead of its projection may seem preferable, we cannot in general
obtain obtain root-n consistent estimates. Even so, we show that we can obtain
oracle rates of convergence in some settings. Both the projection estimator
and the DR-Learner substantially weaken the modeling assumptions used to date
in the literature on estimating the CIIE. Moreover, these methods allow use of
flexible non-parametric and machine-learning methods for estimation while
still obtaining relatively fast rates of convergence.

Our methods, like most, require several identifying assumptions, including
that the mediator-outcome relationship is unconfounded. A natural question is
whether our estimators are sensitive to violations of this assumption. We
therefore consider bounds on the CIIE and show that we can use both the
projection-based approach and the DR-Learner to estimate these quantities.
These methods naturally extend to allow for estimating bounds on the average
effects, and we show that root-n consistent and asymptotically normal
estimates of these bounds are possible under some conditions. Existing
approaches frequently focus on the natural rather than interventional effects,
and are often tied to strong parametric modeling assumptions (see, e.g., Park,
Qin and Lee (2020), Park and Esterling (2021), Imai, Keele and Tingley
(2010)). Moreover, sensitivity analyses for conditional estimands is less
seldom discussed (though see Lindmark, de Luna and Eriksson (2018)).

Finally, we demonstrate these methods using an application previously
considered in Rubinstein, Haviland and Breslau (2023). This study sought to
quantify the extent to which COVID-19 vaccines reduced self-reported
depression via changes in social isolation versus worries about health among
the COVID-19 Trends and Impact Survey (CTIS) respondents in February 2021. The
authors only examined effect heterogeneity across discrete subgroups;
moreover, they did not conduct a sensitivity analysis with respect to the
interventional effect estimates. We revisit this analysis and model how the
vote share for Joe Biden in the 2020 US presidential election in each
respondent’s county of residence moderated the interventional effects. We then
demonstrate our sensitivity analysis for the average and the conditional
interventional effects.

This paper proceeds as follows. In Section 2 we review interventional effects,
the required identifying assumptions, and efficient estimation. In Section 3
we introduce the CIIE and propose the projection-based estimator and the DR-
Learner. We establish conditions required for asymptotic normality and root-n
consistency of the projection estimator, and for obtaining oracle rates of
convergence for the DR-Learner. Section 4 contains a simulation study
demonstrating that these theoretic properties hold in practice. Section 5
proposes our sensitivity analysis, Section 6 contains our application, and
Section 7 contains a discussion of these results.

##  2 Review

We define the average interventional effects, the causal assumptions required
to tie the causal targets to observed data, and efficient estimation of the
observed data functionals. We largely summarize material covered in
Vansteelandt and Daniel (2017) and Benkeser and Ran (2021), and refer to those
papers for more details. We begin by outlining the setup and notation that we
will use throughout.

###  2.1 Setup and notation

Assume that we observe $n$ i.i.d. samples of observations
$Z_{i}=(V_{i},W_{i},M_{1i},M_{2i},A_{i},Y_{i})$, where $Y$ represents some
outcome of interest (either continuous or discrete), $A$ represents a binary
intervention, and $M_{1}$ and $M_{2}$ represent discrete-valued mediators.
Finally, we let $[V,W]$ represent a matrix of either discrete or continuous
covariates, where $W$ may be empty, and which we jointly denote as $X$.111We
distinguish between $V$ and $W$ because we will estimate effects conditional
on $V=v$, which may or may not include all elements of $X$. The figure below
illustrates the assumed relationships between the variables, with an arrow
indicating a causal pathway. Importantly, we do not assume that we know the
causal relationship between $M_{1}$ and $M_{2}$.

Figure 1: Assumed data generating processX = [V, W]A$M_{1}$$M_{2}$Y

While this figure helps to motivate the problem, we will primarily rely on
potential outcomes notation to define our assumptions. Specifically, we assume
the potential outcomes $Y^{am_{1}m_{2}}$ under $A=a$, $M_{1}=m_{1}$, and
$M_{2}=m_{2}$, and that $M_{j}^{a}$ represent the counterfactual outcome for
the $j$-th mediator under $A=a$. For any discrete variable $C$ we use $p(c)$
as a short-hand for $P(C=c)$ throughout.

We also define the following functions of the data. Let $\pi_{a}(X)=p(A=a\mid
X)$ and $\mu_{a}(M_{1},M_{2},X)=\mathbb{E}[Y\mid A=a,X,M_{1},M_{2}]$. We
denote the joint mediator probabilities as $p(M_{1},M_{2}\mid X,A)$, and the
marginal probabilities as $p(M_{1}\mid X,A)$ and $p(M_{2}\mid X,A)$,
respectively. Following the notation in Benkeser and Ran (2021), we define the
marginalized outcome models:

$\displaystyle\mu_{a,M_{1}}(M_{2},X)=\sum_{m_{1}}\mu_{a}(m_{1},M_{2},X)p(m_{1}\mid
a,X)$ $\displaystyle\mu_{a,M_{1}\times
M_{2}^{^{\prime}}}(X)=\sum_{m_{1},m_{2}}\mu_{a}(m_{1},m_{2},X)p(m_{1}\mid
a,X)p(m_{2}\mid a^{\prime},X)$

These quantities are defined with respect to marginalizing the outcome
regression over $p(M_{1}\mid a,X)$. We use this same notation to define
quantities marginalized with respect to $p(M_{1}\mid a^{\prime},X)$ (e.g.
$\mu_{a,M_{1}^{\prime}\times M_{2}^{^{\prime}}}(M_{1},X)$).

We also denote sample averages using $\mathbb{P}_{n}(Z)$ and regression
estimates as $\hat{\mathbb{E}}_{n}(Y\mid X)$. For possibly random functions
$f$, we denote $\|f\|^{2}$ as the squared $L_{2}(\mathbb{P})$ norm, $\int
f(z)^{2}dP(z)$. We use the notation $a\lesssim b$ to indicate that $a\leq Cb$
for some universal constant $C$. Finally, for any random function $\hat{f}$
learned on an independent sample of size $n$, $D_{1}^{n}$, we let
$\mathbb{P}\hat{f}(Z)=\int\hat{f}(z)dP(z\mid D_{1}^{n})$. That is,
$\mathbb{P}\hat{f}(Z)$ refers to the expected value of this estimated function
conditional on the training sample, noting that for a fixed function $f$ this
is equivalent to $\mathbb{E}(f)$.

###  2.2 Interventional effects

We can decompose the average effect $\psi=\mathbb{E}[Y^{a}-Y^{a^{\prime}}]$
into the sum of the following quantities (Vansteelandt and Daniel, 2017):

$\displaystyle\psi_{IDE}$
$\displaystyle=\mathbb{E}\left(\sum_{m_{1},m_{2}}\mathbb{E}[Y^{am_{1}m_{2}}-Y^{a^{\prime}m_{1}m_{2}}\mid
X]p(M_{1}^{a^{\prime}}=m_{1},M_{2}^{a^{\prime}}=m_{2}\mid X)\right)$ (1)
$\displaystyle\psi_{M_{1}}$
$\displaystyle=\mathbb{E}\left(\sum_{m_{1},m_{2}}\mathbb{E}[Y^{am_{1}m_{2}}\mid
X]\\{p(M_{1}^{a}=m_{1}\mid X)-p(M_{1}^{a^{\prime}}=m_{1}\mid
X)\\}p(M_{2}^{a^{\prime}}=m_{2}\mid X)\right)$ (2) $\displaystyle\psi_{M_{2}}$
$\displaystyle=\mathbb{E}\left(\sum_{m_{1},m_{2}}\mathbb{E}[Y^{am_{1}m_{2}}\mid
X]\\{p(M_{2}^{a}=m_{2}\mid X)-p(M_{2}^{a^{\prime}}=m_{2}\mid
X)\\}p(M_{1}^{a}=m_{1}\mid X)\right)$ (3) $\displaystyle\psi_{Cov}$
$\displaystyle=\mathbb{E}\left(\sum_{m_{1},m_{2}}\mathbb{E}[Y^{am_{1}m_{2}}\mid
X]\\{p(M_{1}^{a}=m_{1},M_{2}^{a}=m_{2}\mid X)-p(M_{1}^{a}=m_{1}\mid
X)p(M_{2}^{a}=m_{2}\mid X)\right.$ (4)
$\displaystyle\left.-[p(M_{1}^{a^{\prime}}=m_{1},M_{2}^{a^{\prime}}=m_{2}\mid
X)-p(M_{1}^{a^{\prime}}=m_{1}\mid X)p(M_{2}^{a^{\prime}}=m_{2}\mid
X)]\\}\right)$

The interventional direct effect $(\psi_{IDE})$ represents the contrast in
mean potential outcomes when we set $A=a$ for everyone in the population
versus $A=a^{\prime}$, while drawing the mediators randomly for each
individual from their counterfactual joint distribution under $A=a^{\prime}$
given subject-specific covariates $X$. By contrast the interventional indirect
effect via $M_{1}$ $(\psi_{M_{1}})$ holds $A$ fixed at $a$ for all
individuals, and considers the average contrast between giving everyone
subject-specific values of $M_{1}$ drawn randomly from the counterfactual
distribution under $A=a$ versus the distribution under $A=a^{\prime}$ given
covariates $X$, and simultaneously drawing $M_{2}$ from the counterfactual
distribution under $A=a^{\prime}$ given covariates $X$. The interventional
indirect effect through $M_{2}$ $\left(\psi_{M_{2}}\right)$ is defined
analogously. The covariant effect $(\psi_{Cov})$ is the difference between the
total effect and all three of these effects and captures the effect of the
dependence of the mediators on each other. This decomposition also holds
switching the role of $a$ and $a^{\prime}$: Vansteelandt and Daniel (2017)
indexes these estimands by $a$ to make this distinction, while for conceptual
simplicity we do not. Additionally, for the purposes of this paper we focus on
effects via $M_{1}$; however, our proposed methods can be used for any of
these estimands.

###  2.3 Identification

The estimands defined above require knowledge about the potential outcomes
under each treatment and mediator value for each subject. However, for any
individual we do not observe all of these quantities. We therefore make the
following identifying assumptions to connect these causal quantities to the
observed data distribution. First, we assume consistency, where for
$a^{\star}\in\\{a,a^{\prime}\\}$ and any $(m_{1},m_{2})$:

######  Assumption 1.

[Consistency]

$\displaystyle
A=a^{\star}\implies(M_{1},M_{2})=(M_{1}^{a^{\star}},M_{2}^{a^{\star}})$
$\displaystyle A=a^{\star},M_{1}=m_{1},M_{2}=m_{2}\implies
Y=Y^{a^{\star}m_{1}m_{2}}$

Consistency precludes the potential outcomes for any individual from depending
on another individual’s treatment or mediator assignment. We next assume
sequential ignorability:

######  Assumption 2.

[Sequential ignorability]

$\displaystyle Y^{a^{\star}m_{1}m_{2}}\perp A\mid X$ (5)
$\displaystyle(M_{1}^{a^{\star}},M_{2}^{a^{\star}})\perp A\mid X$ (6)
$\displaystyle Y^{a^{\star}m_{1}m_{2}}\perp(M_{1},M_{2})\mid(A=a^{\star},X)$
(7)

Sequential ignorability consists of three assumptions: equations (5)-(6), or
Y-A and M-A ignorability, state that $A$ is independent of $Y^{am_{1}m_{2}}$
and $(M_{1}^{a},M_{2}^{a})$ given $X$. Equation (7), or Y-M ignorability,
states that $(M_{1},M_{2})$ is independent of $Y^{am_{1}m_{2}}$ given $X$ and
$A$. We next assume positivity:

######  Assumption 3.

[Positivity]

$\displaystyle P(\min_{a^{\star}}\pi_{a^{\star}}(X)>\epsilon)=1,\qquad
P(\min_{m_{1},m_{2},a^{\star}}p(m_{1},m_{2}\mid
a^{\star},X)>\epsilon)=1,\qquad\epsilon>0$ (8)

Equation (8) implies that the propensity scores are bounded away from zero and
one and the joint mediator probabilities are bounded away from zero with
probability one.222Equation 8 is technically stronger than necessary. For
example, for any $x$, we only need that $p(m_{1},m_{2}\mid a,x)>0$ whenever
$p(m_{1},m_{2}\mid a^{\prime},x)>0$. Under these assumptions we can write
$\psi_{M_{1}}$ in terms of the observed data distribution.

$\displaystyle\psi_{M_{1}}$
$\displaystyle=\mathbb{E}\left(\sum_{m_{1},m_{2}}\mu_{a}(X,m_{1},m_{2})\\{p(m_{1}\mid
X,a)-p(m_{1}\mid X,a^{\prime})\\}p(m_{2}\mid
X,a^{\prime})\right)=\mathbb{E}\left(\psi_{M_{1}}(X)\right)$ (9)

The functional in (9) reflects other interesting causal parameters under
weaker assumptions. For example, consider the case where (5) holds but (6)
does not. This situation is frequently relevant in cases where we are using
interventional indirect effects to understand disparities and $A$ is an
indicator of some subgroup of interest (for example, Black versus White
individuals). In that case (9) still targets the causal contrast:

$\displaystyle\mathbb{E}\left(\sum_{m_{1},m_{2}}\mathbb{E}[Y^{m_{1}m_{2}}\mid
a,X][p(m_{1}\mid a,X)-p(m_{1}\mid a^{\prime},X)]p(m_{2}\mid
a^{\prime},X)\right)$ (10)

This estimand tells us about how much an intervention on the distribution of
$M_{1}$ could reduce an observed disparity in some outcome of interest
(Vansteelandt and Daniel, 2017). However, in practice we must carefully
consider the relevant conditioning sets when defining these quantities
(Jackson, 2020).

###  2.4 Estimation

Regardless of the targeted causal quantity, (9) reflects a statistical
parameter that is a fixed function of the observed data and we require methods
to estimate this quantity. One natural idea would be to estimate each function
in (9) separately, plug them into that same expression, and take the empirical
average. If we were to use correctly specified parametric models to estimate
the nuisance functions the resulting estimate would be consistent for
$\psi_{M_{1}}$ and converge at $n^{-1/2}$ rates. Unfortunately, this is
unlikely to occur in practice. We could instead estimate these functions
flexibly using non-parametric models; however, the subsequent estimator will
in general inherit the non-parametric rate of convergence of the slowest
estimated nuisance function, a rate generally slower than $n^{-1/2}$ (Kennedy,
2022). A different estimation strategy instead utilizes the so-called
“influence curve” of (9). Influence curves are important quantities related to
statistical functionals that naturally suggest efficient estimators without
parametric modeling assumptions (Kennedy, 2022).

For example, Benkeser and Ran (2021) previously showed that $\psi_{M_{1}}$ has
the efficient influence curve $\varphi(Z;\eta)-\psi_{M_{1}}$, where the
uncentered influence curve $\varphi(Z;\eta)$ takes the form

$\displaystyle\varphi(Z;\eta)$
$\displaystyle=\frac{\mathds{1}(A=a)}{\pi_{a}(X)}\frac{\\{p(M_{1}\mid
a,X)-p(M_{1}\mid a^{\prime},X)\\}p(M_{2}\mid a^{\prime},X)}{p(M_{1},M_{2},\mid
a,X)}(Y-\mu_{a}(M_{1},M_{2},X))$ (11)
$\displaystyle+\frac{\mathds{1}(A=a)}{\pi_{a}(X)}\\{\mu_{a,M_{2}^{^{\prime}}}(M_{1},X)-\mu_{a,M_{1}\times
M_{2}^{^{\prime}}}(X)\\}$
$\displaystyle-\frac{\mathds{1}(A=a^{\prime})}{\pi_{a^{\prime}}(X)}\\{\mu_{a,M_{2}^{^{\prime}}}(M_{1},X)-\mu_{a,M_{1}^{^{\prime}}\times
M_{2}^{^{\prime}}}(X)\\}$
$\displaystyle+\frac{\mathds{1}(A=a^{\prime})}{\pi_{a^{\prime}}(X)}\left(\mu_{a,M_{1}}(M_{2},X)-\mu_{a,M_{1}\times
M_{2}^{\prime}}(X)-(\mu_{a,M_{1}^{^{\prime}}}(M_{2},X)-\mu_{a,M_{1}^{^{\prime}}\times
M_{2}^{^{\prime}}}(X))\right)$ $\displaystyle+\mu_{a,M_{1}\times
M_{2}^{^{\prime}}}(X)-\mu_{a,M_{1}^{^{\prime}}\times M_{2}^{^{\prime}}}(X)$

and where

$\displaystyle\eta=[p(M_{1},M_{2}\mid a,X),p(M_{1},M_{2}\mid
a^{\prime},X),\mu_{a}(M_{1},M_{2},X),\pi_{a}(X)]$ (12)

Using $\varphi(Z)$, we can construct the so-called “one-step” estimator of
$\psi_{M_{1}}$,

$\displaystyle\hat{\psi}^{os}_{M_{1}}=\mathbb{P}_{n}[\varphi(Z;\hat{\eta})]$
(13)

This estimator involves estimating $\eta$, plugging these estimates into (11)
yielding $\varphi(Z;\hat{\eta})$, and taking the empirical mean. As shown by
Benkeser and Ran (2021), the following conditions are sufficient for this
estimator to yield root-n consistent and asymptotically normal estimates.

  1. 1.

The estimates $\hat{\eta}$ are obtained via sample-splitting

  2. 2.

$\|\varphi(Z;\hat{\eta})-\varphi(Z;\eta)\|^{2}=o_{p}(1)$

  3. 3.

$\|\hat{\eta}-\eta\|=o_{p}(n^{-1/4})$

Condition (1) can be enforced in the estimation procedure. Condition (2)
requires that the mean-squared error of the estimated influence-function
converges in probability to zero at any rate. This would require, for example,
that the propensity-scores and their estimates be bounded away from zero and
one, the joint mediator probability $p(m_{1},m_{2}\mid a,x)$ and its estimates
are bounded away from zero, and that the nuisance estimates $\hat{\eta}$ are
consistent at any rate for $\eta$. Finally, condition (3) holds for a variety
of non-parametric estimators of $\eta$ under structural assumptions on the
underlying nuisance functions: for example, on their smoothness or sparsity.

Influence-function based estimators can therefore attain $n^{-1/2}$
convergence rates without parametric modeling assumptions and allowing for
non-parametric estimates of the nuisance parameters. Intuitively, the reason
is that we can essentially ignore the nuisance estimation error – assuming
that we estimate the nuisance functions well enough (condition (3)). The
asymptotics then follow as if we had an oracle that gave us the fixed (but
unknown) function $\varphi(Z;\eta)$ and we took the average (Kennedy, 2022).
This remarkable fact occurs because the error of influence-function based
estimators is a function of the product of errors in the nuisance estimation.
By contrast, standard methods are in general linear in the nuisance
estimation. As a result, influence-function based estimators can tolerate
relatively slow rates of convergence for the nuisance estimates (e.g. rates
achieved by non-parametric methods) while still obtaining faster rates of
convergence for the estimand itself. Analogous to results shown by Kennedy
(2020) for estimating the conditional average treatment effect (CATE), we will
next show that we can also leverage $\varphi(Z;\eta)$ to achieve relatively
fast rates of convergence when estimating the CIIE.

##  3 Conditional effects

In contrast to the average effects $\psi_{M_{1}}$, we consider estimating the
effect at a given point $V=v$, recalling that $X=[W,V]$, so that $V$ is a
subset of the covariates $X$. In a slight abuse of notation,333The
conditioning in the inner expectation is more precisely written as $[W,V=v]$
rather than $X$. we define this estimand as:

$\displaystyle\psi_{M_{1}}(v)$
$\displaystyle=\mathbb{E}\left[\sum_{m_{1},m_{2}}\mathbb{E}[Y^{am_{1}m_{2}}\mid
X]\\{p(M_{1}^{a}=m_{1}\mid X)-p(M_{1}^{a^{\prime}}=m_{1}\mid
X)\\}p(M_{2}^{a^{\prime}}=m_{2}\mid X)\mid V=v\right]$ (14)

Under Assumptions (1)-(3), these parameters are identified in the observed
data as

$\displaystyle\psi_{M_{1}}(v)$
$\displaystyle=\mathbb{E}\left(\sum_{m_{1},m_{2}}\mu_{a}(X,m_{1},m_{2})\\{p(m_{1}\mid
X,a)-p(m_{1}\mid X,a^{\prime})\\}p(m_{2}\mid X,a^{\prime})\mid V=v\right)$
(15)

A natural question is how well we can estimate these effects. Noting that:

$\displaystyle\mathbb{E}[\varphi(Z;\eta)\mid V=v]=\psi_{M_{1}}(v)$

we can think of an “oracle” influence-function based estimator as providing a
benchmark for comparison for any other CIIE estimate:

$\displaystyle\hat{\mathbb{E}}_{n}[\varphi(Z;\eta)\mid V=v]$ (16)

Just as the “oracle” estimate of $\psi_{M_{1}}$ would be an empirical average
of the true influence function, we can think of (16) as a local average of the
true influence function around the point $V=v$. As long as Assumption 3 holds,
we expect that the rate of convergence of (16) provides a valid, though
possibly unachievable, lower bound for the rate of convergence for any
estimator of the CIIE. This follows from noting that the convergence rate of
this estimator is equivalent, up to constants, of replacing
$\mu_{a}(X,m_{1},m_{2})$ with the potential outcomes $Y^{am_{1}m_{2}}$ in the
expression for $\psi_{M_{1}}(X)$, and regressing this quantity onto $V$. While
the remaining discussion compares our estimators against the oracle estimate,
we expect that the oracle rates provide, in some settings, the best possible
(minimax optimal) rates of convergence.444In fact this is true for the
proposed projection estimator. We show later that this quantity is estimable
at root-n rates in some settings, as with the average effect. However, this
statement is solely conjecture for the proposed DR-Learner.

At a high-level, both ideas we propose – the projection-estimator and the DR-
Learner – substitute the estimated influence function $\varphi(Z;\hat{\eta})$
for $\varphi(Z;\eta)$ into a regression model $\hat{\mathbb{E}}_{n}$. A key
difference between these approaches is that the projection-estimator uses a
parametric model for the regression and targets a projection of the CIIE,
while the DR-Learner instead uses a non-parametric model and targets the CIIE
itself. For either approach, we show the conditions where the corresponding
oracle rates are attainable, analogous to results derived in Kennedy (2020).
Importantly, the theoretic results assume that the nuisance estimates
$\hat{\eta}$ are estimated on an independent sample from the regression
estimate $\hat{\mathbb{E}}_{n}$. We therefore refer to this as a “second-
stage” regression, reflecting the ordering of these two procedures. We discuss
both procedures more below.

###  3.1 Projection estimator

We first consider the case where the second-stage regression estimate
$\hat{\mathbb{E}}_{n}(\cdot\mid V=v)$ is given by the parametric model
$g(v;\beta)$, where $\beta$ is a finite-dimensional parameter that minimizes
some loss function. Specifically,

$\displaystyle\beta=\arg\min_{\tilde{\beta}}\mathbb{E}[w(X)\ell\\{\psi_{M_{1}}(X)-g(V;\tilde{\beta})\\}]$
(17)

Importantly, we need not assume $\psi_{M_{1}}(X)=g(V;\beta)$ for $g(v;\beta)$
to represent a meaningful target of inference. Under no assumptions about
$\psi_{M_{1}}(X)$, $\beta$ represents a population parameter that
characterizes a projection of $\psi_{M_{1}}(X)$ onto $g(V;\beta)$. For
simplicity, we focus on the case where $\beta$ minimizes the squared-error
loss ($\ell(z)=z^{2}$). The weights $w(X)$ can vary to prioritize different
parts of the covariate space when defining the projection, though in the
simplest case we can take them to be uniform. This projection, while not
necessarily representing the true parameter $\psi_{M_{1}}(v)$, can
nevertheless represent a useful summary of this parameter. The interpretation
of parametric models as defining projections is well-known, though perhaps
underappreciated, in applied statistics, and is frequently used in more
general settings when attempting to summarize characteristics of unknown data-
generating processes (see, e.g., Angrist and Pischke (2009), Buja et al.
(2019)).

To estimate this projection, we assume standard regularity conditions and
differentiate (17) with respect to $\beta$ to obtain the following moment
condition:

$\displaystyle\mathbb{E}\left[\frac{\partial
g(V;\beta)}{\partial\beta}w(X)\left\\{\psi_{M_{1}}(X)-g(V;\beta)\right\\}\right]:=\Psi(\beta;\mathbb{P})=0$
(18)

As with the average effects, our estimation approach is again based on the
influence curve of $\Psi(\beta;\mathbb{P})$. This yields an estimating-
equation stated formally in Proposition 1.

######  Proposition 1.

Under a non-parametric model, the uncentered efficient influence curve for the
moment condition $\Psi(\beta^{\star})$ at any fixed $\beta^{\star}$ is given
by

$\displaystyle\phi(Z;\beta^{\star},\eta)$ $\displaystyle=\frac{\partial
g(V;\beta^{\star})}{\partial\beta}w(X)\left(\varphi(Z;\eta)-g(V;\beta^{\star})\right)$
(19)

This then suggests the estimator $\hat{\beta}$ that satisfies:

$\displaystyle\mathbb{P}_{n}\left[\frac{\partial
g(V;\hat{\beta})}{\partial\beta}w(X)\left(\varphi(Z;\hat{\eta})-g(V;\hat{\beta})\right)\right]=0$
(20)

Theorem 1 shows the conditions required to obtain root-n consistent and
asymptotically normal parameter estimates.

######  Theorem 1.

Consider the moment condition $\mathbb{E}[\phi(Z;\beta_{0},\eta_{0})]=0$
evaluated at the true parameters $(\beta_{0},\eta_{0})$. Now consider the
estimator $\hat{\beta}$ that satisfies
$\mathbb{P}_{n}[\phi(Z;\hat{\beta},\hat{\eta})]=0$, where $\hat{\eta}$ is
estimated on an independent sample. Assume that:

  * •

The function class $\\{\phi(Z;\beta,\eta):\beta\in\mathbb{R}^{p}\\}$ is
Donsker in $\beta$ for any fixed $\eta$

  * •

$\|\phi(Z;\hat{\beta},\hat{\eta})-\phi(Z;\beta_{0},\eta_{0})\|=o_{p}(1)$

  * •

The map $\beta\to\mathbb{P}[\phi(Z;\beta,\eta)]$ is differentiable at
$\beta_{0}$ uniformly in the true $\eta$, with non-singular derivative matrix
$\frac{\partial}{\partial\beta}\mathbb{P}\\{\phi(Z;\beta,\eta)\\}\mid_{\beta=\beta_{0}}=M(\beta_{0},\eta)$,
where $M(\beta_{0},\hat{\eta})\to^{p}M(\beta_{0},\eta_{0})$

Then

$\displaystyle\hat{\beta}-\beta$
$\displaystyle=-M^{-1}[\mathbb{P}_{n}-\mathbb{P}]\phi(Z;\beta_{0},\eta_{0})+\mathcal{O}_{p}(T_{1n}+T_{2n}+T_{3n}+T_{4n})$

where

$\displaystyle
T_{1n}=\|\hat{\mu}_{a}(M_{1},M_{2},X)-\mu_{a}(M_{1},M_{2},X)\|\|\pi_{a}(X)-\hat{\pi}_{a}(X)\|$
$\displaystyle
T_{2n}=\|\hat{\mu}_{a}(M_{1},M_{2},X)-\mu_{a}(M_{1},M_{2},X)\|[\|p(M_{1},M_{2}\mid
a,X)-\hat{p}(M_{1},M_{2}\mid a,X)\|$ $\displaystyle+\|p(M_{1}\mid
a,X)-\hat{p}(M_{1}\mid a,X)\|+\|p(M_{2}\mid a^{\prime},X)-\hat{p}(M_{2}\mid
a^{\prime},X)\|$ $\displaystyle+\|p(M_{1}\mid a^{\prime},X)-\hat{p}(M_{1}\mid
a^{\prime},X)\|]$ $\displaystyle
T_{3n}=\|\pi_{a}(X)-\hat{\pi}_{a}(X)\|[\|p(M_{1}\mid a,X)-\hat{p}(M_{1}\mid
a,x)\|$ $\displaystyle+\|p(M_{1}\mid a^{\prime},X)-\hat{p}(M_{1}\mid
a^{\prime},x)\|+\|p(M_{2}\mid a^{\prime},X)-\hat{p}(M_{2}\mid
a^{\prime},X)\|]$ $\displaystyle T_{4n}=\|p(M_{2}\mid
a^{\prime},X)-\hat{p}(M_{2}\mid a^{\prime},X)\|[\|p(M_{1}\mid
a,X)-\hat{p}(M_{1}\mid a,X)\|$ $\displaystyle+\|p(M_{1}\mid
a^{\prime},X)-\hat{p}(M_{1}\mid a^{\prime},X)\|]$

Suppose further that $T_{1n}+T_{2n}+T_{3n}+T_{4n}=o_{p}(n^{-1/2})$. Then the
proposed estimator attains the non-parametric efficiency bound and is
asymptotically normal with

$\displaystyle\sqrt{n}(\hat{\beta}-\beta)\to^{d}\mathcal{N}(0,M^{-1}\mathbb{E}[\phi\phi^{\top}]M^{-\top})$
(21)

This also implies that for any fixed value of $V=v$:

$\displaystyle\sqrt{n}(g(v;\hat{\beta})-g(v;\beta))\to^{d}\mathcal{N}\left(0,\left(\frac{\partial
g(v;\beta)}{\partial\beta}\right)^{\top}M^{-1}\mathbb{E}[\phi\phi^{\top}]M^{-\top}\frac{\partial
g(v;\beta)}{\partial\beta}\right)$ (22)

######  Remark 1.

The expression $T_{1n}+T_{2n}+T_{3n}+T_{4n}$ is expressed in terms of separate
estimates for the marginal probability $p(M_{1}\mid a,X)$ and the joint
probability $p(M_{1},M_{2}\mid a,X)$. However, in practice the estimate of
$p(M_{1}\mid a,X)$ may come from a marginalized estimate of $p(M_{1},M_{2}\mid
a,X)$; similarly, the estimates of $p(M_{1}\mid a^{\prime},X)$ and
$p(M_{2}\mid a^{\prime},X)$ may come from marginalized estimates of
$p(M_{1},M_{2}\mid a^{\prime},X)$. In this case these estimates of the
marginal probabilities will inherit the same rate of convergence as the
estimate of these joint probabilities evaluated at the worst-case value of
$m_{2}$ for $p(M_{1}\mid a,X)$ (and similarly for $p(M_{1}\mid a^{\prime},X)$
and $p(M_{2}\mid a^{\prime},X)$), allowing us to simplify the second-order
expressions above. Such an estimation approach is reasonable if we believe
that these functions have the same underlying complexity. On the other hand,
if we believe that the marginal probabilities are less complex than the joint
probabilities, we may instead wish to estimate each marginal probability
separately.

Theorem 1 illustrates that if we can estimate the components of $\eta$
(defined in (12)) quickly enough $\hat{\beta}$ will be root-n consistent for
$\beta$ and asymptotically normal. It also implies that we can obtain point-
wise confidence intervals for $g(v;\hat{\beta})$ by using any consistent
estimate of the variance of (22) and standard normal quantiles. One sufficient
condition for this to hold would be that we are able to estimate all elements
of $\eta$ at rates of at least $o_{p}(n^{-1/4})$. This is attainable by many
non-parametric regression models under some conditions (Tsybakov, 2004).
Alternatively, we could allow for slower rates for some nuisance components
while requiring faster rates for others. For example, we could allow that
$\|\mu_{a}-\hat{\mu}_{a}\|=o_{p}(n^{-1/6})$ and all other components to be
estimated at $o_{p}(n^{-1/3})$. Regardless, as we saw when reviewing
estimating the average effect $\psi_{M_{1}}$ in Section 2.4, these conditions
would imply that the regression of $\varphi(Z;\hat{\eta})$ onto $g(V;\beta)$
is asymptotically equivalent to the oracle regression of $\varphi(Z;\eta)$
onto $g(V;\beta)$. Intuitively, this is because the error induced by the
nuisance estimation is decreasing at rates faster than $n^{-1/2}$. As with
estimating $\psi_{M_{1}}$, this in turn follows because the bias of
$\hat{\beta}$ is a function of the product of errors in the nuisance
estimation, as shown in the statement of Theorem 1.

###  3.2 DR-Learner

In some applications we may not be satisfied with a projection, and may
instead wish to directly estimate $\psi_{M_{1}}(v)$. We propose estimating
this quantity using a non-parametric second-stage regression model, which we
call a DR-Learner following Kennedy (2020). Algorithm 1 provides specific
details. We then analyze the DR-Learner and derive results analogous to those
found in Kennedy (2020), giving model-free error bounds for arbitrary first-
stage estimators which reveal that under some conditions, the DR-Learner is as
efficient as an oracle estimator that regresses $\varphi(Z;\eta)$ onto $V$
directly. Our results are similar to Kennedy (2020): the primary difference is
that we must consider the sums of products of errors between several more sets
of nuisance functions.

######  Algorithm 1.

Let $(D_{1}^{n},D_{2}^{n})$ denote two independent samples of $n$ observations
of $Z_{i}$.

  * •

Step 1: Nuisance training. Construct estimates of $\eta$ using $D_{1}^{n}$

  * •

Step 2: Pseudo-outcome regression. Construct the pseudo-outcomes
$\varphi(Z;\hat{\eta})$ and regress it onto covariates $V$ in the test sample
$D_{2}^{n}$, giving

$\displaystyle\hat{\psi}_{M_{1}}^{dr}(v)=\hat{\mathbb{E}}_{n}\\{\varphi(Z;\hat{\eta})\mid
V=v\\}$ (23)

  * •

Step 3: Cross-fitting (optional). Repeat Steps 1-2, swapping the roles of
$D_{1}^{n}$ and $D_{2}^{n}$. Use the average of the resulting two estimates as
a final estimate of $\psi_{M_{1}}(v)$.

######  Remark 2.

In practice, when implementing either the DR-Learner or the projection
estimator, one may wish to use sample-split estimates of $\eta$ and regress
them onto $V$ using the entire sample, rather than averaging two separate
estimates as suggested in Step 3 above. However, our results do not provide
theoretic guarantees for this approach.

Proposition 1 from Kennedy (2020) establishes general conditions where the
error of a pseudo-outcome regression of $\hat{f}$ onto $V$ and an oracle
regression of $f$ onto $V$ are asymptotically equivalent. This result relies
on an assumption on the “stability” of the second-stage estimator with respect
to a distance measure $d$ and the convergence in probability of $\hat{f}$ to
$f$ with respect to $d$. Intuitively, estimator stability requires that the
second-stage error between the pseudo-outcome regression and the oracle
estimator converges in probability to the conditional bias of the pseudo-
outcome estimates at a rate determined by root-mean squared error of the
oracle estimator.555We provide the formal definition of stability in Appendix
A.1. Theorem 2 of Kennedy (2020) shows the conditions for the oracle
efficiency of a DR-Learner for the CATE under a direct application of
Proposition 1 from Kennedy (2020), where $f$ is the uncentered influence
function for the average treatment effect (ATE). We use the same approach for
the CIIE to obtain analogous results, formalized in Corollary 1.

######  Corollary 1.

Define $\hat{b}^{\star}(x)=\mathbb{E}\\{\hat{\varphi}(Z)-\varphi(Z)\mid
D_{1}^{n},X=x\\}$; in other words, $\hat{b}^{\star}(x)$ is the conditional
bias of the estimated influence function at $X=x$ conditional on the training
data. Assume

  1. 1.

$\hat{\mathbb{E}}_{n}$ is stable with respect to distance metric $d$

  2. 2.

$d(\hat{\varphi},\varphi)\to^{p}0$

Let $\tilde{\psi}_{M_{1}}(v)$ denote an oracle estimator from a regression of
the true efficient influence function onto $V$ and $K_{n}^{\star}(v)$ denote
the oracle root mean square error,
$\mathbb{E}[\\{\tilde{\psi}_{M_{1}}(v)-\psi_{M_{1}}(v)\\}^{2}]^{1/2}$. Then

$\displaystyle\hat{\psi}^{dr}_{M_{1}}(v)-\tilde{\psi}_{M_{1}}(v)$
$\displaystyle=\hat{\mathbb{E}}_{n}\\{\hat{b}^{\star}(X)\mid
V=v\\}+o_{p}\left(K_{n}^{\star}(v)\right)$ (24)

We provide an expression for $\hat{b}^{\star}(x)$ in Appendix B. Moreover,
$\hat{b}^{\star}(x)\lesssim\hat{b}(x)$, where:

$\displaystyle\hat{b}(x)$
$\displaystyle=T_{1n}(x)+T_{2n}(x)+T_{3n}(x)+T_{4n}(x)$ (25)

and

$\displaystyle T_{1n}(x)$
$\displaystyle=(\hat{\pi}_{a}(x)-\pi_{a}(x))\sum_{m_{1},m_{2}}(\mu_{a}(m_{1},m_{2},x)-\hat{\mu}_{a}(m_{1},m_{2},x))$
$\displaystyle T_{2n}(x)$
$\displaystyle=\sum_{m_{1},m_{2}}(\mu_{a}(m_{1},m_{2},x)-\hat{\mu}_{a}(m_{1},m_{2},x))(p(m_{1},m_{2}\mid
a,x)-\hat{p}(m_{1},m_{2}\mid a,x)$ $\displaystyle+(p(m_{1}\mid
a,x)-\hat{p}(m_{1}\mid a,x))+(p(m_{1}\mid a^{\prime},x)-\hat{p}(m_{1}\mid
a^{\prime},x))+(p(m_{2}\mid a^{\prime},x)-\hat{p}(m_{2}\mid a^{\prime},x)))$
$\displaystyle T_{3n}(x)$
$\displaystyle=(\pi_{a}(x)-\hat{\pi}_{a}(x))(p(m_{1}\mid
a,x)-\hat{p}(m_{1}\mid a,x))+(p(m_{1}\mid a^{\prime},x)-\hat{p}(m_{1}\mid
a^{\prime},x))$ $\displaystyle+(p(m_{2}\mid a^{\prime},x)-\hat{p}(m_{2}\mid
a^{\prime},x)))$ $\displaystyle T_{4n}(x)$
$\displaystyle=\sum_{m_{1},m_{2}}(p(m_{2}\mid a^{\prime},x)-\hat{p}(m_{2}\mid
a^{\prime},x))((p(m_{1}\mid a,x)-\hat{p}(m_{1}\mid a,x))+p(m_{1}\mid
a^{\prime},x)-\hat{p}(m_{1}\mid a^{\prime},x))$

Therefore, the DR-Learner is oracle efficient if
$\hat{\mathbb{E}}_{n}\\{\hat{b}(X)\mid
V=v\\}=o_{p}\left(K_{n}^{\star}(v)\right)$.

######  Remark 3.

As with the second-order expression in Theorem 1, the expression for
$\hat{b}(x)$ depends on products of errors with the marginal mediator
probabilities; however, these estimates, as described in Algorithm 1, come
from estimates of the joint mediator probabilities. The error in the marginal
probabilities is therefore of the same order as the sums of the errors in the
joint probabilities across the relevant mediator values. As noted in Remark 1,
we may wish to estimate the marginal mediator probabilities separately if we
believe that the underlying complexity of the marginal probabilities is
simpler than the joint mediator probabilities, though such a situation may
seem unlikely to occur in practice.

One interesting implication of Corollary 1 is that the rate of convergence is
a function of the cardinality of the joint mediator values $k$. We may
eliminate this dependence by invoking the following assumption:

######  Assumption 4.

The smoothed product of errors between mediator probabilities and/or the
outcome model are of the same order for any values of ($m_{1},m_{2}$). For
example, for any $(m_{1},m_{1}^{\prime})$ and $(m_{2},m_{2}^{\prime})$:

$\displaystyle\hat{\mathbb{E}}_{n}\\{\left(p(m_{1}\mid a,X)-\hat{p}(m_{1}\mid
a,X)\right)\left(p(m_{2}\mid a^{\prime},X)-\hat{p}(m_{2}\mid
a^{\prime},X)\right)\mid V=v\\}=o_{p}(a_{n})$

and

$\displaystyle\hat{\mathbb{E}}_{n}\\{\left(p(m_{1}^{\prime}\mid
a,X)-\hat{p}(m_{1}^{\prime}\mid a,X)\right)\left(p(m_{2}^{\prime}\mid
a^{\prime},X)-\hat{p}(m_{2}^{\prime}\mid a^{\prime},X)\right)\mid
V=v\\}=o_{p}(a_{n})$

Assumption 4 would be reasonable if we do not believe the form of functional
form of the mediator probabilities or outcome models varies in underlying
complexity across different values of the mediators.

We next consider the form of the second-stage regression
$\hat{\mathbb{E}}_{n}$. Proposition 2 from Kennedy (2020) implies that when
$\hat{\mathbb{E}}_{n}$ is a linear smoother of the form
$\sum_{i}w_{i}(v;V^{n})f(Z_{i})$ and $\sum_{i}\lvert
w_{i}(v;V^{n})\rvert=\mathcal{O}_{p}(a_{n})$, and $\hat{b}(x)$ can be
expressed in the form of $\hat{b}_{1}(x)\hat{b}_{2}(x)$, then

$\displaystyle\hat{\mathbb{E}}_{n}\\{\hat{b}(X)\mid
V=v\\}=\mathcal{O}_{p}(a_{n}\|\hat{b}_{1}\|_{w,2}\|\hat{b}_{2}\|_{w,2})$

where

$\displaystyle\|f\|_{w,2}=\left[\sum_{i}\left\\{\frac{\lvert
w_{i}(v;V^{n})\rvert}{\sum_{j}\lvert w_{j}(v;V^{n})\rvert}\right\\}\lvert
f(Z_{i})\rvert^{2}\right]^{1/2}$

Corollary 2 applies this result to Corollary 1.

######  Corollary 2.

Assume the conditions of Corollary 1 and that $\hat{\mathbb{E}}_{n}$ is a
minimax optimal linear smoother with $\sum_{i}\lvert
w_{i}(v;V^{n})\rvert=\mathcal{O}_{p}(1)$. Notice that

$\displaystyle\hat{b}(x)=\sum_{j}\hat{b}_{j1}(x)\hat{b}_{j2}(x)$

where $j$ indexes all of the error products in (25). Therefore:

$\displaystyle\hat{\mathbb{E}}_{n}\\{\hat{b}^{\star}(X)\mid V=v\\}$
$\displaystyle\lesssim\hat{\mathbb{E}}_{n}\\{\sum_{j}\hat{b}_{j1}(X)\hat{b}_{j2}(X)\mid
V=v\\}$

By Proposition 2 of Kennedy (2020) we then obtain that

$\displaystyle\hat{\mathbb{E}}_{n}\\{\hat{b}(X)\mid V=v\\}$
$\displaystyle=\sum_{j}\mathcal{O}_{p}(\|\hat{b}_{j1}\|_{w,2}\|\hat{b}_{j2}\|_{w,2})\asymp\max_{j}\mathcal{O}_{p}(\|\hat{b}_{j1}\|_{w,2}\|\hat{b}_{j2}\|_{w,2})$

To make Corollary 2 concrete consider the case where
$K_{n}^{\star}(v)=n^{-\theta}$. This result implies that when the second-stage
regression estimator is a linear smoother, the DR-learner will achieve the
corresponding oracle rate when, for example, all of nuisance errors are at
least $o_{p}(n^{-\theta/2})$ in the $\|\cdot\|_{w,2}$ norm. More generally,
the DR-Learner is oracle efficient as long as the highest order error product
in (25) is $o_{p}(n^{-\theta})$ in this same norm. Whether these rates are
actually attainable in a given application depends on the underlying
complexity of the nuisance functions.

##  4 Simulations

We verify that the expected performance of these estimation strategies
corresponds with the theory outlined above using a simulation study. First, we
outline the data generating process; second, we demonstrate the performance of
our proposed approaches on samples of size $n=1000$ when estimating the
nuisance functions using SuperLearner (Van der Laan, Polley and Hubbard,
2007); finally, we compare the convergence rates of the DR-Learner versus a
plugin approach while controlling the convergence rates of the nuisance
estimation.

###  4.1 Setup

To illustrate our proposed approach, we conduct a simulation study with a one-
dimensional covariate $X$ (so that $V=X$). Figure 2 illustrates the simulated
nuisance functions as a function of $X$. One aspect of this setup is that the
outcome models and the propensity score models have complexity unlikely to be
fully captured by simple generalized linear models, motivating our use of non-
parametric methods. A second aspect is that the implied function
$\psi_{M_{1}}(X)$, illustrated in Figure 3, is less complex than these
functions. We provide all details about the data generating processes for our
simulations, including the functions illustrated in Figure 2, in Appendix C.

Figure 2: | Simulation: selected nuisance functions. Nuisance function specifications for simulation study

Figure 3 also illustrates the implied curves of $\psi_{M_{2}}(X)$ and
$\psi(X)$, as well as the curves for the proportion mediated via each mediator
(e.g. $\psi_{M_{1}}(X)/\psi(X)$). The effects are entirely mediated via
$M_{1}$ and $M_{2}$,666Additionally, the covariant effects due to the
dependence of the mediators on each other is close to zero throughout. and the
proportion mediated via $M_{1}$ increases with $X$.

Figure 3: | Estimands. Total effects, indirect effects, and proportion mediated as a function of $X$

###  4.2 Estimation: SuperLearner

We evaluate the performance of each estimator across 1000 simulations on test
samples of size 1000. To estimate the nuisance parameters we use SuperLearner,
using both the “SL.glm” and “SL.ranger” libraries.777In practice it is
desirable to use as many libraries as possible when running SuperLearner;
however, for the sake of computation time we only use these two libraries for
our simulation study. These model our nuisance functions as a weighted
combination of predictions from logistic regression and random forests models.
After estimating the nuisance parameters using samples of size 1000, we use a
separate test sample to estimate the DR-Learner and projection estimators. We
then predict the points at $X=0$ and $X=2$.

While we expect both estimates to be consistent, the mean square error at each
point should differ by constants: this is due to differing inverse weights in
the expression for $\varphi(Z)$. Figure 8 in Appendix 12 illustrates the
maximum possible inverse weight as a function of X: at $X=0$, the maximum
weight is lowest while at $X=2$ the maximum weight is highest. These two
points arguably reflect the easiest and hardest parts of the covariate space
to estimate, with the point where $X=2$ reflecting a “worst-case scenario” in
our simulation. As long as our assumptions hold, these weights do not affect
the asymptotic results. However, they can affect their performance in finite
samples, with higher variance estimates where the inverse weights are large.
Additionally, confidence intervals may have under-coverage, since their
validity is also based on asymptotic approximations. Consequently, we expect
the simulations to show better results when estimating the CIIE at $X=0$
compared to $X=2$, also with possibly better coverage in these regions. More
generally, this comparison highlights a key limitation of our proposed
approach: in an actual sample it may be challenging to estimate conditional
effects where the inverse weights are extreme.

Table 1: Projection estimators: simulation performance, $n=1000$
(averaged over 1000 simulations) Point | Strategy | Projection | Truth | Bias | Std | RMSE | Coverage
---|---|---|---|---|---|---|---
0 | Plugin | Linear | 0.07 | 0.00 | 0.03 | 0.03 | 3.10
2 | Plugin | Linear | 0.11 | -0.04 | 0.02 | 0.04 | 1.00
0 | Plugin | Quadratic | 0.07 | -0.00 | 0.03 | 0.03 | 3.10
2 | Plugin | Quadratic | 0.11 | -0.04 | 0.02 | 0.05 | 0.80
0 | Efficient | Linear | 0.07 | 0.00 | 0.06 | 0.06 | 95.50
2 | Efficient | Linear | 0.11 | -0.00 | 0.08 | 0.08 | 93.80
0 | Efficient | Quadratic | 0.07 | 0.00 | 0.06 | 0.06 | 95.10
2 | Efficient | Quadratic | 0.11 | 0.00 | 0.10 | 0.10 | 93.50

Table 1 considers the projection estimates, and displays the bias, RMSE, and
confidence interval coverage associated with our proposed approach
(“Efficient”) and with a plugin approach that regresses plugin estimates of
$\psi_{M_{1}}(x)$ on the same model. Specifically, the plugin estimates
involve estimating each component in the expression for $\psi_{M_{1}}(x)$ and
regressing these estimates, rather than estimates of $\varphi(Z;\eta)$, onto
$g(X;\beta)$. We predict the projection at the points $X\in\\{0,2\\}$ using
either a linear or quadratic projection. While the plugin estimator has lower
RMSE, the confidence interval coverage is close to zero. This reflects that
the bias associated with the nuisance estimation does not converge quickly
enough to zero, and we therefore cannot ignore the estimation error in the
second-stage regression when conducting inference. By contrast, we obtain
close to nominal coverage rates for the efficient estimator. Finally, as
expected, the point $X=0$ is easier to estimate than $X=2$, reflected by the
fact that the RMSE is higher for estimates at $X=2$ than $X=0$.

Table 2 displays analogous results when using the DR-Learner to target the
true CIIE rather than its projection.888Because the implied curves are
relatively easy to approximate using a linear or quadratic model, we see that
the “Truth” column in Table 2 is identical to the “Truth” column in Table 1.
In fact these are only identical rounded to the second decimal place, but this
illustrates that the projections provide good estimates of the true quantities
in our simulation. We use smoothing splines with the default tuning parameters
for the second-stage regression,999The tuning parameters are chosen by default
using generalized cross-validation. We technically should account for post-
selection inference; however, this does not seem to make a difference in
practice in our simulations. and use the variance estimates from the smoothing
matrix and assume that the distribution of the estimates is asymptotically
normal to generate confidence intervals. Table 2 shows that this procedure
yields approximately nominal coverage rates.

Table 2: Non-parametric estimators: simulation performance, $n=1000$
(averaged over 1000 simulations) Point | Strategy | Truth | Bias | Std | RMSE | Coverage
---|---|---|---|---|---|---
0 | DR-Learner | 0.07 | 0.00 | 0.08 | 0.08 | 94.3
0 | Plugin | 0.07 | -0.00 | 0.03 | 0.03 | -
2 | DR-Learner | 0.11 | -0.00 | 0.11 | 0.11 | 92.6
2 | Plugin | 0.11 | -0.04 | 0.03 | 0.05 | -

As with the projection approach, the corresponding plugin approach has lower
RMSE than the DR-Learner. This is likely a function of the inverse-probability
weights associated with the DR-Learner, which could cause this result for a
fixed sample size.

###  4.3 Estimation: convergence rates

We conclude by examining the convergence rates of the DR-Learner versus a
plugin estimator by specifying the convergence rates of the nuisance
estimation. Roughly, we add $\mathcal{N}(C/n^{\alpha},C/n^{2\alpha})$ to the
true values of the nuisance parameters on the logistic scale to simulate
estimates. Using these results we construct “estimated” influence functions
and regress them onto $X$ using smoothing splines. In contrast to the
simulation study above, we use these simulations to estimate the integrated
mean square error across the entire domain of $X$. Figure 4 displays the
results.

Figure 4: | Convergence of DR-Learner versus Plugin and Oracle estimators. Scaled RMSE of each estimation strategy as a function of sample size. “Slow” nuisance function is estimated $\mathcal{O}_{p}(n^{-1/10})$ rates, remaining at $\mathcal{O}_{p}(n^{-1/2})$ rates.

The y-axis displays the RMSE of each estimator scaled by $\sqrt{n}$, while the
x-axis displays the sample sizes. The top left panel considers the case where
we set $\alpha=0.5$ for all nuisance parameters. The other panels instead set
$\alpha=0.1$ for the function indicated in the panel title. As expected, when
all nuisance functions are estimated at the same rate the plugin estimator
converges at the same rate as the DR-Learner. However, once one of the
nuisance functions is estimated slowly, the plugin estimator converges more
slowly (illustrated by the diverging green lines), while the DR-Learner
appears to attain the oracle rates of convergence in all settings. We again
observe that despite the slower convergence rates, the plugin estimator has
lower RMSE than either the DR-Learner or oracle estimators in some settings.

In Appendix C.3 we present additional results where we estimate the CIIE as a
proportion of the corresponding CATE. We outline two general approaches to
this problem: first, where we estimate the CIIE and the CATE and take the
ratio of these estimates; second, where we derive the influence function for
the mean of the ratio and regress this onto $V$ (this is similar to Cuellar
and Kennedy (2020), who consider estimating a conditional risk-ratio). Our
simulations show that this quantity is quite difficult to estimate well due to
the high-variance of the estimators, although we are able to construct
confidence intervals with approximately nominal coverage rates using either
approach.

##  5 Sensitivity analysis

We consider estimating $\psi_{M_{1}}(v)$ when the assumption that the
potential outcomes $Y^{am_{1}m_{2}}$ are independent of the mediators given
the covariates and that $A=a$ does not hold. This might occur, for example, if
there were a post-treatment confounder $L$ that occurs prior to $M$ but after
$A$; or, a pre-treatment confounder $U$ that affects the Y-M relationship but
not the A-M or A-Y relationships. We first outline a general sensitivity
framework to generate bounds on the conditional or average effects while
specifying the degree of these types of violations, where our approach builds
from work in Luedtke, Diaz and van der Laan (2015). We extend the projection
estimator and DR-Learner to estimate bounds on the conditional effects, though
our proposed method naturally also suggests influence function based
estimators of the bounds on the average effects. These analyses can help an
analyst assess how much inferences may change given a specified degree of
confounding. We first briefly introduce additional assumptions and notation to
ease exposition.

###  5.1 Setup and notation

First, we assume for simplicity that $V=X$; that is, that the conditioning set
is identical to the observed confounders, so that our target estimand is
$\psi_{M_{1}}(x)$. Second, to construct a bound for $\psi_{M_{1}}(x)$, it will
be helpful to write
$\psi_{M_{1}}(x)=\psi_{M_{1},a}(x)-\psi_{M_{1},a^{\prime}}(x)$, where:

$\displaystyle\psi_{M_{1},a}(x)$
$\displaystyle=\sum_{m_{1},m_{2}}\mathbb{E}[Y^{m_{1}m_{2}}\mid
a,x][p(m_{1}\mid a,x)]p(m_{2}\mid a^{\prime},x)$
$\displaystyle\psi_{M_{1},a^{\prime}}(x)$
$\displaystyle=\sum_{m_{1},m_{2}}\mathbb{E}[Y^{m_{1}m_{2}}\mid
a,x][p(m_{1}\mid a^{\prime},x)]p(m_{2}\mid a^{\prime},x)$

Third, for any $(m_{1}^{\prime},m_{1},m_{2}^{\prime},m_{2})$ we let

$\displaystyle\mathbb{E}[Y^{m_{1}m_{2}}\mid
m_{1}^{\prime},m_{2}^{\prime},a,x]=\mu_{am_{1}m_{2}}^{\star}(m_{1}^{\prime},m_{2}^{\prime},x)$

Assuming that Y-M ignorability holds when it does not, we define the biased
target of our estimator of $\psi_{M_{1}}(x)$:

$\displaystyle\bar{\psi}_{M_{1}}(x)$
$\displaystyle=\bar{\psi}_{M_{1},a}(x)-\bar{\psi}_{M_{1},a^{\prime}}(v)$
$\displaystyle=\sum_{m_{1},m_{2}}\mu_{a}(m_{1},m_{2},x)[p(m_{1}\mid
a,x)-p(m_{1}\mid a^{\prime},x)]p(m_{2}\mid a^{\prime},x)$

Finally, to ease notation, we let $(\cdot)$ indicate the arguments
$(m_{1},m_{2},x)$.

While we focus the remaining discussion on the case where $V=X$, all of these
results also hold at a point $V=v$ by averaging the relevant quantities over
the distribution $p(W\mid V=v)$, recalling that $X=[V,W]$.

###  5.2 Sensitivity framework

The bias $\bar{\psi}_{M_{1}}(x)-\psi_{M_{1}}(x)$ occurs because for any
$(m_{1},m_{2},x)$, in general $\mathbb{E}[Y^{m_{1}m_{2}}\mid
a,x]\neq\mu_{a}(\cdot)$. We first consider this bias of the outcome
regression. Proposition 2 shows that we can bound this bias as a function of
$\mu_{a}(\cdot)$ and a sensitivity parameter $\tau$. We consider a general
framework where the meaning of $\tau$ changes based on the chosen sensitivity
analysis; however, we describe the interpretation of this parameter under each
assumption below.

######  Proposition 2.

Assume that we know some functions $b_{l}(\cdot;\tau)$ and $b_{u}(\cdot;\tau)$
parameterized by $\tau$ such that for every $(m_{1},m_{2},x)$:

$\displaystyle
b_{u}(\cdot;\mu_{a},\tau)\geq\mu_{am_{1}m_{2}}^{\star}(M_{1}\neq
m_{1},M_{2}\neq m_{2},x)-\mu_{a}(\cdot)\geq b_{l}(\cdot;\mu_{a},\tau)$

This implies the following bounds:

$\displaystyle b_{u}[\cdot;\mu_{a},\tau][1-p(m_{1},m_{2}\mid
a,x)]\geq\mathbb{E}[Y^{m_{1}m_{2}}\mid a,x]-\mu_{a}(\cdot)\geq
b_{l}[\cdot;\mu_{a},\tau][1-p(m_{1},m_{2}\mid a,x)]$

Different assumptions on the selection process can motivate different
functions $b_{l}$ and $b_{u}$. For example, let $\tau(m_{1},m_{2},x)\in[0,1]$.
Consider the following three assumptions for any $(m_{1}^{\prime}\neq m_{1})$
and $(m_{2}^{\prime}\neq m_{2})$:

$\displaystyle\tau(\cdot)\geq\rvert\mu_{am_{1}m_{2}}^{\star}(m_{1}^{\prime},m_{2}^{\prime},x)-\mu_{a}(\cdot)\lvert$
(26)
$\displaystyle\tau(\cdot)+\mu_{a}(\cdot)[1-\tau(\cdot)]\geq\mu_{am_{1}m_{2}}^{\star}(m_{1}^{\prime},m_{2}^{\prime},x)\geq(1-\tau(\cdot))\mu_{a}(\cdot)$
(27) $\displaystyle
1/(1-\tau(\cdot))\geq\mu_{am_{1}m_{2}}^{\star}(m_{1}^{\prime},m_{2}^{\prime},x)/\mu_{a}(\cdot)\geq(1-\tau(\cdot))$
(28)

Under assumption 26, $\tau$ bounds the absolute value of the difference
between the regression functions
$\mu_{am_{1}m_{2}}^{\star}(m_{1}^{\prime},m_{2}^{\prime},x)$ and
$\mu_{a}(\cdot)$ for each level of $(m_{1},m_{2},x)$. Under assumption 28,
$\tau$ parameterizes deviations of these same regression functions on the
risk-ratio scale: below by $1-\tau$, and above by
$\frac{1}{1-\tau}$.101010Often a sensitivity parameter $\Gamma$, which equals
$\frac{1}{1-\tau}$, is used instead in this framework. We choose $\tau$ here
to maintain consistency with the other two possible approaches. Finally, the
meaning of $\tau$ changes for the upper and lower bound under assumption 27.
First, the lower bound is equivalent to the lower bound in assumption 28, and
$\tau$ retains an equivalent meaning in this case. However, the upper bound
instead specifies that
$(1-\mu_{am_{1}m_{2}}^{\star}(m_{1}^{\prime},m_{2}^{\prime},x))/(1-\mu_{a}(\cdot))\geq(1-\tau(\cdot))$.
Under this assumption $\tau$ parameterizes the risk-ratio of the regression
function when the event $Y$ did not occur.111111This bound can be used when
$Y$ is binary, but more generally when $Y$ is bounded and rescaled to fall
within zero and one. We discuss the trade-offs between these assumptions in
greater detail below.

While these assumptions provide bounds on the true outcome model, they imply,
but are not equivalent to, bounds on $\psi_{M_{1}}(x)$. Proposition 3 provides
a generic form of these bounds.

######  Proposition 3.

Consider assumptions (26)-(28) and a sensitivity parameter $\tau(x)$ that is
valid for any value of $(m_{1},m_{2})$ at the point $X=x$. All $[b_{l},b_{u}]$
implied by these assumptions can be expressed as:

$\displaystyle[b_{l},b_{u}]=[(c_{l}\mu_{a}(\cdot)+t_{l})f_{l}(\tau(x)),(c_{u}\mu_{a}(\cdot)+t_{u})f_{u}(\tau(x))]$

for constants $(c_{l},c_{u},t_{l},t_{u})\in\\{0,1\\}^{4}$ and functions
$f_{l}$ and $f_{u}$. At a point $X=x$, this yields the following bounds on
$\psi_{M_{1}}(x)$:

$\displaystyle\psi_{M_{1},ub}(x;\tau)=$ (29)
$\displaystyle[\bar{\psi}_{M_{1}}(x)+\bar{\psi}_{M_{1},a}(x)f_{u}(\tau)c_{u}-\bar{\psi}_{M_{1},a^{\prime}}(x)f_{l}(\tau)c_{l}+t_{u}f_{u}(\tau)-t_{l}f_{l}(\tau)$
$\displaystyle-
f_{u}(\tau)\sum_{m_{1},m_{2}}[c_{u}\mu_{a}(m_{1},m_{2},x)+t_{u}]p(m_{1},m_{2}\mid
a,x)p(m_{1}\mid a,x)p(m_{2}\mid a^{\prime},x)$
$\displaystyle+f_{l}(\tau)\sum_{m_{1},m_{2}}[c_{l}\mu_{a}(m_{1},m_{2},x)+t_{l}]p(m_{1},m_{2}\mid
a,x)p(m_{1}\mid a^{\prime},x)p(m_{2}\mid a^{\prime},x)]$
$\displaystyle\psi_{M_{1},lb}(x;\tau)=$ (30)
$\displaystyle[\bar{\psi}_{M_{1}}(x)+\bar{\psi}_{M_{1},a}(x)f_{l}(\tau)c_{l}-\bar{\psi}_{M_{1},a^{\prime}}(x)f_{u}(\tau)c_{u}+t_{l}f_{l}(\tau)-t_{u}f_{u}(\tau)$
$\displaystyle-
f_{l}(\tau)\sum_{m_{1},m_{2}}[c_{l}\mu_{a}(m_{1},m_{2},x)+t_{l}]p(m_{1},m_{2}\mid
a,x)p(m_{1}\mid a,x)p(m_{2}\mid a^{\prime},x)$
$\displaystyle+f_{u}(\tau)\sum_{m_{1},m_{2}}[c_{u}\mu_{a}(m_{1},m_{2},x)+t_{u}]p(m_{1},m_{2}\mid
a,x)p(m_{1}\mid a^{\prime},x)p(m_{2}\mid a^{\prime},x)]$

######  Remark 4.

If we desired bounds at the point $V=v$, we could choose a $\tau$ valid for
any realization of $W$ at the point $V=v$ and average these expressions over
the conditional distribution of $W$ given $V=v$. If we desired bounds on the
average effect, we could choose a $\tau$ valid for all $(x,m_{1},m_{2})$ and
marginalize the above expressions over the entire covariate distribution.

The assumptions outlined in equations (26)-(28) yield different bounds. While
(26) is perhaps most intuitive, for a given $\tau$ the width of the implied
bounds on $\psi_{M_{1}}$ in equations (29) and (30) are twice as large as
those from (27), and are thus perhaps less useful in practice than the others.
Comparing the assumptions in equations (27) and (28) is difficult: for a fixed
$\tau$ the scale of the confounding for the upper bounds of
$\mu_{am_{1}m_{2}}^{\star}$ in these equations is simply different. One
benefit of (28) is that it only requires reasoning about the risk-ratio
$\mu_{am_{1}m_{2}}^{\star}(m_{1}^{\prime},m_{2}^{\prime},x)/\mu_{a}(\cdot)$.
By contrast (27) requires additional reasoning about the risk-ratio bias in
the estimates that the event $Y$ did not occur, demanding more thought from
the user. On the other hand, when using a binary outcome, equation (28) may
result in an upper bound on $\mu_{am_{1}m_{2}}^{\star}$ greater than one,
while equation (27), and the resulting bound on $\psi_{M_{1}}(x)$, will always
respect the parameter space. Finally, for a fixed value of $\tau$ it is
unclear whether the bounds yielded by equations (27) or (28) will be wider.
However, for fixed $(m_{1},m_{2},x)$, the upper bound on
$\mu_{am_{1}m_{2}}^{\star}$ in (28) will always be larger than the bounds in
(27) when $\mu_{a}(\cdot)\geq 0.5$. Heuristically, we may therefore expect the
bounds yielded by (28) to be narrower than those from (27) when $\mu_{a}$
tends to be small across values of $(m_{1},m_{2},x)$. As a final point, the
bound given by (27) is only useful for binary outcomes, or bounded outcomes
rescaled to fall within 0 and 1, so that the sensitivity analysis given by
(28) is more general.

Finally, we can extend this general approach in several ways. For example,
Assumptions (26)-(28) yield similar bounds for $\psi_{M_{2}}$, $\psi_{Cov}$,
$\psi_{IDE}$ by averaging over the relevant distributions. We discuss these
extensions in Appendix E. We could also specify $\tau$ as a function that
varies across $(m_{1},m_{2},x)$ to arrive at a slightly different expression
for the bounds. However, specifying this function would be challenging in
practice.

###  5.3 Illustration

Figure 5 uses simulated data to illustrate the estimand, the biased target,
and the bounds as a function of $x$ choosing $\tau=0.1$ under (27) and
$\tau=0.15$ under (28). These parameters reflect the true maximal values of
$\tau$ guaranteed to hold under these assumptions in our simulation. We obtain
biased estimates for our outcome model based on (27) and a parameter $\tau$
that varies between $0.1*\\{1/3,2/3,1\\}$ depending on the value of $x$. We
describe the selection process in greater detail in Appendix C.4. The upper
and lower bounds are depicted in purple and red, while the orange and green
lines reflects $\bar{\psi}_{M_{1}}(x)$ and $\psi_{M_{1}}(x)$, respectively. As
$x$ increases, the bounds given by (28) are at first narrower and eventually
wider than the bounds given by (27). This is generally expected as the values
of $\mu_{a}(m_{1},m_{2},x)$ tend to increase with $x$ (see Figure 10 in
Appendix C). These bounds are also quite conservative: assuming (27), the
bounds are only guaranteed to hold for all $x$ for $\tau=0.1$. However, even
$\tau=0.02$ provides valid upper and lower bounds across the entire domain in
our simulation. Of course, the gap between the value of $\tau$ guaranteed to
hold and the minimum $\tau$ that actually does may be smaller for other data
distributions; however, it does suggest that these bounds may be conservative
in practice.

Figure 5: | Bounds on CIIE. Target estimand in green, biased target of inference in orange. Purple and red lines reflect upper and lower bounds that differ in terms of $\tau$ specification.

###  5.4 Alternative approach

Tchetgen and Shpitser (2012) and VanderWeele and Chiba (2014) considered
similar approaches for bounds on natural effects under the assumption that M-A
and Y-A ignorability holds but that Y-M ignorability does not. Both proposals
assume a known selection function that holds with equality rather than
inequality. We could modify our proposed approach in a similar spirit. For
example, we could assume that for all $(m_{1}^{\prime}\neq m_{1})$,
$(m_{2}^{\prime}\neq m_{2})$:

$\displaystyle\mu_{am_{1}m_{2}}^{\star}(m_{1}^{\prime},m_{2}^{\prime},x)-\mu_{a}(m_{1},m_{2},x)=f(\tau)[c\mu_{a}(m_{1},m_{2},x)+t]$

We could then recover $\psi_{M_{1}}(x)$ (and any averages of it) exactly as:

$\displaystyle\psi_{M_{1}}(x)$
$\displaystyle=\bar{\psi}_{M_{1}}(x)[1+cf(\tau)]$
$\displaystyle-c\sum_{m_{1},m_{2}}\mu_{a}(m_{1},m_{2},x)p(m_{1},m_{2}\mid
a,x)[p(m_{1}\mid a,x)-p(m_{1}\mid a^{\prime},x)]p(m_{2}\mid a^{\prime},x)$
$\displaystyle-tf(\tau)\sum_{m_{1},m_{2}}p(m_{1},m_{2}\mid a,x)[p(m_{1}\mid
a,x)-p(m_{1}\mid a^{\prime},x)]p(m_{2}\mid a^{\prime},x)$

While such an assumption would allow us to point identify $\psi_{M_{1}}(x)$,
the concept requires knowledge about a selection function that we are unlikely
to have. Despite the conservative inferences, we therefore prefer our proposed
approach.

###  5.5 Estimation

We extend all of the above methods to estimate the bounds on
$\psi_{M_{1}}(v)$. Theorem 3 in Appendix A provides the expressions for the
influence functions of the bounds on the average effect $\psi_{M_{1}}$.
Equipped with this expression, we can again use a projection estimator or the
DR-Learner to estimate the bounds.121212Moreover, the choice of $\tau$ need
only be valid across the domain of $x$ where the weights $w(v;V^{n})$ in the
second-stage regression are non-zero. Intuitively, these approaches share the
property that the upper bound on the convergence rates in the estimation is
governed by the products of errors in the nuisance estimation. Corollaries 3
and 4 in Appendix A give formal statements of these results using the lower
bound as an example, although we can derive an upper bound analogously. We
also provide expressions for the influence function for the bounds of
$\psi_{M_{2}}$, $\psi_{Cov}$, and $\psi_{IDE}$ in Appendix E and conjecture
that results analogous to Corollaries 3 and 4 can be derived for these
estimands. Finally, for a fixed $\tau$, Theorem 4 in Appendix A establishes
the conditions where the one-step estimator for the bounds on the average
effects is root-n consistent and asymptotically normal. We illustrate this
estimation procedure in the application in Section 6.

##  6 Application

To demonstrate these methods we revisit the data and application considered in
Rubinstein, Haviland and Breslau (2023), who examined the effect of COVID-19
vaccinations on depression, social isolation, and worries about health during
February 2021 using the COVID-19 Vaccine Trends and Impact Survey (CTIS). The
Delphi group at Carnegie Mellon University (CMU) conducted the CTIS from April
2020 through June 2022 in collaboration with the Facebook Data for Good group
(Salomon et al., 2021). Using this data, Rubinstein, Haviland and Breslau
(2023) posit a model that COVID-19 vaccinations affect depression via a direct
path, social isolation, and worries about health. Using the decomposition from
Vansteelandt and Daniel (2017), they found that pathways via social isolation
were more important than pathways via worries about health in explaining the
effect of COVID-19 vaccinations on depression. We refer to that paper for
details on the data and the limitations of this analysis. While this study
examined effect heterogeneity, the authors only examined heterogeneity within
discrete subgroups and primarily focused on the outcomes analysis. Moreover,
the authors did not find substantial effect heterogeneity with respect to the
mediation analysis among the specified subgroups.

We examine the decomposition of the total effect within the following subset
of CTIS respondents: employed, non-Hispanic, White respondents aged 25-54 with
at least a college degree, no chronic health conditions, who work outside the
home, and who had previously received an influenza vaccination. This included
a total of 13,764 individuals. Table 3 displays the average effect estimates
using influence-function based estimators, where the nuisance parameters were
estimated using twenty stacked XGBoost models with different hyperparameter
settings on the full dataset. While restricted to a much smaller subgroup,
these results are qualitatively comparable to the average estimates in
Rubinstein, Haviland and Breslau (2023).

Table 3: Average effect estimates on CTIS subset in February 2021, $N=13,764$ Estimand | Estimate | Lower 95% CI | Upper 95% CI
---|---|---|---
Total effect | -4.61 | -6.10 | -3.12
IIE - M1 | -1.86 | -2.45 | -1.28
IIE - M2 | -0.46 | -0.71 | -0.20
IIE - Cov | 0.08 | -0.15 | 0.30
IDE | -2.37 | -3.81 | -0.92

We next compare whether the interventional effects differed among those who
live in counties where Trump led Biden by 50 percentage points in the 2020
election (“Trump counties”), and those where Biden led Trump by 50 percentage
points (“Biden counties”). By limiting our sample to the subgroup defined
above, effect heterogeneity across the Biden vote share may proxy for how
social factors may moderate the mediated effects.131313Since we are unable to
fully control for socio-demographic variables, this variable may also pick up
on these moderating influence of these omitted factors that vary with the
Biden vote share. Specifically, we hypothesize that these relatively educated,
vaccine-accepting, and health conscious respondents who live in Trump-voting
counties may have lower total effects than those who live in Biden areas due
to the added stress of living in areas that generally took relatively fewer
COVID precautions. We similarly hypothesize that the effects via worries about
health might be lower in Trump-voting counties than Biden-voting counties for
this same reason. Figure 6 displays the results using both the DR-Learner and
projection estimators at these two points, where we use a simple linear model
for the projection.141414While we use sample-splitting to estimate $\eta$ and
construct influence function value estimates, we run the second-stage
regression on the entire sample instead of averaging two separate estimates.
Figure 11 in Appendix D display the entire estimated curves.

Figure 6: | Application results. Conditional total, indirect, and direct effects in Biden (Pct Dem Lead = 50) versus Trump (Pct Dem Lead = -50) counties

The total effect estimates are comparable in the Trump counties relative to
Biden counties; however, the projection estimates are slightly lower in Trump
relative to Biden counties while the DR-Learner suggests these effects may be
slightly higher.151515The uncertainty estimates for the proportion mediated
are obtained via the delta method, and the DR-Learner uncertainty estimates
are only valid assuming positive dependence between the errors in the models.
The uncertainty estimates for the DR-Learner also do not account for post-
selection inference and are therefore likely to be anti-conservative. On the
other hand, the effects via worries about health and social isolation are
nearly identical for both the the projection-estimator and DR-Learner. As seen
in Figure 11 in Appendix D, the chosen smoothing parameter ends up essentially
fitting a linear model for all functions other than the total effect.
Regardless, the point estimates are consistent with our expectations, where
effects via worries about health are lower in Trump counties relative to Biden
counties. Meanwhile, effects via isolation appear slightly larger in Trump
counties relative to Biden counties. However, all observed differences in
these effects are small relative to the uncertainty estimates and we are
unable to draw statistically significant conclusions.

###  6.1 Sensitivity analysis

We conduct a sensitivity analysis for $\psi_{M_{1}}$ both on average and as a
function of Biden’s vote share. Figure 7 displays the results assuming (28)
and where the conditional bounds are estimated using the DR-Learner.

We find that that our average effect estimates are robust to a $\tau$ as large
as 0.05, where $\tau$ parameterizes the deviations of the unobserved
counterfactual regression function to the observed regression function on the
risk-ratio scale (see equation 28). In other words, if this ratio were less
than $0.95$ ($1-0.05$), or greater than $1.05$ ($\frac{1}{1-0.05}$) for any
value of $(x,m_{1},m_{2})$, our bounds would include a null effect. Our
conditional effect estimates are less robust, in part due to the greater
uncertainty estimates. For example, our estimates for Trump counties is robust
only up to $\tau$ of 0.01 and for Biden counties is robust to $\tau$ of 0.03.

As a point of comparison, if we assumed that no unmeasured confounding held
conditional on $X$, but we failed to control for any covariates, across all
values of $(m_{1},m_{2})$ we would calculate a maximal $\tau=0.95$. To be
precise, we estimate that:

$\displaystyle\frac{1}{1-0.95}$
$\displaystyle\geq\frac{\mathbb{E}[Y^{m_{1}m_{2}}\mid A=a,M_{1}\neq
m_{1},M_{2}\neq m_{2}]}{\mathbb{E}[Y^{m_{1}m_{2}}\mid
A=a,M_{1}=m_{1},M_{2}=M_{2}]}$
$\displaystyle=\frac{\mathbb{E}[\mathbb{E}[Y\mid
A=a,M_{1}=m_{1},M_{2}=m_{2},X]]}{\mathbb{E}[Y\mid
A=a,M_{1}=m_{1},M_{2}=M_{2}]}\geq(1-0.95)$

where the equality holds via assuming no unmeasured confounding conditional on
$X$ and consistency. Therefore, a set of unmeasured confounders with
comparable association with the potential outcome regression would easily
explain away our estimated effects, as we find that we would be unable to rule
out a null effect at $\tau=0.05$. In other words, our significant effect would
disappear if there were some unmeasured confounder $U$ that were at least
approximately 5% (0.05 / 0.95) as associative with the outcome as our entire
observed covariate set. However, this comparison might be best thought of as a
“worst-case scenario,” as we estimate $\tau$ using all measured confounders
and our covariate set is quite rich. Interesting future work would be to
estimate different values of $\tau$ under different covariate subsets to
obtain possibly less conservative ranges of $\tau$. Sensitivity results for
the remaining parameters are available in Appendix D.

Figure 7: | Bounds for application. Bounds for average and conditional interventional indirect effects via social isolation as a function of $\tau$

##  7 Discussion

We propose two methods for estimating conditional average interventional
indirect effects: a semi-parametric projection-based approach and a fully non-
parametric approach. These procedures are conceptually simple: regress an
estimate of the uncentered influence function for the average parameter onto
the desired covariates. The projection-based estimator uses a parametric
regression model, and therefore targets a projection of the CIIE, while the
DR-Learner uses a fully non-parametric for this regression, and therefore
targets the CIIE itself. Our primary contribution is to establish the
conditions where the convergence rates of these estimators is equivalent to
that of an oracle regression of the true influence function onto these same
models. As with estimating the CATE, the error of these estimators is a
function of the product of errors in the nuisance estimation. However, unlike
the CATE, we must consider the sums of several products of nuisance functions,
which is in general a function of the cardinality of the joint mediators.
While our discussion focused primarily on estimating the effect via $M_{1}$,
this approach can be extended to estimate other interventional effects,
mediated effects, and likely a broad class of causal estimands.

As a second contribution we propose a sensitivity analysis for the conditional
effects that allows for mediator-outcome confounding. While the resulting
bounds may be quite wide in practice, they make only weak assumptions on the
underlying confounding mechanisms. Moreover, if one is willing to make
stronger assumptions on the selection mechanism, tighter bounds can be
obtained using a slight variant of our approach. We propose a general approach
to estimating these bounds using the projection estimators or DR-Learner,
where our results are again not tied to any particular estimation method. Our
methods also easily extend to estimating bounds on the average effects,
allowing for root-n consistent and asymptotically normal estimates under some
standard conditions.

Our proposed methods have several limitations: first, we only consider two
discrete mediators and a binary treatment. However, we could broaden this
general approach for more complex settings. For example, we could likely allow
for several mediators by regressing the corresponding influence functions
derived by Benkeser and Ran (2021) onto $V$. Similarly, we could likely extend
our results to allow for continuous mediators. This would require additional
assumptions, including, for example, the boundedness of the joint mediator
density. A complete treatment of this topic would be an interesting area for
future research. On the other hand, allowing for a continuous treatment would
be a more challenging problem as the causal estimands themselves would have to
be redefined, and an influence-function for the average effect does not exist.
A second limitation of our proposed method is that our sensitivity analysis
provides bounds that may be conservative. This is in part a function of the
fact that the methods we considered are all with respect to worst-case
scenarios that may occur infrequently in practice. A third limitation is that
we do not study any number of other possible non-parametric estimation
methods, such as an extension of the R-Learner proposed by Nie and Wager
(2021). Finally, we do not explore the minimax optimal rates for CIIE
estimation or propose estimators that might achieve these rates. Valuable
future work could explore any of these questions.

Acknowledgments: The authors would like to thank Amelia Haviland for helpful
discussions as this work developed. The authors would also like to thank the
two anonymous reviews and the Associate Editor for helpful comments,
questions, and suggestions that improved the quality of this manuscript.

Funding information: Authors state no funding involved.

Conflict of interest: Authors state no conflict of interest.

Data availability: The data used in the application for this study are not
publicly available, but are available on request at this link. Code is
available online at https://github.com/mrubinst757/ciie.

## References

  * Angrist and Pischke (2009) [author] Angrist, Joshua DJ. D. and Pischke, Jörn-SteffenJ.-S. (2009). Mostly harmless econometrics: An empiricist’s companion. Princeton university press.
  * Benkeser and Ran (2021) [author] Benkeser, DavidD. and Ran, JialuJ. (2021). Nonparametric inference for interventional effects with multiple mediators. Journal of Causal Inference 9 172–189.
  * Buja et al. (2019) [author] Buja, AndreasA., Brown, LawrenceL., Berk, RichardR., George, EdwardE., Pitkin, EmilE., Traskin, MikhailM., Zhang, KaiK. and Zhao, LindaL. (2019). Models as approximations I: Consequences illustrated with linear regression. Statistical Science 34 523–544.
  * Cuellar and Kennedy (2020) [author] Cuellar, MariaM. and Kennedy, Edward HE. H. (2020). A non-parametric projection-based estimator for the probability of causation, with application to water sanitation in Kenya. Journal of the Royal Statistical Society: Series A (Statistics in Society) 183 1793–1818.
  * Didelez, Dawid and Geneletti (2006) [author] Didelez, VV., Dawid, APA. and Geneletti, SS. (2006). Proceedings of the Twenty-Second Conference on Uncertainty in Artificial Intelligence.
  * Imai, Keele and Tingley (2010) [author] Imai, KosukeK., Keele, LukeL. and Tingley, DustinD. (2010). A general approach to causal mediation analysis. Psychological methods 15 309\.
  * Jackson (2020) [author] Jackson, John WJ. W. (2020). Meaningful causal decompositions in health equity research: definition, identification, and estimation through a weighting framework. Epidemiology 32 282–290.
  * Kennedy (2020) [author] Kennedy, Edward H.E. H. (2020). Towards optimal doubly robust estimation of heterogeneous causal effects. 10.48550/ARXIV.2004.14497
  * Kennedy (2022) [author] Kennedy, Edward HE. H. (2022). Semiparametric doubly robust targeted double machine learning: a review. arXiv preprint arXiv:2203.06469.
  * Kennedy, Balakrishnan and Wasserman (2021) [author] Kennedy, Edward HE. H., Balakrishnan, SivaramanS. and Wasserman, LarryL. (2021). Semiparametric counterfactual density estimation. arXiv preprint arXiv:2102.12034.
  * Kennedy, Lorch and Small (2019) [author] Kennedy, Edward HE. H., Lorch, ScottS. and Small, Dylan SD. S. (2019). Robust causal inference with continuous instruments using the local instrumental variable curve. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 81 121–143.
  * Lindmark, de Luna and Eriksson (2018) [author] Lindmark, AnitaA., de Luna, XavierX. and Eriksson, MarieM. (2018). Sensitivity analysis for unobserved confounding of direct and indirect effects using uncertainty intervals. Statistics in medicine 37 1744–1762.
  * Loh et al. (2020) [author] Loh, Wen WeiW. W., Moerkerke, BeatrijsB., Loeys, TomT. and Vansteelandt, StijnS. (2020). Heterogeneous indirect effects for multiple mediators using interventional effect models. Epidemiologic Methods 9\.
  * Luedtke, Diaz and van der Laan (2015) [author] Luedtke, Alexander RA. R., Diaz, IvanI. and van der Laan, Mark JM. J. (2015). The statistics of sensitivity analyses.
  * Miles (2022) [author] Miles, Caleb HC. H. (2022). On the Causal Interpretation of Randomized Interventional Indirect Effects. arXiv preprint arXiv:2203.00245.
  * Nie and Wager (2021) [author] Nie, XinkunX. and Wager, StefanS. (2021). Quasi-oracle estimation of heterogeneous treatment effects. Biometrika 108 299–319.
  * Park and Esterling (2021) [author] Park, SoojinS. and Esterling, Kevin MK. M. (2021). Sensitivity analysis for pretreatment confounding with multiple mediators. Journal of Educational and Behavioral Statistics 46 85–108.
  * Park, Qin and Lee (2020) [author] Park, SoojinS., Qin, XuX. and Lee, ChiounC. (2020). Estimation and sensitivity analysis for causal decomposition in health disparity research. Sociological Methods & Research 00491241211067516\.
  * Rubinstein, Haviland and Breslau (2023) [author] Rubinstein, MaxM., Haviland, AmeliaA. and Breslau, JoshuaJ. (2023). The effect of COVID-19 vaccinations on self-reported depression and anxiety during February 2021. Statistics and Public Policy 0 1-24. 10.1080/2330443X.2023.2190008
  * Salomon et al. (2021) [author] Salomon, Joshua AJ. A., Reinhart, AlexA., Bilinski, AlyssaA., Chua, Eu JingE. J., La Motte-Kerr, WichadaW., Rönn, Minttu MM. M., Reitsma, Marissa BM. B., Morris, Katherine AK. A., LaRocca, SarahS., Farag, Tamer HT. H. et al. (2021). The US COVID-19 Trends and Impact Survey: Continuous real-time measurement of COVID-19 symptoms, risks, protective behaviors, testing, and vaccination. Proceedings of the National Academy of Sciences 118 e2111454118.
  * Tchetgen and Shpitser (2012) [author] Tchetgen, Eric J TchetgenE. J. T. and Shpitser, IlyaI. (2012). Semiparametric theory for causal mediation analysis: efficiency bounds, multiple robustness, and sensitivity analysis. Annals of statistics 40 1816\.
  * Tsybakov (2004) [author] Tsybakov, Alexandre BA. B. (2004). Introduction to nonparametric estimation, 2009. URL https://doi. org/10.1007/b13794. Revised and extended from the 9\.
  * Van der Laan, Polley and Hubbard (2007) [author] Van der Laan, Mark JM. J., Polley, Eric CE. C. and Hubbard, Alan EA. E. (2007). Super learner. Statistical applications in genetics and molecular biology 6\.
  * VanderWeele and Chiba (2014) [author] VanderWeele, Tyler JT. J. and Chiba, YasutakaY. (2014). Sensitivity analysis for direct and indirect effects in the presence of exposure-induced mediator-outcome confounders. Epidemiology, biostatistics, and public health 11\.
  * Vansteelandt and Daniel (2017) [author] Vansteelandt, StijnS. and Daniel, Rhian MR. M. (2017). Interventional effects for mediation analysis with multiple mediators. Epidemiology (Cambridge, Mass.) 28 258\.

##  Appendix A Other definitions, theorems, and corollaries

###  A.1 Definitions and theorems from Kennedy (2020)

For completeness we reproduce key results from Kennedy (2020) that we use in
our proofs. In particular, we reproduce the definition of estimator stability,
a theorem that linear smoothers are stable (Theorem 1 from Kennedy (2020)),
and a proposition about the asymptotics of pseudo-outcome regression where the
second-stage regression is a linear smoother (Proposition 1 from Kennedy
(2020)). The proofs can be found in Kennedy (2020).

######  Definition 1 (Stability).

Suppose $D^{n}=(Z_{01},...,Z_{0n})$ and $Z^{n}=(Z_{1},...,Z_{n})$ are
independent test and training samples, respectively, with covariates
$X_{i}\subset Z_{i}$ and $V_{i}\subseteq X_{i}$. Let

  1. 1.

$\hat{f}(z)=\hat{f}(z;D^{n})$ be an estimate of the function $f$ using
training data $D^{n}$

  2. 2.

$\hat{b}(x)=\mathbb{E}\\{\hat{f}(Z)-f(Z)\mid X=x,D^{n}\\}$ be the conditional
bias of the estimator $\hat{f}$

  3. 3.

$\hat{\mathbb{E}}_{n}(Y\mid V=v)$ denote a generic regression estimator that
regresses outcomes $(Y_{1},...,Y_{n})$ on the covariates $(V_{1},...,V_{n})$
in the test sample $Z^{n}$.

Then the regression is defined as “stable” with respect to distance metric $d$
if

$\displaystyle\frac{\hat{\mathbb{E}}_{n}\\{\hat{f}(Z)\mid
V=v\\}-\hat{\mathbb{E}}_{n}\\{f(Z)\mid
V=v\\}-\hat{\mathbb{E}}_{n}\\{\hat{b}(X)\mid
V=v\\}}{\sqrt{\mathbb{E}\left(\left[\hat{\mathbb{E}}_{n}\\{f(Z)\mid
V=v\\}-\mathbb{E}\\{f(Z)\mid V=v\\}\right]^{2}\right)}}\to^{p}0$ (31)

whenever $d(\hat{f},f)\to^{p}0$.

######  Theorem 2 (Theorem 1 from Kennedy (2020)).

Linear smoothers of the form $\hat{\mathbb{E}}_{n}(\hat{f}(Z)\mid
V=v)=\sum_{i=1}^{n}w(v;V^{n})\hat{f}(Z_{i})$ are stable with respect to the
distance

$\displaystyle
d(\hat{f},f)=\|\hat{f}-f\|_{w^{2}}=\sum_{i=1}^{n}\frac{w_{i}(v;V^{n})^{2}}{\sum_{j=1}^{n}w_{j}(v;V^{n})^{2}}\int\left\\{\hat{f}(z)-f(z)\right\\}^{2}dP(z\mid
V_{i})$ (32)

whenever $1/\|\sigma\|_{w^{2}}=\mathcal{O}_{p}(1)$ for
$\sigma^{2}(v)=Var(f(Z)\mid V=v)$.

######  Proposition 4 (Proposition 2 from Kennedy (2020)).

If $\hat{b}(x)=\hat{b}_{1}(x)\hat{b}_{2}(x)$, and $\hat{\mathbb{E}}_{n}$ is a
linear smoother with $\sum_{i}\lvert
w_{i}(v;V^{n})\rvert=\mathcal{O}_{p}(c_{n})$, then

$\displaystyle\hat{\mathbb{E}}_{n}\\{\hat{b}(X)\mid
V=v\\}=\mathcal{O}_{p}(c_{n}\|\hat{b}_{1}\|_{w,p}\|\hat{b}_{2}\|_{w,q})$ (33)

for the norm $\|f\|_{w,p}=\left[\sum_{i}\left\\{\frac{\lvert
w_{i}(v;V^{n})\rvert}{\sum_{j}\lvert w_{j}(v;V^{n})\rvert}\right\\}\lvert
f(X_{i})^{p}\rvert\right]^{1/p}$ for $1/p+1/q=1$ and $p,q>1$.

###  A.2 Other theorems and corollaries

This section provides a statement of Theorem 3, which provides an expression
for the efficient influence function of the bounds on $\psi_{M_{1}}$. We then
state the implied one-step estimator for the bounds on the average effect, and
establish the conditions for the root-n consistency and asymptotic normality
of these estimates in Theorem 4. We then extend the results of Theorem 1
applied to the projection estimator of the bounds in Corollary 3, and the
results of Corollary 1 for the DR-Learner in Corollary 4. Proofs are included
in Section B.

######  Theorem 3 (Efficient influence function of bounds).

Define

$\displaystyle\Gamma_{1,a}$
$\displaystyle=\mathbb{E}\left\\{\sum_{m_{1},m_{2}}\mu_{a}(m_{1},m_{2},X)p(m_{1},m_{2}\mid
a,X)p(m_{1}\mid a,X)p(m_{2}\mid
a^{\prime},X)\right\\}=\mathbb{E}\left\\{\zeta_{1,a}(X)]\right\\}$
$\displaystyle\Gamma_{1,a^{\prime}}$
$\displaystyle=\mathbb{E}\left\\{\sum_{m_{1},m_{2}}\mu_{a}(m_{1},m_{2},X)p(m_{1},m_{2}\mid
a,X)p(m_{1}\mid a^{\prime},X)p(m_{2}\mid
a^{\prime},X)\right\\}=\mathbb{E}\left\\{\zeta_{1,a^{\prime}}(X)]\right\\}$
$\displaystyle\Gamma_{2,a}$
$\displaystyle=\mathbb{E}\left\\{\sum_{m_{2},m_{2}}p(m_{1},m_{2}\mid
a,X)p(m_{1}\mid a,X)p(m_{2}\mid
a^{\prime},X)\right\\}=\mathbb{E}\left\\{\zeta_{2,a}(X)]\right\\}$
$\displaystyle\Gamma_{2,a^{\prime}}$
$\displaystyle=\mathbb{E}\left\\{\sum_{m_{2},m_{2}}p(m_{1},m_{2}\mid
a,X)p(m_{1}\mid a^{\prime},X)p(m_{2}\mid
a^{\prime},X)\right\\}=\mathbb{E}\left\\{\zeta_{2,a^{\prime}}(X)]\right\\}$

Let $\phi_{1,a}(Z)$, $\phi_{1,a^{\prime}}(Z)$, $\phi_{2,a}(Z)$, and
$\phi_{2,a^{\prime}}(Z)$ be the uncentered efficient influence functions for
$\Gamma_{1,a}$, $\Gamma_{1,a^{\prime}}$, $\Gamma_{2,a}$, and
$\Gamma_{2,a^{\prime}}$, respectively. The efficient influence function for
each quantity is given below:

$\displaystyle\phi_{1,a}(Z)$
$\displaystyle=\frac{\mathds{1}(A=a)}{\pi_{a}(X)}[Yp(M_{1}\mid a,X)p(M_{2}\mid
a^{\prime},X)-\zeta_{1,a}(X)]$
$\displaystyle+\frac{\mathds{1}(A=a)}{\pi_{a}(X)}\left[\sum_{m_{2}}\mu_{a}(M_{1},m_{2},X)p(M_{1},m_{2}\mid,a,X)p(m_{2}\mid
a^{\prime},X)-\zeta_{1,a}(X)\right]$
$\displaystyle+\frac{\mathds{1}(A=a^{\prime})}{\pi_{a^{\prime}}(X)}\left[\sum_{m_{1}}\mu_{a}(m_{1},M_{2},X)p(m_{1},M_{2}\mid
a,X)p(m_{1}\mid a,X)-\zeta_{1,a}(X)\right]+\zeta_{1,a}(X)$
$\displaystyle\phi_{1,a^{\prime}}(Z)$
$\displaystyle=\frac{\mathds{1}(A=a)}{\pi_{a}(X)}[Yp(M_{1}\mid
a^{\prime},X)p(M_{2}\mid a^{\prime},X)-\zeta_{1,a^{\prime}}(X)]+$
$\displaystyle+\frac{\mathds{1}(A=a^{\prime})}{\pi_{a^{\prime}}(X)}\left[\sum_{m_{2}}\mu_{a}(M_{1},m_{2},X)p(M_{1},m_{2}\mid
a,X)p(m_{2}\mid a^{\prime},X)-\zeta_{1,a^{\prime}}(X)\right]$
$\displaystyle+\frac{\mathds{1}(A=a^{\prime})}{\pi_{a^{\prime}}(X)}\left[\sum_{m_{1}}\mu_{a}(m_{1},M_{2},X)p(m_{1},M_{2}\mid
a,X)p(m_{1}\mid
a^{\prime},X)-\zeta_{1,a^{\prime}}(X)\right]+\zeta_{1,a^{\prime}}(X)$
$\displaystyle\phi_{2,a}(Z)$
$\displaystyle=\frac{\mathds{1}(A=a)}{\pi_{a}(X)}[p(M_{1}\mid a,X)p(M_{2}\mid
a^{\prime},X)-\zeta_{2,a}(X)]$
$\displaystyle+\frac{\mathds{1}(A=a)}{\pi_{a}(X)}\left[\sum_{m_{2}}p(m_{2}\mid
a^{\prime},X)p(M_{1},m_{2}\mid,a,x)-\zeta_{2,a}(X)\right]$
$\displaystyle+\frac{\mathds{1}(A=a^{\prime})}{\pi_{a^{\prime}}(X)}\left[\sum_{m_{1}}p(m_{1}\mid
a,X)p(m_{1},M_{2}\mid a,x)-\zeta_{2,a}(X)\right]+\zeta_{2,a}(X)$
$\displaystyle\phi_{2,a^{\prime}}(Z)$
$\displaystyle=\frac{\mathds{1}(A=a)}{\pi_{a}(X)}\left[p(M_{1}\mid
a^{\prime},X)p(M_{2}\mid a^{\prime},X)-\zeta_{2,a^{\prime}}(X)\right]$
$\displaystyle+\frac{\mathds{1}(A=a^{\prime})}{\pi_{a^{\prime}}(X)}\left[\sum_{m_{2}}p(m_{2}\mid
a^{\prime},X)p(M_{1},m_{2}\mid,a,x)-\zeta_{2,a^{\prime}}(X)\right]$
$\displaystyle+\frac{\mathds{1}(A=a^{\prime})}{\pi_{a^{\prime}}(X)}\left[\sum_{m_{1}}p(m_{1}\mid
a^{\prime},X)p(m_{1},M_{2}\mid
a,x)-\zeta_{2,a^{\prime}}(X)\right]+\zeta_{2,a^{\prime}}(X)$

Recalling the expressions for the bounds in equations (29)-(30), the
(uncentered) efficient influence function for the parameters
$\psi_{M_{1},lb}(\tau)$ and $\psi_{M_{1},ub}(\tau)$ are therefore equal to

$\displaystyle\Xi_{ub}(Z;\eta,\tau)$
$\displaystyle=\varphi(Z;\eta)+\varphi_{a}(Z;\eta)f_{u}(\tau)c_{u}-\varphi_{a^{\prime}}(Z)f_{l}(\tau)c_{l}+t_{u}f_{u}(\tau)-t_{l}f_{l}(\tau)$
$\displaystyle-
f_{u}(\tau)[c_{u}\phi_{1,a}+t_{u}\phi_{1,a^{\prime}}]+f_{l}(\tau)[c_{l}\phi_{1,a^{\prime}}+t_{l}\phi_{2,a^{\prime}}]$
$\displaystyle\Xi_{lb}(Z;\eta,\tau)$
$\displaystyle=\varphi(Z;\eta)+\varphi_{a}(Z;\eta)f_{l}(\tau)c_{l}-\varphi_{a^{\prime}}(Z)f_{u}(\tau)c_{u}+t_{l}f_{l}(\tau)-t_{u}f_{u}(\tau)$
$\displaystyle-
f_{l}(\tau)[c_{l}\phi_{1,a}+t_{l}\phi_{1,a^{\prime}}]+f_{u}(\tau)[c_{u}\phi_{1,a^{\prime}}+\phi_{2,a^{\prime}}]$

######  Theorem 4 (Asymptotic normality of estimates of bounds).

Given $\tau$, consider the one-step estimator of the bounds:

$\displaystyle[\hat{\psi}_{M_{1},lb}(\tau),\hat{\psi}_{M_{1},ub}(\tau)]^{\top}$
$\displaystyle=\left[\mathbb{P}_{n}[\Xi_{lb}(Z;\hat{\eta},\tau)],\mathbb{P}_{n}[\Xi_{ub}(Z;\hat{\eta},\tau)]\right]^{\top}]]$
$\displaystyle=\mathbb{P}_{n}\Xi(Z;\hat{\eta},\tau)=\hat{\psi}_{M_{1},b}$ (34)

Let $\hat{\psi}_{M_{1},b}$ denote the vector of estimates in (34). Assume that

  1. 1.

Nuisance parameters $\eta$ are estimated using sample-splitting

  2. 2.

$\|\hat{\Xi}-\Xi\|^{2}=o_{p}(1)$

  3. 3.

$\|\hat{\eta}-\eta\|=o_{p}(n^{-1/4})$

We then obtain the following limiting distribution for the vector of
estimates:

$\displaystyle\sqrt{n}\mathbb{P}_{n}[\hat{\psi}_{M_{1},b}(\tau)-\psi_{M_{1},b}(\tau)]\to^{d}\mathcal{N}(0,\mathbb{E}[\Xi\Xi^{\top}])$
(35)

######  Corollary 3 (Projection estimator of bounds).

For fixed parameter $\tau$, consider the function $f(V;\alpha)$, the moment
condition,

$\displaystyle\mathbb{E}\left[w(X)\frac{\partial
f(V;\alpha)}{\partial\alpha}[\psi_{M_{1},lb}(Z;\tau)-f(V;\alpha)]\right]=0$

the associated influence curve at any $\alpha^{\star}$,

$\displaystyle\phi_{lb}(Z;\alpha^{\star},\eta,\tau)$
$\displaystyle=\frac{\partial
f(V;\alpha)}{\partial\alpha}w(X)\\{\Xi_{lb}(Z;\eta,\tau)-f(V;\alpha^{\star})\\}$

and the associated influence-function based estimating equation for
$\hat{\alpha}$:

$\displaystyle\mathbb{P}_{n}\left[\frac{\partial
f(V;\hat{\alpha})}{\partial\alpha}w(X)\\{\hat{\Xi}_{lb}(Z;\hat{\eta},\tau)-f(V;\hat{\alpha})\\}\right]=0$

Consider the moment condition
$\mathbb{E}[\phi_{lb}(Z;\alpha_{0},\eta_{0},\tau)]=0$ evaluated at the true
parameters $(\alpha_{0},\eta_{0})$. Now consider the estimator $\hat{\alpha}$
that satisfies $\mathbb{P}_{n}[\phi_{lb}(Z;\hat{\alpha},\hat{\eta},\tau)]=0$,
where $\hat{\eta}$ is estimated on an independent sample. Assume that:

  * •

The function class $\\{\phi_{lb}(Z;\alpha,\eta):\alpha\in\mathbb{R}^{p}\\}$ is
Donsker in $\alpha$ for any fixed $\eta$

  * •

$\|\phi_{lb}(Z;\hat{\alpha},\hat{\eta},\tau)-\phi(Z;\alpha_{0},\eta_{0},\tau)\|=o_{p}(1)$

  * •

The map $\alpha\to\mathbb{P}[\phi_{lb}(Z;\alpha,\eta,\tau)]$ is differentiable
at $\alpha_{0}$ uniformly in the true $\eta$, with non-singular derivative
matrix
$\frac{\partial}{\partial\alpha}\mathbb{P}\\{\phi_{lb}(Z;\alpha,\eta,\tau\\}\mid_{\alpha=\alpha_{0}}=G(\alpha_{0},\eta,\tau)$,
where $G(\alpha_{0},\hat{\eta},\tau)\to^{p}M(\alpha_{0},\eta_{0})$

Then

$\displaystyle\hat{\beta}-\beta$
$\displaystyle=-G^{-1}[\mathbb{P}_{n}-P]\phi_{lb}(Z;\alpha_{0},\eta_{0},\tau)+\mathcal{O}_{p}(T_{1n}+T_{2n}+T_{3n}+T_{4n}$
$\displaystyle+\underbrace{\|p(M_{1},M_{2}\mid a,X)-\hat{p}(M_{1},M_{2}\mid
a,X)\|\|p(M_{1}\mid a^{\prime},X)-\hat{p}(M_{1}\mid a^{\prime},X)\|}_{(xi)}$
$\displaystyle+\underbrace{\|p(M_{1},M_{2}\mid a,X)-\hat{p}(M_{1},M_{2}\mid
a,X)\|\|p(M_{2}\mid a^{\prime},X)-\hat{p}(M_{2}\mid a^{\prime},X)\|}_{(xii)}$
$\displaystyle+\underbrace{\|p(M_{1},M_{2}\mid a,X)-\hat{p}(M_{1},M_{2}\mid
a,X)\|\|p(M_{1}\mid a,X)-\hat{p}(M_{1}\mid a,X)\|}_{(xiii)}$
$\displaystyle+\underbrace{\|p(M_{1},M_{2}\mid a,X)-\hat{p}(M_{1},M_{2}\mid
a,X)\|\|\pi_{a}(X)-\hat{\pi}_{a}(X)\|}_{(xiv)})+o_{p}(1/\sqrt{n})$

where the terms $T_{1n}$-$T_{4n}$ are defined in Theorem 1. Moreover, if
$Z_{n}=T_{1n}+T_{2n}+T_{3n}+T_{4n}+(x)-(iv)$ are $o_{p}(n^{1/2})$, then

$\displaystyle\sqrt{n}(\hat{\alpha}-\alpha)\to^{d}\mathcal{N}(0,G^{-1}\mathbb{E}[\phi_{lb}\phi_{lb}^{\top}]G^{-\top})$

where
$G=\frac{\partial}{\partial\alpha}\mathbb{P}\\{\phi_{lb}(Z;\alpha,\eta,\tau)\\}_{\alpha=\alpha_{0}}$.

######  Remark 5.

There are 21 possible combinations of second-order terms among nuisance
elements $[\pi_{a}(X),\mu_{a}(M_{1},M_{2},X),p(M_{1}\mid a,X),p(M_{1}\mid
a^{\prime},X),p(M_{2}\mid a^{\prime},X),p(M_{1},M_{2}\mid a,X)]$.
$\hat{\psi}^{dr}_{M_{1}}(x)$ requires considering ten of these error products.
Estimating $\hat{\psi}^{dr}_{M_{1},lb}(x)$ requires considering the four
additional terms noted above.

On the other hand, if the estimate of $p(M_{1}\mid a,X)$ comes from an
estimate of $p(M_{1},M_{2}\mid a,X)$ (and similarly $p(M_{1}\mid
a^{\prime},x)$ and $p(M_{2}\mid a^{\prime},x)$ come from an estimate of
$p(M_{1},M_{2}\mid a^{\prime},X)$), then these expressions simplify. In that
case there are only 10 total possible second-order terms. $\psi_{M_{1}}$
requires estimating seven of the ten and the bounds only require one
additional error product that comes from (xiii), which is on the same order of
the square of the L2-norm of the error of the estimated joint mediator
probability $p(M_{1},M_{2}\mid a,X)$.

######  Remark 6.

This result is agnostic as to which of equations (26)-(28) are assumed for the
sensitivity analysis. While not all second-order remainder terms from
$\phi_{1,a}$, $\phi_{1,a^{\prime}}$, $\phi_{2,a}$, $\phi_{2,a^{\prime}}$ will
appear in each expression, as the constants on these functions may be zero
depending on the assumption used, these same terms will appear regardless.
This can be seen from the Proof of Theorem 4, where we derive the second-order
terms for

$\displaystyle
P\left[[\hat{\phi}_{1,a},\hat{\phi}_{1,a^{\prime}},\hat{\phi}_{2,a},\hat{\phi}_{2,a^{\prime}}]^{\top}\right]-P\left[[\phi_{1,a},\phi_{1,a^{\prime}},\phi_{2,a},\phi_{2,a^{\prime}}]^{\top}\right]$

######  Remark 7.

It is easy to show an analogous result for a projection estimator of
$\psi_{M_{1},ub}(v)$.

######  Corollary 4 (DR-Learner for bounds).

Define Algorithm 2 as Algorithm 1, where Step 2 instead constructs the
estimated pseudo-outcomes $\Xi_{lb}(Z;\hat{\eta},\tau)$. Assume that the
conditions of Corollary 1 hold and that $\tau\neq 0$.

Denote $\tilde{\psi}_{M_{1},lb}(v;\tau)$ as an oracle estimator from a
regression of $\Xi_{lb}(Z;\eta,\tau)$ onto $V$ and let $Q_{n}^{2\star}(v)$
denote the oracle risk at the point $V=v$. Then

$\displaystyle[\hat{\psi}_{M_{1},lb}(\tau)-\psi_{M_{1},lb}(\tau)]\lesssim\hat{\mathbb{E}}_{n}(\hat{q}_{lb}(X)\mid
V=v)+o_{p}(Q^{\star}_{n}(v))$

where

$\displaystyle\hat{q}_{lb}(x)$
$\displaystyle=\hat{b}(x)+(\pi_{a}(x)-\hat{\pi}_{a}(x))\sum_{m_{1},m_{2}}(p(m_{1},m_{2}\mid
a,x)-\hat{p}(m_{1},m_{2}\mid a,x))$
$\displaystyle+\sum_{m_{1},m_{2}}(p(m_{1},m_{2}\mid
a,x)-\hat{p}(m_{1},m_{2}\mid a,x))(p(m_{1}\mid a,x)-\hat{p}(m_{1}\mid a,x))$
$\displaystyle+\sum_{m_{1},m_{2}}(p(m_{1},m_{2}\mid
a,x)-\hat{p}(m_{1},m_{2}\mid a,x))(p(m_{1}\mid a^{\prime},x)-\hat{p}(m_{1}\mid
a^{\prime},x))$ $\displaystyle+\sum_{m_{1},m_{2}}(p(m_{1},m_{2}\mid
a,x)-\hat{p}(m_{1},m_{2}\mid a,x))(p(m_{2}\mid a^{\prime},x)-\hat{p}(m_{2}\mid
a^{\prime},x))$

Moreover, if $\hat{\mathbb{E}}_{n}$ is a linear smoother with weights
$\sum_{i}\lvert w_{i}(v;V^{n})\rvert=\mathcal{O}_{p}(a_{n})$, then

$\displaystyle\hat{\mathbb{E}}_{n}(\hat{q}_{lb}(X)\mid
V=v)=\max_{j}\mathcal{O}_{p}\left(a_{n}\|\hat{q}_{lb,j1}\|_{w,2}\|\hat{q}_{lb,j2}\|_{w,2}\right)$

where $\hat{q}_{lb,j}$ represents the $j$-th second-order term in the
expression of $\hat{q}_{lb}(x)$ above.

######  Remark 8.

The comments in Remarks 5 and 6 apply here as well. The result for an estimate
of $\psi_{M_{1},ub}(v)$ is asymptotically equivalent.

##  Appendix B Proofs

This section contains proofs of all theorems, corollaries, and propositions in
the main paper and in Section A.2. We divide the proofs into two sections: the
first pertaining to conditional effect estimation; the second pertaining to
the bounds on $\psi_{M_{1}}$, both average and conditional. We defer tedious
algebraic derivations to Section F.

####  B.0.1 Conditional effect estimation

The following proofs pertain to estimating the conditional effects. This
section contains four proofs:

  * •

Proposition 1: efficient influence function for the projection estimator

  * •

Theorem 1: root-n consistency and asymptotic normality of projection estimator

  * •

Corollaries 1 and 2: DR-Learner rates of convergence

###### Proof of Proposition 1.

We consider the functional $\Psi(\beta;\mathbb{P})$ defined in
(LABEL:eqn:projestimator) and derive its efficient influence curve. This
result follows directly from the result previously derived in Cuellar and
Kennedy (2020), though we reproduce it here for completeness. We treat the
data as discrete throughout to simplify the derivation. Let
$h(x)=\frac{\partial g(v;\beta)}{\partial\beta}w(x)$. We first define our
parameter:

$\displaystyle\Psi(\beta;P)$
$\displaystyle=\mathbb{E}[h(X)[\psi_{M_{1}}(X)-g(V;\beta)]]$
$\displaystyle=\sum_{x}h(x)[\psi_{M_{1}}(x)-g(v;\beta)]p(x)$

We can use the product rule to obtain that:

$\displaystyle\text{IF}[\Psi(\beta)]$
$\displaystyle=\sum_{x}h(x)\text{IF}[\psi_{M_{1}}(x)]p(x)+h(x)[\psi_{M_{1}}(x)-g(v;\beta)]\text{IF}[p(x)]$
$\displaystyle=\sum_{x}h(x)\text{IF}[\psi_{M_{1}}(x)]p(x)+h(x)[\psi_{M_{1}}(x)-g(v;\beta)][\mathds{1}(X=x)-p(x)]$
$\displaystyle=h(X)[\varphi(Z)-g(V;\beta)]-\Psi(\beta;P)$
$\displaystyle=\phi(Z)-\Psi(\beta;P)$

where $\varphi(Z)$ takes the form defined in (11). ∎

###### Proof of Theorem 1.

This result follows directly from Lemma 3 from Kennedy, Balakrishnan and
Wasserman (2021). We further make the assumptions outlined in Theorem 1, which
correspond to the same assumptions of Lemma 3. The following term drives our
ability to obtain root-n consistency and asymptotically normal estimates:

$\displaystyle\mathbb{P}[\phi(Z;\beta,\hat{\eta})-\phi(Z;\beta,\eta)]$

This is simply equal to:

$\displaystyle\mathbb{P}[\varphi(Z;\hat{\eta})-\varphi(Z;\eta)]$

We have shown that this expression is bounded by the product of errors of the
nuisance functions (these products are reproduced and used to prove Corollary
1). These terms then characterize the rate $R_{n}=T_{1n}+T_{2n}+T_{3n}+T_{4n}$
in the $\mathcal{O}_{p}(R_{n})$ term in Lemma 3. An immediate consequence is
that if $R_{n}=o_{p}(n^{-1/4})$ then this term is of smaller order and thus

$\displaystyle\hat{\beta}-\beta=-M^{-1}[\mathbb{P}_{n}-\mathbb{P}]\phi(Z;\beta,\eta)+o_{p}(n^{-1/2})$

Because this expression is simply a centered sample mean of a fixed function,
the asymptotic normality of the difference scaled by $\sqrt{n}$ follows via
the Central Limit Theorem with variance equal to
$M^{-1}\mathbb{E}[\phi\phi^{\top}]M^{-\top}$. ∎

###### Proof of Corollary 1.

Let $\varphi(Z;\eta)$ be the (uncentered) efficient influence function for
$\psi_{M_{1}}$. Let
$\hat{\psi}_{M_{1}}=n^{-1}\sum_{i=1}^{n}\varphi(Z;\hat{\eta})$. Recall that
$P\\{f(Z)\\}=\int f(z)dP(z\mid D_{1}^{n})$ where $D_{1}^{n}$ is an independent
sample of size $n$. To ease notation we omit the conditioning on the training
samples for the remainder of this proof; however, the consequence is that we
derive the results thinking of the estimated functions as fixed and not
dependent on the data. We show in Section F that
$P[\hat{\psi}_{M_{1}}-\psi_{M_{1}}]$ decomposes into the following
form:161616Benkeser and Ran (2021) show a similar derivation; our
decomposition largely agrees though we find a few additional second-order
terms. They do not provide a proof of their result in their paper, and we
therefore cannot compare our derivations; however, we include ours here for
completeness. The differences between our decomposition and theirs do not
substantively affect our results.

$\displaystyle P[\hat{\psi}_{M_{1}}-\psi_{M_{1}}]$
$\displaystyle=\mathbb{E}\left[\frac{\pi_{a}(X)}{\hat{\pi}_{a}(X)}\sum_{m_{1},m_{2}}\frac{(\hat{p}(m_{1}\mid
a,X)-\hat{p}(m_{1}\mid a^{\prime},X))\hat{p}(m_{2}\mid
a^{\prime},X)}{\hat{p}(m_{1},m_{2}\mid a^{\prime},X)}\right.$ (36)
$\displaystyle\left.\times\left(p(m_{1},m_{2}\mid a,X)-\hat{p}(m_{1},m_{2}\mid
a,X))(\mu_{a}(m_{1},m_{2},X)-\hat{\mu}_{a}(m_{1},m_{2},X))\right)\right]$
$\displaystyle+\mathbb{E}\left[\left(\frac{\pi_{a}(X)-\hat{\pi}_{a}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}\hat{p}(m_{2}\mid
a^{\prime},X)[\hat{p}(m_{1}\mid a,X)-\hat{p}(m_{1}\mid a^{\prime},X)]\right.$
(37)
$\displaystyle\left.\times(\mu_{a}(m_{1},m_{2},X)-\hat{\mu}_{a}(m_{1},m_{2},X))\right]$
$\displaystyle+\mathbb{E}\left[\left(\frac{\pi_{a}(X)-\hat{\pi}_{a}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},X)\hat{p}(m_{2}\mid
a^{\prime},X)[(p(m_{1}\mid a,X)-\hat{p}(m_{1}\mid a,X))\right.$ (38)
$\displaystyle\left.+(p(m_{1}\mid a^{\prime},X)-\hat{p}(m_{1}\mid
a^{\prime},X))\right]$
$\displaystyle+\mathbb{E}\left[\left(\frac{\pi_{a^{\prime}}(X)-\pi_{a^{\prime}}(X)}{\hat{\pi}_{a^{\prime}}(X)}\right)\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},X)\right.$
(39) $\displaystyle\left.\times[\hat{p}(m_{1}\mid a,X)-\hat{p}(m_{1}\mid
a^{\prime},X)][p(m_{2}\mid a,X)-\hat{p}(m_{2}\mid a^{\prime},X)\right]$
$\displaystyle-\mathbb{E}\left[\sum_{m_{1},m_{2}}\hat{p}(m_{2}\mid
a^{\prime},X)(\mu_{a}(m_{1},m_{2},X)-\hat{\mu}_{a}(m_{1},m_{2},X))\right.$
(40) $\displaystyle\left.\times((p(m_{1}\mid a,X)-\hat{p}(m_{1}\mid
a,X))-(p(m_{1}\mid a^{\prime},X)-\hat{p}(m_{1}\mid a^{\prime},X))\right]$
$\displaystyle-\mathbb{E}\left[\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},X)[p(m_{2}\mid
a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X)]\right.$ (41)
$\displaystyle\left.\times((p(m_{1}\mid a,X)-\hat{p}(m_{1}\mid
a,X))-(p(m_{1}\mid a^{\prime},X)-\hat{p}(m_{1}\mid a^{\prime},X))\right]$
$\displaystyle-\mathbb{E}\left[\sum_{m_{1},m_{2}}(p(m_{1}\mid a,X)-p(m_{1}\mid
a^{\prime},X))\right.$ (42)
$\displaystyle\left.\times(\mu_{a}(m_{1},m_{2},X)-\hat{\mu}_{a}(m_{1},m_{2},X))(p(m_{2}\mid
a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X)\right]$

We use this result to bound the error of the DR-Learner at a given point
$V=v$. We first define the conditional bias at $X=x$ as
$\hat{b}^{\star}(x)=\mathbb{E}\\{\hat{\varphi}(Z)-\varphi(Z)\mid
X=x,D_{1}^{n}\\}$. Following similar logic we obtain the following expression
for $\hat{b}^{\star}(x)$:

$\displaystyle\hat{b}^{\star}(x)=\frac{\pi_{a}(x)}{\hat{\pi}_{a}(x)}\left[\sum_{m_{1},m_{2}}\frac{(\hat{p}(m_{1}\mid
a,x)-\hat{p}(m_{1}\mid a^{\prime},x))\hat{p}(m_{2}\mid
a^{\prime},x)}{\hat{p}(m_{1},m_{2}\mid a^{\prime},x)}\right.$
$\displaystyle\left.\times\left(p(m_{1},m_{2}\mid a,x)-\hat{p}(m_{1},m_{2}\mid
a,x))(\mu_{a}(m_{1},m_{2},x)-\hat{\mu}_{a}(m_{1},m_{2},x))\right)\right]$
$\displaystyle+\left(\frac{\pi_{a}(x)-\hat{\pi}_{a}(x)}{\hat{\pi}_{a}(x)}\right)\left[\sum_{m_{1},m_{2}}\hat{p}(m_{2}\mid
a^{\prime},x)[\hat{p}(m_{1}\mid a,x)-\hat{p}(m_{1}\mid a^{\prime},x)]\right.$
$\displaystyle\left.\times(\mu_{a}(m_{1},m_{2},x)-\hat{\mu}_{a}(m_{1},m_{2},x))\right]$
$\displaystyle+\left[\left(\frac{\pi_{a}(x)-\hat{\pi}_{a}(x)}{\hat{\pi}_{a}(x)}\right)\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},x)\hat{p}(m_{2}\mid
a^{\prime},x)[p(m_{1}\mid a,x)-\hat{p}(m_{1}\mid a,x)\right.$
$\displaystyle\left.+(p(m_{1}\mid a^{\prime},x)-\hat{p}(m_{1}\mid
a^{\prime},x))\right]$
$\displaystyle+\left(\frac{\pi_{a^{\prime}}(x)-\hat{\pi}_{a^{\prime}}(x)}{\hat{\pi}_{a^{\prime}}(x)}\right)\left[\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},x)\right.$
$\displaystyle\left.\times[\hat{p}(m_{1}\mid a,x)-\hat{p}(m_{1}\mid
a^{\prime},x)][p(m_{2}\mid a,x)-\hat{p}(m_{2}\mid a^{\prime},x)\right]$
$\displaystyle-\left[\sum_{m_{1},m_{2}}\hat{p}(m_{2}\mid
a^{\prime},x)(\mu_{a}(m_{1},m_{2},x)-\hat{\mu}_{a}(m_{1},m_{2},x))\right.$
$\displaystyle\left.\times(p(m_{1}\mid a,x)-\hat{p}(m_{1}\mid
a,x)-(p(m_{1}\mid a^{\prime},x)-\hat{p}(m_{1}\mid a^{\prime},x))\right]$
$\displaystyle-\left[\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},x)[p(m_{2}\mid
a^{\prime},x)-\hat{p}(m_{2}\mid a^{\prime},x)]\right.$
$\displaystyle\left.\times(p(m_{1}\mid a,x)-\hat{p}(m_{1}\mid
a,x)-(p(m_{1}\mid a^{\prime},x)-\hat{p}(m_{1}\mid a^{\prime},x))\right]$
$\displaystyle-\left[\sum_{m_{1},m_{2}}(p(m_{1}\mid a,x)-p(m_{1}\mid
a^{\prime},x))\times(\mu_{a}(m_{1},m_{2},x)-\hat{\mu}_{a}(m_{1},m_{2},x))(p(m_{2}\mid
a^{\prime},x)-\hat{p}(m_{2}\mid a^{\prime},x)\right]$

The expression for $\hat{b}(x)$ in Corollary 1 follows by the boundedness of
the propensity-scores, the joint mediator probabilities, and their estimates.
The result then follows from the definition of estimator stability, defined in
(1). ∎

###### Proof of Corollary 2.

This result follows directly from Proposition 4 from Kennedy (2020). ∎

####  B.0.2 Bounds on interventional indirect effect

This section contains proofs pertaining to estimating the proposed bounds on
$\psi_{M_{1}}$.

###### Proof of Proposition 2.

Using the law of iterated expectations we obtain:

$\displaystyle\mathbb{E}[Y^{m_{1}m_{2}}\mid a,x]-\mu_{a}(m_{1},m_{2},x)$ (43)
$\displaystyle=\sum_{m_{1}^{\prime}\neq m_{1},m_{2}^{\prime}\neq
m_{2}}\underbrace{[\mu_{am_{1}m_{2}}^{\star}(m_{1}^{\prime},m_{2}^{\prime},x)-\mu_{a}(m_{1},m_{2},x)]}_{(i)}p(m_{1}^{\prime},m_{2}^{\prime}\mid
a,x)$

The Proposition follows by replacing the term $(i)$ with $b_{l}(\cdot)$ or
$b_{u}(\cdot)$. ∎

###### Proof of Proposition 3.

Recall that we assume throughout that $\tau(x)$ is valid for all values of
$(m_{1},m_{2})$ at $X=x$. The assumptions in (26)-(28) imply:

$\displaystyle(\ref{eqn:a1})\implies b_{l}(\cdot)=-\tau(x),\qquad\qquad
b_{u}(\cdot)=\tau(x)$ $\displaystyle(\ref{eqn:a2})\implies
b_{l}(\cdot)=-\tau(x)\mu_{a}(\cdot),\qquad
b_{u}(\cdot)=[1-\mu_{a}(\cdot)]\tau(x)$ $\displaystyle(\ref{eqn:a3})\implies
b_{l}(\cdot)=-\tau(x)\mu_{a}(\cdot),\qquad
b_{u}(\cdot)=\mu_{a}(\cdot)\tau(x)/(1-\tau(x))$

Notice that all functions $b_{l}$ and $b_{u}$ that take the form:

$\displaystyle h(\cdot)=[c\mu_{a}(\cdot)+t]f(\tau(x))$

for constants $(c,t)\in\\{0,1\\}^{2}$. Replacing $b_{l}$ and $b_{u}$ with
expressions of these form, combining this with the result from Proposition 2,
and then averaging this expression over $[p(m_{1}\mid a,x)-p(m_{1}\mid
a^{\prime},x)]p(m_{2}\mid a^{\prime},x)$ for all values of $(m_{1},m_{2})$
gives the desired result. ∎

###### Proof of Theorem 3.

We focus on the terms $\Gamma_{1,a}$ and $\Gamma_{2,a}$ noting that the other
derivations follow analogously. We treat the data as discrete throughout to
simplify the derivations (see, e.g., Kennedy (2022)).

$\displaystyle\Gamma_{2,a}$
$\displaystyle=\sum_{x,m_{1},m_{2}}p(m_{1},m_{2}\mid a,x)p(m_{1}\mid
a,x)p(m_{2}\mid a^{\prime},x)p(x)$ $\displaystyle\Gamma_{1,a}$
$\displaystyle=\sum_{x,m_{1},m_{2}}\mu_{a}(x,m_{1},m_{2})p(m_{1},m_{2}\mid
a,x)p(m_{1}\mid a,x)p(m_{2}\mid a^{\prime},x)p(x)$

We can use the chain rule to obtain the result:

$\displaystyle\text{IF}(\Gamma_{2,a})$
$\displaystyle=\sum_{x,m_{1},m_{2}}\text{IF}[p(m_{1},m_{2}\mid
a,x)]p(m_{1}\mid a,x)p(m_{2}\mid a^{\prime},x)p(x)$
$\displaystyle+p(m_{1},m_{2}\mid a,x)\text{IF}[p(m_{1}\mid a,x)]p(m_{2}\mid
a^{\prime},x)p(x)$ $\displaystyle+p(m_{1},m_{2}\mid a,x)p(m_{1}\mid
a,x)\text{IF}[p(m_{2}\mid a^{\prime},x)]p(x)$ $\displaystyle+p(m_{1},m_{2}\mid
a,x)p(m_{1}\mid a,x)p(m_{2}\mid a^{\prime},x)\text{IF}[p(x)]$
$\displaystyle=\sum_{x,m_{1},m_{2}}\frac{\mathds{1}(A=a,X=x)}{p(a\mid
x)}[M_{1}M_{2}-p(m_{1},m_{2}\mid a,x)]p(m_{1}\mid a,x)p(m_{2}\mid
a^{\prime},x)$ $\displaystyle+\frac{\mathds{1}(A=a,X=x)}{p(a\mid
x)}[M_{1}-p(m_{1}\mid a,x)]p(m_{1},m_{2}\mid a,x)p(m_{2}\mid a^{\prime},x)$
$\displaystyle+\frac{\mathds{1}(A=a^{\prime},X=x)}{p(a^{\prime}\mid
x)}[M_{2}-p(m_{2}\mid a^{\prime},x)]p(m_{1}\mid a^{\prime},x)p(m_{1},m_{2}\mid
a,x)$ $\displaystyle+p(m_{1},m_{2}\mid a,x)p(m_{1}\mid a,x)p(m_{2}\mid
a^{\prime},x)[\mathds{1}(X=x)-p(x)]$ $\displaystyle=[p(M_{1}\mid
a,X)p(M_{2}\mid a^{\prime},X)-\zeta_{2,a}(X)]$
$\displaystyle+\frac{\mathds{1}(A=a)}{\pi_{a}(X)}[\sum_{m_{2}}p(M_{1},m_{2}\mid
a,X)p(m_{2}\mid a^{\prime},x)-\zeta_{2,a}(X)]$
$\displaystyle+\frac{\mathds{1}(A=a^{\prime})}{\pi_{a^{\prime}}(X)}\sum_{m_{1}}[p(m_{1},M_{2}\mid
a,X)p(m_{1}\mid a,X)-\zeta_{2,a}(X)]+\zeta_{2,a}-\Gamma_{2,a}$

We now show the derivation for $\Gamma_{1,a}$ noting that
$\Gamma_{1,a}=\sum_{x,m_{1},m_{2}}\mu_{a}(x,m_{1},m_{2})\zeta_{2,a}(x)p(x)$.
Therefore,

$\displaystyle\text{IF}(\Gamma_{1,a})$
$\displaystyle=\sum_{x,m_{1},m_{2}}\text{IF}[\mu_{a}(x,m_{1},m_{2})]\zeta_{2,a}(x)p(x)$
$\displaystyle+\mu_{a}(x,m_{1},m_{2})\text{IF}[\zeta_{2,a}(x)]p(x)+\mu_{a}(x,m_{1},m_{2})\zeta_{2,a}(x)\text{IF}[p(x)]$
$\displaystyle=\sum_{x,m_{1},m_{2}}\frac{\mathds{1}(A=a,M_{1}=m_{1},M_{2}=m_{2})}{p(m_{1},m_{2}\mid
a,x)p(a\mid x)}[Y-\mu_{a}(m_{1},m_{2},x)]$ $\displaystyle\times[p(m_{1}\mid
a,x)p(m_{2}\mid a^{\prime},x)p(m_{1},m_{2}\mid a,x)-\zeta_{1,a}]$
$\displaystyle+\mu_{a}(m_{1},m_{2},x)[\frac{\mathds{1}(A=a,X=x)}{p(a\mid
x)}[M_{1}M_{2}-p(m_{1},m_{2}\mid a,x)]p(m_{1}\mid a,x)p(m_{2}\mid
a^{\prime},x)$ $\displaystyle+\frac{\mathds{1}(A=a,X=x)}{p(a\mid
x)}[M_{1}-p(m_{1}\mid a,x)]p(m_{1},m_{2}\mid a,x)p(m_{2}\mid a^{\prime},x)$
$\displaystyle+\frac{\mathds{1}(A=a^{\prime},X=x)}{p(a^{\prime}\mid
x)}[M_{2}-p(m_{2}\mid a^{\prime},x)]p(m_{1}\mid a^{\prime},x)p(m_{1},m_{2}\mid
a,x)$ $\displaystyle+p(m_{1},m_{2}\mid a,x)p(m_{1}\mid a,x)p(m_{2}\mid
a^{\prime},x)[\mathds{1}(X=x)-p(x)]$
$\displaystyle+\mu_{a}(m_{1},m_{2},x)\zeta_{2,a}(x)[\mathds{1}(X=x)-p(x)]$
$\displaystyle=\frac{\mathds{1}(A=a)}{\pi_{a}(X)}[Yp(M_{1}\mid a,X)p(M_{2}\mid
a^{\prime},X)-\zeta_{1,a}(X)]$
$\displaystyle+\frac{\mathds{1}(A=a)}{\pi_{a}(X)}\sum_{m_{2}}[\mu_{a}(M_{1},m_{2},X)p(M_{1},m_{2}\mid,a,X)p(m_{2}\mid
a^{\prime},X)-\zeta_{1,a}(X)$
$\displaystyle+\frac{\mathds{1}(A=a^{\prime})}{\pi_{a^{\prime}}(X)}\left[\sum_{m_{1}}\mu_{a}(m_{1},M_{2},X)p(m_{1},M_{2}\mid
a,X)p(m_{1}\mid a,X)-\zeta_{1,a}(X)\right]$
$\displaystyle+\zeta_{1,a}(X)-\Gamma_{1,a}$

The derivations for $\Gamma_{1,a^{\prime}}$ and $\Gamma_{2,a^{\prime}}$ are
analogous. ∎

###### Proof of Theorem 4.

Let $\hat{\Xi}=\Xi(Z;\hat{\eta},\tau)$. We begin by considering the
decomposition in Kennedy (2022):

$\displaystyle\mathbb{P}_{n}\hat{\Xi}-P\Xi$
$\displaystyle=\underbrace{[\mathbb{P}_{n}-P][\hat{\Xi}(Z)-\Xi(Z)]}_{T_{1}}+\underbrace{[\mathbb{P}_{n}-P]\Xi(Z)}_{T_{2}}+\underbrace{P[\hat{\Xi}(Z)-\Xi(Z)]}_{T_{3}}$

$T_{1}$ is $o_{p}(1/\sqrt{n})$ as long as $\hat{\Xi}$ is consistent for $\Xi$
in the $L_{2}(\mathbb{P})$ norm – i.e. $\|\hat{\Xi}-\Xi\|=o_{p}(1)$ (see,
e.g., Kennedy (2022)). $T_{2}$ is a sample average of a fixed function with
zero mean and therefore by the Central Limit Theorem converges in distribution
to $\mathcal{N}(0,\mathbb{E}[\phi(Z)\phi(Z)^{\top}]$.

It remains to analyze $T_{3}$. We proceed in typical fashion by showing that
this term can be expressed as the product of errors in the nuisance
estimation. We previously showed in the proof of Corollary 1 that
$P[\varphi(Z;\hat{\eta})-\varphi(Z;\eta)]$ is second-order in the nuisance
estimation. It suffices to show that

$\displaystyle
P\\{[\hat{\phi}_{1,a}(Z),\hat{\phi}_{1,a^{\prime}}(Z),\hat{\phi}_{2,a}(Z),\hat{\phi}_{2,a^{\prime}}(Z)]^{\top}-[\phi_{1,a}(Z),\phi_{1,a^{\prime}}(Z),\phi_{2,a}(Z),\phi_{2,a^{\prime}}(Z)]^{\top}\\}$

is also second-order in the nuisance estimation. We show this for each term
separately. We provide the final results below and include the derivations in
Section F.

$\displaystyle P[\hat{\phi}_{1,a}(Z)-\phi_{1,a}(Z)]$
$\displaystyle=\mathbb{E}\left[\frac{\pi_{a}(X)}{\hat{\pi}_{a}(X)}\sum_{m_{1},m_{2}}(p(m_{1},m_{2}\mid
a,x)-\hat{p}(m_{1},m_{2}\mid
a,x))(\mu_{a}(m_{1},m_{2},X)-\hat{\mu}_{a}(m_{1},m_{2},X))\right.$
$\displaystyle\left.\times\hat{p}(m_{1}\mid a,X)\hat{p}(m_{2}\mid
a^{\prime},X)\right]$
$\displaystyle+\left(\frac{\pi_{a}(X)-\hat{\pi}_{a}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}(\mu_{a}(m_{1},m_{2},X)-\hat{\mu}_{a}(m_{1},m_{2},X))\hat{p}(m_{1}\mid
a,X)$ $\displaystyle\times\hat{p}(m_{2}\mid
a^{\prime},X)\hat{p}(m_{1},m_{2}\mid a,X)$
$\displaystyle-\frac{\pi_{a}(X)}{\hat{\pi}_{a}(X)}\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2}\mid
a,X)[p(m_{1},m_{2}\mid a,X)-\hat{p}(m_{1},m_{2}\mid a,X)]$
$\displaystyle\times[p(m_{2}\mid a^{\prime},X)(p(m_{1}\mid
a,X)-\hat{p}(m_{1}\mid a,X))+\hat{p}(m_{1}\mid a,X)(p(m_{2}\mid
a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X))]$
$\displaystyle+\left(\frac{\pi_{a}(X)-\hat{\pi}_{a}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},X)(p(m_{1},m_{2}\mid
a,X)-\hat{p}(m_{1},m_{2}\mid a,X))\hat{p}(m_{1}\mid a,X)\hat{p}(m_{2}\mid
a^{\prime},X)$
$\displaystyle+\left(\frac{\pi_{a}(X)-\hat{\pi}_{a}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},X)\hat{p}(m_{1},m_{2}\mid
a,X)[p(m_{1}\mid a,X)-\hat{p}(m_{1}\mid a,X)]\hat{p}(m_{2}\mid a^{\prime},X)$
$\displaystyle+\left(\frac{\pi_{a^{\prime}}(X)-\hat{\pi}_{a^{\prime}}(X)}{\hat{\pi}_{a^{\prime}}(X)}\right)\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},X)\hat{p}(m_{1},m_{2}\mid
a,X)(p(m_{2}\mid a^{\prime},X)-\hat{p}(m_{2}\mid
a^{\prime},X))\hat{p}(m_{1}\mid a,X)$
$\displaystyle-\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2}\mid
a,X)(p(m_{1},m_{2}\mid a,X)-\hat{p}(m_{1},m_{2}\mid a,X))[(p(m_{1}\mid
a,x)-\hat{p}(m_{1}\mid a,X))]\hat{p}(m_{2}\mid a^{\prime},X)$
$\displaystyle-\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},X)\hat{p}(m_{1},m_{2}\mid
a,X)[p(m_{2}\mid a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X)]$
$\displaystyle\times\\{\hat{p}(m_{1},m_{2}\mid a,X)[p(m_{1}\mid
a,X)-\hat{p}(m_{1}\mid a,X)]+p(m_{1}\mid a,X)[p(m_{1},m_{2}\mid
a,X)-\hat{p}(m_{1},m_{2}\mid a,X)]\\}$
$\displaystyle-\sum_{m_{1},m_{2}}[\mu_{a}(m_{1},m_{2},X)-\hat{\mu}_{a}(m_{1},m_{2},X)]\\{p(m_{1},m_{2}\mid
a,X)p(m_{2}\mid a^{\prime},X)[p(m_{1}\mid a,X)-\hat{p}(m_{1}\mid a,X)]$
$\displaystyle+p(m_{1},m_{2}\mid a,X)\hat{p}(m_{1}\mid a,X)[p(m_{2}\mid
a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X)]$
$\displaystyle+\hat{p}(m_{1}\mid a,X)\hat{p}(m_{2}\mid
a^{\prime},X)[p(m_{1},m_{2}\mid a,X)-\hat{p}(m_{1},m_{2}\mid a,x)]\\}]$

$\displaystyle P[\hat{\phi}_{1,a^{\prime}}(Z)-\phi_{1,a^{\prime}}(Z)]$
$\displaystyle=\mathbb{E}\left[\frac{\pi_{a}(X)}{\hat{\pi}_{a}(X)}\sum_{m_{1},m_{2}}(p(m_{1},m_{2}\mid
a,x)-\hat{p}(m_{1},m_{2}\mid
a,x))(\mu_{a}(m_{1},m_{2},X)-\hat{\mu}_{a}(m_{1},m_{2},X))\right.$
$\displaystyle\left.\times\hat{p}(m_{1}\mid a^{\prime},X)\hat{p}(m_{2}\mid
a^{\prime},X)\right]$
$\displaystyle+\left(\frac{\pi_{a}(X)-\hat{\pi}_{a}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}(\mu_{a}(m_{1},m_{2},X)-\hat{\mu}_{a}(m_{1},m_{2},X))[\hat{p}(m_{1}\mid
a^{\prime},X)]$ $\displaystyle\times\hat{p}(m_{2}\mid
a^{\prime},X)\hat{p}(m_{1},m_{2}\mid a,X)$
$\displaystyle-\frac{\pi_{a}(X)}{\hat{\pi}_{a}(X)}\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2}\mid
a,X)[p(m_{1},m_{2}\mid a,X)-\hat{p}(m_{1},m_{2}\mid a,X)]$
$\displaystyle\times[p(m_{2}\mid a^{\prime},X)(p(m_{1}\mid
a^{\prime},X)-\hat{p}(m_{1}\mid a^{\prime},X))+\hat{p}(m_{1}\mid
a^{\prime},X)(p(m_{2}\mid a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X))]$
$\displaystyle-\left(\frac{\pi_{a}(X)-\hat{\pi}_{a}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},X)(p(m_{1},m_{2}\mid
a,X)-\hat{p}(m_{1},m_{2}\mid a,X))\hat{p}(m_{1}\mid
a^{\prime},X)\hat{p}(m_{2}\mid a^{\prime},X)$
$\displaystyle-\left(\frac{\pi_{a}(X)-\hat{\pi}_{a}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},X)\hat{p}(m_{1},m_{2}\mid
a,X)[(p(m_{1}\mid a^{\prime},X)-\hat{p}(m_{1}\mid
a^{\prime},X))]\hat{p}(m_{2}\mid a^{\prime},X)$
$\displaystyle+\left(\frac{\pi_{a^{\prime}}(X)-\hat{\pi}_{a^{\prime}}(X)}{\hat{\pi}_{a^{\prime}}(X)}\right)\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},X)\hat{p}(m_{1},m_{2}\mid
a,X)(p(m_{2}\mid a^{\prime},X)-\hat{p}(m_{2}\mid
a^{\prime},X))\hat{p}(m_{1}\mid a^{\prime},X)$
$\displaystyle-\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2}\mid
a,X)(p(m_{1},m_{2}\mid a,X)-\hat{p}(m_{1},m_{2}\mid a,X))[p(m_{1}\mid
a^{\prime},x)-\hat{p}(m_{1}\mid a^{\prime},X)]\hat{p}(m_{2}\mid a^{\prime},X)$
$\displaystyle-\sum_{m_{1},m_{2}}\hat{\mu}_{a}(m_{1},m_{2},X)\hat{p}(m_{1},m_{2}\mid
a,X)[p(m_{2}\mid a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X)]$
$\displaystyle\times\\{\hat{p}(m_{1},m_{2}\mid a,X)[p(m_{1}\mid
a^{\prime},X)-\hat{p}(m_{1}\mid a^{\prime},X)]+p(m_{1}\mid
a^{\prime},X)[p(m_{1},m_{2}\mid a,X)-\hat{p}(m_{1},m_{2}\mid a,X)]\\}$
$\displaystyle-\sum_{m_{1},m_{2}}[\mu_{a}(m_{1},m_{2},X)-\hat{\mu}_{a}(m_{1},m_{2},X)]\\{p(m_{1},m_{2}\mid
a,X)p(m_{2}\mid a^{\prime},X)[p(m_{1}\mid a^{\prime},X)-\hat{p}(m_{1}\mid
a^{\prime},X)]$ $\displaystyle+p(m_{1},m_{2}\mid a,X)\hat{p}(m_{1}\mid
a^{\prime},X)[p(m_{2}\mid a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X)]$
$\displaystyle+\hat{p}(m_{1}\mid a^{\prime},X)\hat{p}(m_{2}\mid
a^{\prime},X)[p(m_{1},m_{2}\mid a,X)-\hat{p}(m_{1},m_{2}\mid a,x)]\\}]$

$\displaystyle P[\hat{\phi}_{2,a}(Z)-\phi_{2,a}(Z)]$
$\displaystyle=\mathbb{E}\left[\left(\frac{\pi_{a}(X)-\hat{\pi}_{a}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}[p(m_{1}\mid
a,X)-\hat{p}(m_{1}\mid a,X)]\hat{p}(m_{2}\mid
a^{\prime},X)\hat{p}(m_{1},m_{2}\mid a,X)\right.$
$\displaystyle\left.+\left(\frac{\pi_{a^{\prime}}(X)-\hat{\pi}_{a^{\prime}}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}[p(m_{2}\mid
a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X)]\hat{p}(m_{1}\mid
a,X)\hat{p}(m_{1},m_{2}\mid a,X)\right.$
$\displaystyle\left.+\left(\frac{\pi_{a}(X)-\hat{\pi}_{a}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}[p(m_{1},m_{2}\mid
a,X)-\hat{p}(m_{1},m_{2}\mid a,X)]\hat{p}(m_{1}\mid a,X)\hat{p}(m_{2}\mid
a^{\prime},X)\right.$
$\displaystyle\left.-\sum_{m_{1},m_{2}}\hat{p}(m_{1},m_{2}\mid
a,X)[p(m_{1}\mid a,X)-\hat{p}(m_{1}\mid a,X)][p(m_{2}\mid
a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X)]\right.$
$\displaystyle\left.-\sum_{m_{1},m_{2}}\hat{p}(m_{1}\mid
a,x)[p(m_{1},m_{2}\mid a,X)-\hat{p}(m_{1},m_{2}\mid a,X)][p(m_{2}\mid
a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X)]\right.$
$\displaystyle\left.-\sum_{m_{1},m_{2}}p(m_{2}\mid
a^{\prime},X)[p(m_{1},m_{2}\mid a,X)-\hat{p}(m_{1},m_{2}\mid a,X)][p(m_{1}\mid
a,X)-\hat{p}(m_{1}\mid a,X)]\right]$

$\displaystyle P[\hat{\phi}_{2,a^{\prime}}(Z)-\phi_{2,a^{\prime}}(Z)]$
$\displaystyle=\mathbb{E}\left[\left(\frac{\pi_{a^{\prime}}(X)-\hat{\pi}_{a^{\prime}}(X)}{\hat{\pi}_{a^{\prime}}(X)}\right)\sum_{m_{1},m_{2}}[p(m_{1}\mid
a^{\prime},X)-\hat{p}(m_{1}\mid a^{\prime},X)]\hat{p}(m_{2}\mid
a^{\prime},X)\hat{p}(m_{1},m_{2}\mid a,X)\right.$
$\displaystyle\left.+\left(\frac{\pi_{a^{\prime}}(X)-\hat{\pi}_{a^{\prime}}(X)}{\hat{\pi}_{a^{\prime}}(X)}\right)\sum_{m_{1},m_{2}}[p(m_{2}\mid
a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X)]\hat{p}(m_{1}\mid
a^{\prime},X)\hat{p}(m_{1},m_{2}\mid a,X)\right.$
$\displaystyle\left.+\left(\frac{\pi_{a}(X)-\hat{\pi}_{a}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}[p(m_{1},m_{2}\mid
a,X)-\hat{p}(m_{1},m_{2}\mid a,X)]\hat{p}(m_{1}\mid
a^{\prime},x)\hat{p}(m_{2}\mid a^{\prime},X)\right.$
$\displaystyle\left.-\sum_{m_{1},m_{2}}\hat{p}(m_{1}\mid
a^{\prime},x)[p(m_{1},m_{2}\mid a,X)-\hat{p}(m_{1},m_{2}\mid a,X)][p(m_{2}\mid
a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X)]\right.$
$\displaystyle\left.-\sum_{m_{1},m_{2}}\hat{p}(m_{1},m_{2}\mid
a,x)[p(m_{1}\mid a^{\prime},X)-\hat{p}(m_{1}\mid a^{\prime},X)][p(m_{2}\mid
a^{\prime},X)-\hat{p}(m_{2}\mid a^{\prime},X)]\right.$
$\displaystyle\left.-\sum_{m_{1},m_{2}}p(m_{2}\mid
a^{\prime},X)[p(m_{1},m_{2}\mid a,X)-\hat{p}(m_{1},m_{2}\mid a,X)][p(m_{1}\mid
a^{\prime},X)-\hat{p}(m_{1}\mid a^{\prime},X)]\right]$

Each of the terms above can be shown to be bounded by the product of nuisance
estimation times a constant, having assumed that the propensity scores, the
joint mediator probabilities, and their estimates are bounded away from zero
(and one for the propensity scores) by $\epsilon$. For example, consider the
first term in the decomposition of $P[\hat{\phi}_{2,a}(Z)-\phi_{2,a}(Z)]$:

$\displaystyle\mathbb{E}\left[\left(\frac{\pi_{a}(X)-\hat{\pi}_{a}(X)}{\hat{\pi}_{a}(X)}\right)\sum_{m_{1},m_{2}}[p(m_{1}\mid
a,X)-\hat{p}(m_{1}\mid a,X)]\hat{p}(m_{2}\mid
a^{\prime},X)\hat{p}(m_{1},m_{2}\mid a,X)\right]$

By the boundedness of the propensity score estimates, we obtain that this
term:

$\displaystyle\leq\frac{1}{\epsilon}P\left[\frac{\mathds{1}(A=a)[\pi_{a}(X)-\hat{\pi}_{a}(X)][p(M_{1}\mid
a,X)-\hat{p}(M_{1}\mid a,X)]\hat{p}(M_{2}\mid
a^{\prime},X)\hat{p}(M_{1},M_{2}\mid a,X)}{p(A\mid X)p(M_{1},M_{2}\mid
A,X)}\right]$
$\displaystyle\leq\frac{1}{\epsilon^{3}}\|\pi_{a}(X)-\hat{\pi}_{a}(X)\|\|p(M_{1}\mid
a,X)-\hat{p}(M_{1}\mid a,X)\|$

where the final line holds by the boundedness of the propensity scores,
mediator probabilities, and their estimates, and by Cauchy-Schwarz. We can
make similar derivations for the remaining terms. Therefore, as long as
$Z_{n}=o_{p}(n^{-1/2})$, where

$\displaystyle Z_{n}$
$\displaystyle=\|\hat{\mu}_{a}(M_{1},M_{2},X)-\mu_{a}(M_{1},M_{2},X)\|\left[\|\hat{\pi}(X)_{a}-\pi_{a}(X)\|\right.$
$\displaystyle\left.+\|\hat{p}(M_{1}\mid a,X)-p(M_{1}\mid
a,X)\|+\|\hat{p}(M_{1},M_{2}\mid a,X)-p(M_{1},M_{2}\mid a,X)\|\right.$
$\displaystyle\left.+\|\hat{p}(M_{1}\mid a^{\prime},X)-p(M_{1}\mid
a^{\prime},X)\|+\|\hat{p}(M_{2}\mid a^{\prime},X)-p(M_{2}\mid
a^{\prime},X)\|\right]$
$\displaystyle+\|\hat{\pi}_{a}(X)-\pi_{a}(X)\|\left[\|\hat{p}(M_{1}\mid
a,X)-p(M_{1}\mid a,X)\|+\|\hat{p}(M_{1}\mid a^{\prime},X)-p(M_{1}\mid
a^{\prime},X)\|\right.$ $\displaystyle\left.+\|\hat{p}(M_{2}\mid
a^{\prime},X)-p(M_{2}\mid a^{\prime},X)\|+\|\hat{p}(M_{1},M_{2}\mid
a,X)-p(M_{1},M_{2}\mid a,X)\|\right]$ $\displaystyle+\|\hat{p}(M_{1},M_{2}\mid
a,X)-p(M_{1},M_{2}\mid a,X)\|\left[\|\hat{p}(M_{1}\mid a,X)-p(M_{1}\mid
a,X)\|\right.$ $\displaystyle\left.+\|\hat{p}(M_{1}\mid
a^{\prime},X)-p(M_{1}\mid a^{\prime},X)\|+\|\hat{p}(M_{2}\mid
a^{\prime},X)-p(M_{2}\mid a^{\prime},X)\|\right]$
$\displaystyle+\|\hat{p}(M_{2}\mid a^{\prime},X)-p(M_{2}\mid
a^{\prime},X)\|[\|\hat{p}(M_{1}\mid a,X)-p(M_{1}\mid a,X)\|$
$\displaystyle+\|\hat{p}(M_{1}\mid a^{\prime},X)-p(M_{1}\mid a^{\prime},X)\|]$

then the result follows. This would be satisfied, if, for example, each
element of $\eta$ were estimated at a rate of $o_{p}(n^{-1/4})$. ∎

######  Remark 9.

Estimating each element of $\eta$ at rates of $o_{p}(n^{-1/4})$ suffices to
guarantee both asymptotic normality of the point estimate of $\psi_{M_{1}}$
and the bounds above.

###### Proof of Corollary 3.

The proof of this corollary for either the projection estimator of the bounds
is identical to the proof of Theorem 1, where the result follows from the fact
that:

$\displaystyle P[\phi(Z;\hat{\eta},\tau)-\phi(Z;\eta,\tau)]$

is second-order in the nuisance estimation, with the second-order expression
provided by the term $Z_{n}$ in the statement of Corollary 3. ∎

###### Proof of Corollary 4.

The proof of this corollary relies primarily on the fact that the conditional
bias term at $X=x$ is second-order in the nuisance estimation. This follows
using the decomposition from the proof of Theorem 4 and following the steps of
Corollaries 1 and 2.

We first derive the following expressions:

$\displaystyle\hat{k}_{1,a}^{\star}(x)=\mathbb{E}[\hat{\phi}_{1,a}-\phi_{1,a}\mid
X=x,D_{1}^{n}]$
$\displaystyle\hat{k}_{1,a^{\prime}}^{\star}(x)=\mathbb{E}[\hat{\phi}_{1,a^{\prime}}-\phi_{1,a^{\prime}}\mid
X=x,D_{1}^{n}]$