#  Semiparametric Estimation for Causal Mediation Analysis with Multiple
Causally Ordered Mediators††thanks: Direct all correspondence to Xiang Zhou,
Department of Sociology, Harvard University, 33 Kirkland Street, Cambridge MA
02138; email<EMAIL_ADDRESS>The author thanks the Editor, the
Associate Editor, two anonymous reviewers, two reviewers from the Alexander
and Diviya Magaro Peer Pre-Review Program, and Aleksei Opacic for helpful
comments.

Xiang Zhou
Harvard University

###### Abstract

Causal mediation analysis concerns the pathways through which a treatment
affects an outcome. While most of the mediation literature focuses on settings
with a single mediator, a flourishing line of research has examined settings
involving multiple mediators, under which path-specific effects (PSEs) are
often of interest. We consider estimation of PSEs when the treatment effect
operates through $K(\geq 1)$ causally ordered, possibly multivariate
mediators. In this setting, the PSEs for many causal paths are not
nonparametrically identified, and we focus on a set of PSEs that are
identified under Pearl’s nonparametric structural equation model. These PSEs
are defined as contrasts between the expectations of $2^{K+1}$ potential
outcomes and identified via what we call the generalized mediation functional
(GMF). We introduce an array of regression-imputation, weighting, and “hybrid”
estimators, and, in particular, two $K+2$-robust and locally semiparametric
efficient estimators for the GMF. The latter estimators are well suited to the
use of data-adaptive methods for estimating their nuisance functions. We
establish the rate conditions required of the nuisance functions for
semiparametric efficiency. We also discuss how our framework applies to
several estimands that may be of particular interest in empirical
applications. The proposed estimators are illustrated with a simulation study
and an empirical example.

Keywords: causal inference, mediation, path-specific effects, multiple
robustness, semiparametric efficiency

###  1 Introduction

Causal mediation analysis aims to disentangle the pathways through which a
treatment affects an outcome. While traditional approaches to mediation
analysis have relied on linear structural equation models, along with their
stringent parametric assumptions, to define and estimate direct and indirect
effects (e.g., Baron and Kenny 1986), a large body of research has emerged
within the causal inference literature that disentangles the tasks of
definition, identification, and estimation in the study of causal mechanisms.
Using the potential outcomes framework (Neyman 1923; Rubin 1974), this body of
research has provided model-free definitions of direct and indirect effects
(Robins and Greenland 1992; Pearl 2001), established the assumptions needed
for nonparametric identification (Robins and Greenland 1992; Pearl 2001;
Robins 2003; Petersen _et al._ 2006; Imai _et al._ 2010; Hafeman and
VanderWeele 2011; VanderWeele 2015), and developed an array of imputation,
weighting, and multiply robust methods for estimation (e.g., Goetgeluk _et
al._ 2009; Albert 2012; Tchetgen Tchetgen and Shpitser 2012; Vansteelandt _et
al._ 2012; Zheng and van der Laan 2012; Tchetgen Tchetgen 2013; VanderWeele
2015; Wodtke and Zhou 2020).

While the bulk of the causal mediation literature focuses on settings with a
single mediator (or a set of mediators considered as a whole), a flourishing
line of research has studied settings that involve multiple causally dependent
mediators, under which a set of path-specific effects (PSEs) are often of
interest (Avin _et al._ 2005; Albert and Nelson 2011; Shpitser 2013;
VanderWeele and Vansteelandt 2014; VanderWeele _et al._ 2014; Daniel _et al._
2015; Lin and VanderWeele 2017; Miles _et al._ 2017; Steen _et al._ 2017;
Vansteelandt and Daniel 2017; Miles _et al._ 2020). In particular, Daniel _et
al._ (2015) demonstrated a large number of ways in which the total effect of a
treatment can be decomposed into PSEs, established the assumptions under which
a subset of these PSEs are identified, and provided a parametric method for
estimating these effects (see also Albert and Nelson 2011). More recently, for
a particular PSE in the case of two causally ordered mediators, Miles _et al._
(2020) offered an in-depth discussion of alternative estimation methods, and,
utilizing the efficient influence function of its identification formula,
developed a triply robust and locally semiparametric efficient estimator. This
estimator, by virtue of its multiple robustness, is well suited to the use of
data-adaptive methods for estimating its nuisance functions.

To date, most of the literature on PSEs has focused on the case of two
mediators, and it remains underexplored how the estimation methods developed
in previous studies, such as those in VanderWeele _et al._ (2014) and Miles
_et al._ (2020), generalize to the case of $K(\geq 1)$ causally ordered
mediators. This article aims to bridge this gap. First, we observe that
despite a multitude of ways in which a PSE can be defined for each causal path
from the treatment to the outcome, most of these PSEs are not identified under
Pearl’s nonparametric structural equation model. This observation leads us to
focus on the much smaller set of PSEs that can be nonparametrically
identified. These PSEs are defined as contrasts between the expectations of
$2^{K+1}$ potential outcomes, which, in turn, are identified through a formula
that can be viewed as an extension of Pearl’s (2001) and Daniel _et al._ ’s
(2015) mediation formulae to the case of $K$ causally ordered mediators.
Following Tchetgen Tchetgen and Shpitser (2012), we refer to the
identification formula for these expected potential outcomes as the
generalized mediation functional (GMF).

We then show that the GMF can be estimated via an array of regression,
weighting, and “hybrid” estimators. More important, building on its efficient
influence function (EIF), we develop two multiply robust and locally
semiparametric efficient estimators for the GMF. Both of these estimators are
$K+2$-robust, in the sense that they are consistent provided that one of $K+2$
sets of nuisance functions is correctly specified and consistently estimated.
These multiply robust estimators are well suited to the use of data-adaptive
methods for estimating the nuisance functions. We establish rate conditions
for consistency and semiparametric efficiency when data-adaptive methods and
cross-fitting (Zheng and van der Laan 2011; Chernozhukov _et al._ 2018) are
used to estimate the nuisance functions.

Compared with existing estimators that have been proposed for causal mediation
analysis, the methodology proposed in this article is distinct in its
generality. In fact, the doubly robust estimator for the mean of an incomplete
outcome (Scharfstein _et al._ 1999), the triply robust estimator developed by
Tchetgen Tchetgen and Shpitser (2012) for the mediation functional in the one-
mediator setting (see also Zheng and van der Laan 2012), and the estimator
proposed by Miles _et al._ (2020) for their particular PSE, can all be viewed
as special cases of the $K+2$-robust estimators — when $K=0,1,2$,
respectively. Yet, our framework also encompasses important estimands for
which semiparametric estimators have not been proposed. To demonstrate the
generality of our framework, we show how our multiply robust semiparametric
estimators apply to several estimands that may be of particular interest in
empirical applications, including the natural direct effect (NDE), the
natural/total indirect effect (NIE/TIE), the natural path-specific effect
(nPSE), and the cumulative path-specific effect (cPSE). In Supplementary
Material E, we discuss how our framework can also be employed to estimate
noncausal decompositions of between-group disparities that are widely used in
social science research (Fortin _et al._ 2011).

Before proceeding, we note that in a separate strand of literature, the term
“multiple robustness” has been used to characterize a class of estimators for
the mean of incomplete data that are consistent if one of several working
models for the propensity score or one of several working models for the
outcome is correctly specified (e.g., Han and Wang 2013; Han 2014). In this
paper, we use “$V$-robustness” to characterize estimators that require
modeling multiple parts of the observed data likelihood and are consistent
provided that one of $V$ sets of the corresponding models is correctly
specified, in keeping with the terminology in the causal mediation literature.
This definition of “multiple robustness” does not imply that a “$K+2$-robust”
estimator is necessarily more robust than, for example, a “$K+1$-robust”
estimator. First, they may correspond to different estimands that require
modeling different parts of the likelihood. For example, the doubly robust
estimator of the average treatment effect only involves a propensity score
model and an outcome model; it is thus less demanding than Tchetgen Tchetgen
and Shpitser’s (2012) triply robust estimator of the mediation functional,
which involves an additional model for the mediator. Second, for our
semiparametric estimators of the GMF, the “$K+2$-robustness” property is not
“sharp” because it can be tightened in various special cases. As we
demonstrate in Section 4 and Supplementary Material E, such a tightening may
result in a lower $V$ (as in the case of NDE, NIE/TIE, nPSE, and cPSE), or a
higher $V$ (as in the case of noncausal decompositions of between-group
disparities).

The rest of the paper is organized as follows. In Section 2, we define the
PSEs of interest, lay out their identification assumptions, and introduce the
GMF. In Section 3, we introduce a range of regression-imputation, weighting,
“hybrid,” and multiply robust estimators for the GMF, and present several
techniques that could be used to improve the finite sample performance of the
multiply robust estimators. In Section 4, we discuss how our results apply to
a number of special cases such as the NDE, NIE/TIE, nPSE, and cPSE. A
simulation study and an empirical example are given in Section 5 and Section 6
to illustrate the proposed estimators. Proofs of Theorems 1-4 are given in
Supplementary Materials A, C, and D. Replication data and code for the
simulation study and the empirical example are available at
https://doi.org/10.7910/DVN/5TBUM3.

$X$$A$$M_{2}$$M_{1}$$Y$

(a)$A$$M_{2}$$M_{1}$$Y$(c)$A$$M_{2}$$M_{1}$$Y$

(b)$A$$M_{2}$$M_{1}$$Y$(d)$A$$M_{2}$$M_{1}$$Y$

Figure 1: Causal relationships with two causally ordered mediators.

Note: $A$ denotes the treatment, $Y$ denotes the outcome of interest, $X$
denotes a vector of pretreatment covariates, and $M_{1}$ and $M_{2}$ denote
two causally ordered mediators.

###  2 Notation, Definitions, and Identification

To ease exposition, we start with the case of two causally ordered mediators
before moving onto the general setting of $K$ mediators.

#### 2.1 The Case of Two Causally Ordered Mediators

Let $A$ denote a binary treatment, $Y$ an outcome of interest, and $X$ a
vector of pretreatment covariates. In addition, let $M_{1}$ and $M_{2}$ denote
two causally ordered mediators, and assume $M_{1}$ precedes $M_{2}$. We allow
each of these mediators to be multivariate, in which case the causal
relationships among the component variables are left unspecified. A directed
acyclic graph (DAG) representing the relationships between these variables is
given in the top panel of Figure 1. In this DAG, four possible causal paths
exist from the treatment to the outcome, as shown in the lower panels: (a)
$A\to Y$; (b) $A\to M_{2}\to Y$; (c) $A\to M_{1}\to Y$; and (d) $A\to M_{1}\to
M_{2}\to Y$.

A formal definition of path-specific effects (PSEs) requires the potential-
outcomes notation for both the outcome and the mediators. Specifically, let
$Y(a,m_{1},m_{2})$ denote the potential outcome under treatment status $a$ and
mediator values $M_{1}=m_{1}$ and $M_{2}=m_{2}$, $M_{2}(a,m_{1})$ the
potential value of the mediator $M_{2}$ under treatment status $a$ and
mediator value $M_{1}=m_{1}$, and $M_{1}(a)$ the potential value of the
mediator $M_{1}$ under treatment status $a$. This notation allows us to define
nested counterfactuals in the form of
$Y\big{(}a,M_{1}(a_{1}),M_{2}(a_{2},M_{1}(a_{12}))\big{)}$, where $a$,
$a_{1}$, $a_{2}$, and $a_{12}$ can each take 0 or 1. For example,
$Y\big{(}1,M_{1}(0),M_{2}(0,M_{1}(0))\big{)}$ represents the potential outcome
in the hypothetical scenario where the subject was treated but the mediators
$M_{1}$ and $M_{2}$ were set to values they would have taken if the subject
had not been treated. Further, if we let $Y(a)$ denote the potential outcome
when treatment status is set to $a$ and the mediators $M_{1}$ and $M_{2}$ take
on their “natural” values under treatment status $a$ (i.e., $M_{1}(a)$ and
$M_{2}(a,M_{1}(a))$), we have
$Y(a)=Y\big{(}a,M_{1}(a),M_{2}(a,M_{1}(a))\big{)}$ by construction. This is
sometimes referred to as the “composition” assumption (VanderWeele 2009).

Under the above notation, for each of the causal paths shown in Figure 1, its
PSE can be defined in eight different ways, depending on the reference levels
chosen for $A$ for each of the other three paths (Daniel _et al._ 2015). For
example, the average direct effect of $A$ on $Y$, i.e., the portion of the
treatment effect that operates through the path $A\to Y$, can be defined as

$\tau_{A\to
Y}(a_{1},a_{2},a_{12})=\mathbb{E}[Y\big{(}1,M_{1}(a_{1}),M_{2}(a_{2},M_{1}(a_{12}))\big{)}-Y\big{(}0,M_{1}(a_{1}),M_{2}(a_{2},M_{1}(a_{12}))\big{)}],$

where $a_{1}$, $a_{2}$, and $a_{12}$ can each take 0 or 1. In particular,
$\tau_{A\to Y}(0,0,0)$ corresponds to the natural direct effect (NDE; Pearl
2001) or pure direct effect (PDE; Robins and Greenland 1992) if the mediators
$M_{1}$ and $M_{2}$ are considered as a whole. In a similar vein, the PSEs via
$A\to M_{2}\to Y$, $A\to M_{1}\to Y$, and $A\to M_{1}\to M_{2}\to Y$ can be
defined as

$\displaystyle\tau_{A\to M_{2}\to Y}(a,a_{1},a_{12})$
$\displaystyle=\mathbb{E}[Y\big{(}a,M_{1}(a_{1}),M_{2}(1,M_{1}(a_{12}))\big{)}-Y\big{(}a,M_{1}(a_{1}),M_{2}(0,M_{1}(a_{12}))\big{)}],$
$\displaystyle\tau_{A\to M_{1}\to Y}(a,a_{2},a_{12})$
$\displaystyle=\mathbb{E}[Y\big{(}a,M_{1}(1),M_{2}(a_{2},M_{1}(a_{12}))\big{)}-Y\big{(}a,M_{1}(0),M_{2}(a_{2},M_{1}(a_{12}))\big{)}],$
$\displaystyle\tau_{A\to M_{1}\to M_{2}\to Y}(a,a_{1},a_{2})$
$\displaystyle=\mathbb{E}[Y\big{(}a,M_{1}(a_{1}),M_{2}(a_{2},M_{1}(1))\big{)}-Y\big{(}a,M_{1}(a_{1}),M_{2}(a_{2},M_{1}(0))\big{)}].$

In addition, if we use $A\to M_{1}\rightsquigarrow Y$ to denote the
combination of the causal paths $A\to M_{1}\to Y$ and $A\to M_{1}\to M_{2}\to
Y$, the corresponding PSE for this “composite path” can be defined as

$\tau_{A\to M_{1}\rightsquigarrow
Y}(a,a_{2})=\mathbb{E}[Y\big{(}a,M_{1}(1),M_{2}(a_{2},M_{1}(1))\big{)}-Y\big{(}a,M_{1}(0),M_{2}(a_{2},M_{1}(0))\big{)}].$

This quantity reflects the portion of the treatment effect that operates
through $M_{1}$, regardless of whether it further operates through $M_{2}$ or
not. In particular, $\tau_{A\to M_{1}\rightsquigarrow Y}(0,0)$ is often
referred to as the natural indirect effect (NIE; Pearl 2001) or the pure
indirect effect (PIE; Robins and Greenland 1992) for $M_{1}$, whereas
$\tau_{A\to M_{1}\rightsquigarrow Y}(1,1)$ is sometimes called the total
indirect effect (TIE; Robins and Greenland 1992) for $M_{1}$. Note, however,
that the term NIE has also been used to denote $\tau_{A\to
M_{1}\rightsquigarrow Y}(1,1)$ (e.g., Tchetgen Tchetgen and Shpitser 2012). To
avoid ambiguity, we use NIE and TIE to denote $\tau_{A\to
M_{1}\rightsquigarrow Y}(0,0)$ and $\tau_{A\to M_{1}\rightsquigarrow Y}(1,1)$,
respectively. By definition, these PSEs are identified if the corresponding
expected potential outcomes, i.e.,
$\mathbb{E}[Y\big{(}a,M_{1}(a_{1}),M_{2}(a_{2},M_{1}(a_{12}))\big{)}]$, are
identified. Below, we review the assumptions under which these expected
potential outcomes are identified from observed data.

Following Pearl (2009), we use a DAG to encode a nonparametric structural
equation model (NPSEM) with mutually independent errors. In this framework,
the top panel of Figure 1 implies no unobserved confounding for any of the
treatment-mediator, treatment-outcome, mediator-mediator, and mediator-outcome
relationships. Formally, we invoke the following assumptions.

######  Assumption 1.

Consistency of $A$ on $M_{1}$, $(A,M_{1})$ on $M_{2}$, and $(A,M_{1},M_{2})$
on $Y$: For any unit and any $a,m_{1},m_{2}$, $M_{1}=M_{1}(a)$ if $A=a$;
$M_{2}=M_{2}(a,m_{1})$ if $A=a$ and $M_{1}=m_{1}$; and $Y=Y(a,m_{1},m_{2})$ if
$A=a$, $M_{1}=m_{1}$, and $M_{2}=m_{2}$.

######  Assumption 2.

Conditional independence among treatment and potential outcomes: for any
$a,a_{1},a_{2},m_{1},m_{1}^{*},m_{2}$,
$\big{(}M_{1}(a_{1}),M_{2}(a_{2},m_{1}),Y(a,m_{1},m_{2})\big{)}\perp\\!\\!\\!\perp
A|X$; $\big{(}M_{2}(a_{2},m_{1}),Y(a,m_{1},m_{2})\big{)}\perp\\!\\!\\!\perp
M_{1}(a_{1})|X,A$, and $Y(a,m_{1},m_{2})\perp\\!\\!\\!\perp
M_{2}(a_{2},m_{1}^{*})|X,A,M_{1}$.

######  Assumption 3.

Positivity: $p_{A|X}(a|x)>\epsilon>0$ whenever $p_{X}(x)>0$;
$p_{A|X,M_{1}}(a|x,m_{1})>\epsilon>0$ whenever $p_{X,M_{1}}(x,m_{1})>0$, and
$p_{A|X,M_{1},M_{2}}(a|x,m_{1},m_{2})>\epsilon>0$ whenever
$p_{X,M,M_{2}}(x,m_{1},m_{2})>0$, where $p(\cdot)$ denotes a probability
density/mass function.

Note that Assumption 2 involves conditional independence relationships between
the so-called cross-world counterfactuals, such as
$\big{(}M_{2}(a_{2},m_{1}),Y(a,m_{1},m_{2})\big{)}\perp\\!\\!\\!\perp
M_{1}(a_{1})|X,A$. This assumption is a direct consequence of Pearl’s NPSEM
with mutually independent errors. It implies, but is not implied by, the
sequential ignorability assumption that Robins (2003) invokes in interpreting
causal diagrams (see Robins and Richardson 2010 for an in-depth discussion).
In addition, we note that Assumption 2 does not rule out all forms of
unobserved confounding for the causal effects of $X$ on its descendants. For
example, unobserved variables are permitted (although not shown) in Figure 1
that affect both $X$ and $Y$.

Under Assumptions 1-3, it can be shown that
$\mathbb{E}[Y\big{(}a,M_{1}(a_{1}),M_{2}(a_{2},M_{1}(a_{12}))\big{)}]$ is
identified if and only if $a_{12}=a_{1}$ (Avin _et al._ 2005; Albert and
Nelson 2011; Daniel _et al._ 2015). Consequently, none of the PSEs for the
path $A\to M_{1}\to Y$ is identified because given $a_{12}$, either
$\mathbb{E}[Y\big{(}a,M_{1}(1),M_{2}(a_{2},M_{1}(a_{12}))\big{)}]$ or
$\mathbb{E}[Y\big{(}a,M_{1}(0),M_{2}(a_{2},M_{1}(a_{12}))\big{)}]$ is
unidentified. Similarly, none of the PSEs for the path $A\to M_{1}\to M_{2}\to
Y$ is identified. Interestingly, the PSEs for the composite path $A\to
M_{1}\rightsquigarrow Y$ are all identified, even if $a\neq a_{2}$. These
results echo the recanting witness criterion developed by Avin _et al._
(2005), which implies that the PSE for a (possibly composite) path from $A$ to
$Y$ when $A$ is set to 0 (or 1) for all other paths is identified if and only
if the path of interest contains no “recanting witness” — a variable $W$ that
has an additional path to $Y$ that is not contained in the path of interest.
Thus the PSE $\tau_{A\to M_{1}\to Y}(0,0,0)$ is not identified because $M_{1}$
has an additional path to $Y$ ($M_{1}\to M_{2}\to Y$) that is not contained in
$A\to M_{1}\to Y$, but the PSE $\tau_{A\to M_{1}\rightsquigarrow Y}(0,0)$ is
identified because all possible paths from $M_{1}$ to $Y$ is contained in
$A\to M_{1}\rightsquigarrow Y$.

Because $\mathbb{E}[Y\big{(}a,M_{1}(a_{1}),M_{2}(a_{2},M_{1}(a_{12}))\big{)}]$
is identified if and only if $a_{1}=a_{12}$, we restrict our attention to
cases where $a_{1}=a_{12}$ and use the following notation:

$\psi_{a_{1},a_{2},a}\stackrel{{\scriptstyle\Delta}}{{=}}\mathbb{E}\big{[}Y\big{(}a,M_{1}(a_{1}),M_{2}(a_{2},M_{1}(a_{1}))\big{)}\big{]}.$

Under Assumptions 1-3, $\psi_{a_{1},a_{2},a}$ is identified via the following
formula:

$\displaystyle\psi_{a_{1},a_{2},a}$
$\displaystyle=\iiint\mathbb{E}[Y|x,a,m_{1},m_{2}]dP(m_{2}|x,a_{2},m_{1})dP(m_{1}|x,a_{1})dP(x).$
(1)

For a proof of the above formula, see Daniel _et al._ (2015). Equation (1) can
be seen as an extension of Pearl’s (2001) mediation formula to the case of two
causally ordered mediators.

It should be noted that Assumptions 1-3 constitute a sufficient set of
conditions that allow us to identify $\psi_{a_{1},a_{2},a}$ for arbitrary
combinations of $a_{1}$, $a_{2}$, and $a$. For specific combinations of
$a_{1}$, $a_{2}$, and $a$, Assumption 2 can be relaxed. For example,
$\psi_{100}$ is still identified via equation (1) when unobserved confounding
exists for the $M_{2}$-$Y$ relationship, and $\psi_{010}$ is still identified
via equation (1) when unobserved confounding exists for the $M_{1}$-$Y$
relationship (Shpitser 2013; Miles _et al._ 2020).

#### 2.2 The Case of $K(\geq 1)$ Causally Ordered Mediators

We now generalize the preceding results to the setting where the treatment
effect of $A$ on $Y$ operates through $K$ causally ordered, possibly
multivariate mediators, $M_{1},M_{2},\ldots M_{K}$. We assume that for any
$k<k^{\prime}$, $M_{k}$ precedes $M_{k^{\prime}}$, such that no component of
$M_{k^{\prime}}$ causally affects any component of $M_{k}$. In a DAG that is
consistent with this setup, a directed path from the treatment to the outcome
can pass through any combination of the $K$ mediators, resulting in $2^{K}$
possible paths. Among the $2^{K}$ paths, each can be switched “on” or “off,”
creating $2^{2^{K}}$ potential outcomes. Also, for each of the $2^{K}$ paths,
the corresponding PSE can be defined in $2^{2^{K}-1}$ different ways,
depending on whether each of the other $2^{K}-1$ paths is switched “on” or
“off.” For example, when $K=3$, for each causal path from $A$ to $Y$, its PSE
can be defined in $2^{2^{3}-1}=128$ different ways.

As we will see, despite the exponential growth of possible causal paths and
the double exponential growth of possible PSEs, most of these PSEs are not
identified under the assumptions associated with Pearl’s NPSEM. To fix ideas,
let an overbar denote a vector of variables, so that
$\overline{M}_{k}=(M_{1},M_{2},\ldots M_{k})$,
$\overline{m}_{k}=(m_{1},m_{2},\ldots m_{k})$, and
$\overline{a}_{k}=(a_{1},a_{2},\ldots a_{k})$, where
$\overline{M}_{l}=\overline{m}_{l}=\overline{a}_{l}=\varnothing$ if $l\leq 0$.
In addition, let $[K]$ denote the set $\\{1,2,\ldots K\\}$, and let $a_{K+1}$,
instead of $a$, denote the treatment status set to the path $A\to Y$.
Assumptions 1-3 can now be generalized as below.

######  Assumption 1∗.

Consistency: For any unit, $M_{k}=M_{k}(a_{k},\overline{m}_{k-1})$ if
$A=a_{k}$ and $\overline{M}_{k-1}=\overline{m}_{k-1}$, $\forall k\in[K]$; and
$Y=Y(a_{K+1},\overline{m}_{K})$ if $A=a_{K+1}$ and
$\overline{M}_{K}=\overline{m}_{K}$.

######  Assumption 2∗.

Conditional independence among treatment and potential outcomes:
$\big{(}M_{1}(a_{1}),M_{2}(a_{2},\overline{m}_{1}),\ldots
Y(a_{K+1},\overline{m}_{K})\big{)}\perp\\!\\!\\!\perp A|X$; and
$\big{(}M_{k+1}(a_{k+1},\overline{m}_{k}),\ldots
M_{K}(a_{K},\overline{m}_{K-1}),Y(a_{K+1},\overline{m}_{K})\big{)}\perp\\!\\!\\!\perp
M_{k}(a_{k},\overline{m}_{k-1}^{*})|X,A,\overline{M}_{k-1},\forall k\in[K]$.

######  Assumption 3∗.

Positivity: $p_{A|X}(a|x)>\epsilon>0$ whenever $p_{X}(x)>0$;
$p_{A|X,\overline{M}_{k}}(a|x,\overline{m}_{k})>\epsilon>0$ whenever
$p_{X,\overline{M}_{k}}(x,\overline{m}_{k})>0$, $\forall k\in[K]$.

Before giving the identification results, we introduce the following
notational shorthands:

$\displaystyle\overline{M}_{k}(\overline{a}_{k})$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}\big{(}\overline{M}_{k-1}(\overline{a}_{k-1}),M_{k}(a_{k},\overline{M}_{k-1}(\overline{a}_{k-1}))\big{)},\forall
k\in[K],$ $\displaystyle\psi_{\overline{a}}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}\mathbb{E}[Y(a_{K+1},\overline{M}_{k}(\overline{a}_{k}))],$

where $\overline{M}_{k}(\overline{a}_{k})$ is defined iteratively, with the
assumption that $\overline{M}_{0}(\overline{a}_{0})=\varnothing$. For example,
when $K=3$,

$\psi_{\overline{a}}=\mathbb{E}\big{[}Y\big{(}a_{4},M_{1}(a_{1}),M_{2}(a_{2},M_{1}(a_{1})),M_{3}(a_{3},M_{1}(a_{1}),M_{2}(a_{2},M_{1}(a_{1})))\big{)}\big{]}.$

Theorem 1 states that $\psi_{\overline{a}}$ is identified under Assumptions
1*-3*.

######  Theorem 1.

Under Assumptions 1*-3*, we have

$\psi_{\overline{a}}=\int_{x}\int_{\overline{m}_{K}}\mathbb{E}[Y|x,a_{K+1},\overline{m}_{K}]\big{[}\prod\limits_{k=1}^{K}dP(m_{k}|x,a_{k},\overline{m}_{k-1})\big{]}dP(x).$
(2)

The above equation extends Pearl’s (2001) and Daniel _et al._ ’s (2015)
mediation formula to the case of $K$ causally ordered mediators. Following the
terminology of Tchetgen Tchetgen and Shpitser (2012), we refer to the right-
hand side of equation (1) as the generalized mediation functional (GMF).
Theorem 1 echoes Avin _et al._ ’s (2005) recanting witness criterion: a
potential outcome is identified (in expectation) if the value that a mediator
$M_{k}$ takes, i.e., $M_{k}(a_{k})$, is carried over to all future mediators.
This result leads us to focus on the set of expected potential outcomes and
PSEs that are nonparametrically identified. For example, to assess the
mediating role of $M_{k}$, we focus on the composite causal path $A\to
M_{k}\rightsquigarrow Y$, where, as before, the squiggle arrow encompasses all
possible causal paths from $M_{k}$ to $Y$. An identifiable PSE for this path
can be expressed as

$\tau_{A\to M_{k}\rightsquigarrow
Y}(\overline{a}_{k-1},\cdot,\underline{a}_{k+1})=\psi_{\overline{a}_{k-1},1,\underline{a}_{k+1}}-\psi_{\overline{a}_{k-1},0,\underline{a}_{k+1}},$

where $\underline{a}_{k+1}\stackrel{{\scriptstyle\Delta}}{{=}}(a_{k+1},\ldots
a_{K+1})$. The notation $\psi_{\overline{a}}$ makes it clear that the average
total effect (ATE) of $A$ on $Y$ can be decomposed into $K+1$ identifiable
PSEs corresponding to $A\to Y$ and $A\to M_{k}\rightsquigarrow Y$ ($k\in[K]$):

$\textup{ATE}=\psi_{\overline{1}}-\psi_{\overline{0}}=\underbrace{\psi_{\overline{0}_{K},1}-\psi_{\overline{0}_{K+1}}}_{A\to
Y}+\sum_{k=1}^{K}\underbrace{\big{(}\psi_{\overline{0}_{k-1},\underline{1}_{k}}-\psi_{\overline{0}_{k},\underline{1}_{k+1}}\big{)}}_{A\to
M_{k}\rightsquigarrow Y}.$ (3)

To be sure, equation (3) is not the only way of decomposing the ATE. Depending
on the order in which the paths $A\to Y$ and $A\to M_{k}\rightsquigarrow Y$
($k\in[K]$) are considered, there are $(K+1)!$ different ways of decomposing
the ATE. In the above decomposition,
$\psi_{\overline{0}_{K},1}-\psi_{\overline{0}_{K+1}}$ corresponds to the NDE
if the mediators $\overline{M}_{K}$ are considered as a whole.

###  3 Estimation

In this section, we focus on the estimation of the GMF, i.e., the right-hand
side of equation (2). When Assumptions 1*-3* hold, the GMF is equal to the
causal parameter $\psi_{\overline{a}}$, but otherwise, it is still a well-
defined statistical parameter of potential scientific interest. To distinguish
it from the causal parameter $\psi_{\overline{a}}$, we henceforth denote the
GMF by $\overline{\theta}_{a}$.

#### 3.1 MLE, Regression-Imputation, and Weighting

Equation (2) suggests that $\theta_{\overline{a}}$ can be estimated via
maximum likelihood (MLE) (Miles _et al._ 2017). Specifically, we can fit a
parametric model for each $p(m_{k}|x,a_{k},\overline{m}_{k-1})$ ($k\in[K]$)
and for $\mathbb{E}[Y|x,a_{K+1},\overline{m}_{K}]$, and then estimate the GMF
via the following equation:

$\hat{\theta}_{\overline{a}}^{\textup{mle}}=\mathbb{P}_{n}\big{[}\int_{\overline{m}_{K}}\hat{\mathbb{E}}[Y|X,a_{K+1},\overline{m}_{K}]\big{(}\prod\limits_{k=1}^{K}\hat{p}(m_{k}|x,a_{k},\overline{m}_{k-1})d\nu(m_{k})\big{)}\big{]},$
(4)

where $\mathbb{P}_{n}[\cdot]=n^{-1}\sum_{i}[\cdot]_{i}$ and $\nu(\cdot)$ is an
appropriate dominating measure. This approach works best when the mediators
$M_{1},M_{2},\ldots M_{K}$ are all discrete and the covariates $X$ are low-
dimensional, in which case the working models for
$p(m_{k}|x,a_{k},\overline{m}_{k-1})$ are simply models for the conditional
probabilities of $M_{k}$ that can be reliably estimated. When some of the
mediators are continuous/multivariate or when the covariates $X$ are high-
dimensional, estimates of the corresponding conditional density/probability
functions can be unstable and sensitive to model misspecification. This
problem could be mitigated by imposing highly constrained functional forms on
the conditional means of the mediators and the outcome. For example, when
$\mathbb{E}[M_{k}|x,a_{k},\overline{m}_{k-1}]$ and
$\mathbb{E}[Y|x,a_{K+1},\overline{m}_{K}]$ are all assumed to be linear with
no higher-order or interaction terms,
$\hat{\theta}_{\overline{a}}^{\textup{mle}}$ will reduce to a simple function
of regression coefficients (e.g., Alwin and Hauser 1975). Yet, the assumptions
of linearity and additivity are unrealistic in many applications, which may
lead to biased estimates of $\theta_{\overline{a}}$. Below, we describe
several imputation- and weighting-based strategies for estimating
$\theta_{\overline{a}}$.

First, we observe that the GMF can be written as

$\theta_{\overline{a}}=\mathbb{E}_{X}\big{[}\underbrace{\mathbb{E}_{M_{1}|X,a_{1}}\ldots\underbrace{\mathbb{E}_{M_{K}|X,a_{K},\overline{M}_{K-1}}\underbrace{\mathbb{E}[Y|X,a_{K+1},\overline{M}_{K}]}_{\stackrel{{\scriptstyle\Delta}}{{=}}\mu_{K}(X,\overline{M}_{K})}}_{\stackrel{{\scriptstyle\Delta}}{{=}}\mu_{K-1}(X,\overline{M}_{K-1})}}_{\stackrel{{\scriptstyle\Delta}}{{=}}\mu_{0}(X)}\big{]}.$
(5)

This expression suggests that $\theta_{\overline{a}}$ can be estimated via an
iterated regression-imputation (RI) approach (Zhou and Yamamoto 2020):

  1. 1.

Estimate $\mu_{K}(X,\overline{M}_{K})$ by fitting a parametric model for the
conditional mean of $Y$ given $(X,A,\overline{M}_{K})$ and then setting
$A=a_{K+1}$ for all units;

  2. 2.

For $k=K-1,\ldots 0$, estimate $\mu_{k}(X,\overline{M}_{k})$ by fitting a
parametric model for the conditional mean of $\mu_{k+1}(X,\overline{M}_{k+1})$
and then setting $A=a_{k+1}$ for all units;

  3. 3.

Estimate $\theta_{\overline{a}}$ by averaging the fitted values
$\hat{\mu}_{0}(X)$ among all units:

$\hat{\theta}_{\overline{a}}^{\textup{ri}}=\mathbb{P}_{n}\big{[}\hat{\mu}_{0}(X)\big{]}.$
(6)

The regression-imputation estimator can be seen as an extension of the
imputation strategy proposed by Vansteelandt _et al._ (2012) for estimating
the NDE and NIE in the one-mediator setting. Since this approach requires
modeling only the conditional means of observed/imputed outcomes given
different sets of mediators, it is more flexible to use with
continuous/multivariate mediators than MLE. Nonetheless, because
$\mu_{k}(x,\overline{m}_{k})$ is estimated iteratively, correct specification
of all of the outcome models is required for
$\hat{\theta}_{\overline{a}}^{\textup{ri}}$ to be consistent. Thus, in
practice, when parametric models are used to estimate
$\mu_{k}(x,\overline{m}_{k})$, care should be taken to ensure that the outcome
models used to estimate these functions are mutually compatible. For example,
if $\mu_{1}(X,M_{1})$ follows a linear model that includes $X$ and $X^{2}$ as
predictors, then the model used to estimate
$\mu_{0}(X)=\mathbb{E}[\mu_{1}(X,M_{1})|X,A=a_{1}]$ should also include $X$
and $X^{2}$ in the predictor set.

The GMF can also be written as

$\displaystyle\theta_{\overline{a}}$
$\displaystyle=\mathbb{E}\big{[}\frac{\mathbb{I}(A=a_{K+1})}{p(a_{K+1}|X)}\big{(}\prod\limits_{k=1}^{K}\frac{p(M_{k}|X,a_{k},\overline{M}_{k-1})}{p(M_{k}|X,a_{K+1},\overline{M}_{k-1})}\big{)}Y\big{]}.$

This expression suggests a weighting estimator of $\theta_{\overline{a}}$:

$\hat{\theta}_{\overline{a}}^{\textup{w-m}}=\mathbb{P}_{n}\big{[}\frac{\mathbb{I}(A=a_{K+1})}{\hat{p}(a_{K+1}|X)}\big{(}\prod\limits_{k=1}^{K}\frac{\hat{p}(M_{k}|X,a_{k},\overline{M}_{k-1})}{\hat{p}(M_{k}|X,a_{K+1},\overline{M}_{k-1})}\big{)}Y\big{]}.$
(7)

This estimator can be seen as an extension of the weighting estimator proposed
in VanderWeele _et al._ (2014) for the case of two mediators. It shares a
limitation of $\hat{\theta}_{\overline{a}}^{\textup{mle}}$ in that it requires
estimates of the conditional densities/probabilities of the mediators, which
tend to be noisy if the mediators are continuous or multivariate. This
problem, however, can be sidestepped by recasting the mediator density ratios,
via Bayes’ rule, as odds ratios in terms of the treatment variable:

$\displaystyle\frac{p(M_{k}|X,a_{k},\overline{M}_{k-1})}{p(M_{k}|X,a_{K+1},\overline{M}_{k-1})}$
$\displaystyle=\frac{p(a_{k}|X,\overline{M}_{k})/p(a_{K+1}|X,\overline{M}_{k})}{p(a_{k}|X,\overline{M}_{k-1})/p(a_{K+1}|X,\overline{M}_{k-1})}.$

This observation leads to an alternative weighting estimator based on
estimates of the conditional probabilities of treatment given different sets
of mediators:

$\hat{\theta}_{\overline{a}}^{\textup{w-a}}=\mathbb{P}_{n}\big{[}\frac{\mathbb{I}(A=a_{K+1})}{\hat{p}(a_{1}|X)}\big{(}\prod\limits_{k=1}^{K}\frac{\hat{p}(a_{k}|X,\overline{M}_{k})}{\hat{p}(a_{k+1}|X,\overline{M}_{k})}\big{)}Y\big{]}.$
(8)

In applications where the mediators are continuous/multivariate,
$\hat{\theta}_{\overline{a}}^{\textup{w-a}}$ should be easier to work with
than $\hat{\theta}_{\overline{a}}^{\textup{w-m}}$. Yet, the parameters for
$p(a|x,\overline{m}_{k})$ are not variationally independent across different
values of $k$. As in the case of the regression-imputation estimator, care
should be taken to ensure the compatibility of the models specified for
$p(a|x,\overline{m}_{k})$ (see Miles _et al._ 2020 for some practical
recommendations).

The regression-imputation approach and the weighting approach can be combined
to form various “hybrid estimators” of $\theta_{\overline{a}}$. For example,
in the case of $K=2$, one can use regression-imputation to estimate
$\mu_{2}(x,m_{1},m_{2})$, another regression-imputation step to estimate
$\mu_{1}(x,m_{1})$, and weighting to estimate $\theta_{\overline{a}}$,
yielding an “RI-RI-W” estimator:

$\hat{\theta}_{\overline{a}}^{\textup{ri-
ri-w}}=\mathbb{P}_{n}\big{[}\frac{\mathbb{I}(A=a_{1})}{\hat{p}(a_{1}|X)}\hat{\mu}_{1}(X,M_{1})\big{]}.$
(9)

One can also use regression-imputation to estimate $\mu_{2}(x,m_{1},m_{2})$
and then employ appropriate weights to estimate $\theta_{\overline{a}}$, which
leads to an “RI-W-W” estimator:

$\hat{\theta}_{\overline{a}}^{\textup{ri-
w-w}}=\mathbb{P}_{n}\big{[}\frac{\mathbb{I}(A=a_{2})}{\hat{p}(a_{2}|X)}\frac{\hat{p}(M_{1}|X,a_{1})}{\hat{p}(M_{1}|X,a_{2})}\hat{\mu}_{2}(X,M_{1},M_{2})\big{]}.$
(10)

In fact, with $K$ mediators, there are $2^{K+1}$ different ways to combine
regression-imputation and weighting, each of which involves estimating $K+1$
nuisance functions, which entail a choice between $p(a|x)$ and $\mu_{0}(x)$
and a choice between $p(m_{k}|x,a,\overline{m}_{k-1})$ and
$\mu_{k}(x,\overline{m}_{k})$ for each $k\in[K]$ (see Supplementary Material B
for detailed expressions of these hybrid estimators in the case of $K=2$). As
with $\hat{\theta}_{\overline{a}}^{\textup{mle}}$,
$\hat{\theta}_{\overline{a}}^{\textup{ri}}$,
$\hat{\theta}_{\overline{a}}^{\textup{w-m}}$, and
$\hat{\theta}_{\overline{a}}^{\textup{w-a}}$, each of these hybrid estimators
will be consistent only if the corresponding nuisance functions are all
correctly specified and consistently estimated. In applications where the
pretreatment covariates $X$ and/or the mediators have many components, all of
the above estimators will be prone to model misspecification bias.

#### 3.2 Multiply Robust and Semiparametric Efficient Estimation

Henceforth, let $O=(X,A,\overline{M}_{K},Y)$ denote the observed data, and
$\mathcal{P}_{\textup{np}}$ a nonparametric model over $O$ wherein all laws
$P$ satisfy the positivity assumption described in Section 2.2. In addition,
define $\mu_{k}(X,\overline{M}_{k})$ iteratively as in equation (5):

$\displaystyle\mu_{K}(X,\overline{M}_{K})$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}\mathbb{E}[Y|X,a_{K+1},\overline{M}_{K}]$
$\displaystyle\mu_{k}(X,\overline{M}_{k})$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}\mathbb{E}[\mu_{k+1}(X,\overline{M}_{k+1})|X,a_{k+1},\overline{M}_{k}],\quad
k=K-1,\ldots,0.$

######  Theorem 2.

The efficient influence function (EIF) of $\theta_{\overline{a}}$ in
$\mathcal{P}_{\textup{np}}$ is given by

$\varphi_{\overline{a}}(O)=\sum_{k=0}^{K+1}\varphi_{k}(O),$ (11)

where

$\displaystyle\varphi_{0}(O)$
$\displaystyle=\mu_{0}(X)-\theta_{\overline{a}},$
$\displaystyle\varphi_{k}(O)$
$\displaystyle=\frac{\mathbb{I}(A=a_{k})}{p(a_{k}|X)}\Big{(}\prod_{j=1}^{k-1}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{k},\overline{M}_{j-1})}\Big{)}\big{(}\mu_{k}(X,\overline{M}_{k})-\mu_{k-1}(X,\overline{M}_{k-1})\big{)},\quad
k\in[K],$ $\displaystyle\varphi_{K+1}(O)$
$\displaystyle=\frac{\mathbb{I}(A=a_{K+1})}{p(a_{K+1}|X)}\Big{(}\prod_{j=1}^{K}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{K+1},\overline{M}_{j-1})}\Big{)}\big{(}Y-\mu_{K}(X,\overline{M}_{K})\big{)}.$

The semiparametric efficiency bound for any regular and asymptotically linear
estimator of $\theta_{\overline{a}}$ in $\mathcal{P}_{\textup{np}}$ is
therefore
$\mathbb{E}\big{[}\big{(}\varphi_{\overline{a}}(O)\big{)}^{2}\big{]}$.

We now present two estimators of $\theta_{\overline{a}}$ based on the EIF.
First, consider the factorized likelihood of $O$:
$p(O)=p(X)p(A|X)\Big{(}\prod_{k=1}^{K}p(M_{k}|X,A,\overline{M}_{k-1})\big{)}p(Y|X,A,\overline{M}_{K}).$
Suppose we have estimated $K+2$ nuisance functions, each of which corresponds
to a component of $p(O)$: $\hat{\pi}_{0}(a|x)$ for $p(a|x)$,
$\hat{f}_{k}(m_{k}|x,a,\overline{m}_{k-1})$ for
$p(m_{k}|x,a,\overline{m}_{k-1})$, and $\hat{\mu}_{K}(x,\overline{m}_{K})$ for
$\mathbb{E}[Y|x,a_{K+1},\overline{m}_{K}]$. The GMF can now be estimated as

$\displaystyle\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}=$
$\displaystyle\mathbb{P}_{n}\big{[}\frac{\mathbb{I}(A=a_{K+1})}{\hat{\pi}_{0}(a_{K+1}|X)}\Big{(}\prod_{j=1}^{K}\frac{\hat{f}_{j}(M_{j}|X,a_{j},\overline{M}_{j-1})}{\hat{f}_{j}(M_{j}|X,a_{K+1},\overline{M}_{j-1})}\Big{)}\big{(}Y-\hat{\mu}_{K}(X,\overline{M}_{K})\big{)}$
$\displaystyle+\sum_{k=1}^{K}\frac{\mathbb{I}(A=a_{k})}{\hat{\pi}_{0}(a_{k}|X)}\Big{(}\prod_{j=1}^{k-1}\frac{\hat{f}_{j}(M_{j}|X,a_{j},\overline{M}_{j-1})}{\hat{f}_{j}(M_{j}|X,a_{k},\overline{M}_{j-1})}\Big{)}\big{(}\hat{\mu}_{k}^{\textup{mle}}(X,\overline{M}_{k})-\hat{\mu}_{k-1}^{\textup{mle}}(X,\overline{M}_{k-1})\big{)}$
$\displaystyle+\hat{\mu}_{0}^{\textup{mle}}(X)\big{]},$ (12)

where
$\hat{\mu}_{K}^{\textup{mle}}(X,\overline{M}_{K})=\hat{\mu}_{K}(X,\overline{M}_{K})$
and $\hat{\mu}_{k}^{\textup{mle}}(X,\overline{M}_{k})$ is iteratively
constructed as

$\hat{\mu}_{k}^{\textup{mle}}(X,\overline{M}_{k})=\int\hat{\mu}_{k+1}^{\textup{mle}}(X,\overline{M}_{k},m_{k+1})\hat{f}_{k+1}(m_{k+1}|X,a_{k+1},\overline{M}_{k})d\nu(m_{k+1}),\quad
k=K-1,\ldots 0.$ (13)

When $M_{k+1}$ involves continuous components, equation (13) can be evaluated
via Monte Carlo simulation.

When some of the mediators are continuous/multivariate, it can be difficult to
estimate the conditional distributions $p(m_{k}|x,a,\overline{m}_{k-1})$. In
such cases, it is often preferable to estimate the mediator density ratios
using the corresponding odds ratios of the treatment variable, and estimate
the functions $\mu_{k}(x,\overline{m}_{k})$ using the regression-imputation
approach. Specifically, suppose we have estimated $2(K+1)$ nuisance functions:
$\hat{\pi}_{0}(a|x)$ for $p(a|x)$, $\hat{\pi}_{k}(a|x,\overline{m}_{k})$ for
$p(a|x,\overline{m}_{k})$ ($k\in[K]$), and $\hat{\mu}_{k}(x,\overline{m}_{k})$
for $\mu_{k}(x,\overline{m}_{k})$ ($k\in\\{0,1,\ldots K\\}$), where for $k<K$,
$\mu_{k}(x,\overline{m}_{k})$ is estimated iteratively by fitting a model for
the conditional mean of $\hat{\mu}_{k+1}(X,\overline{M}_{k+1})$ given
$(X,A,\overline{M}_{k})$ and then setting $A=a_{k+1}$ for all units. The GMF
can then be estimated as

$\displaystyle\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}=$
$\displaystyle\mathbb{P}_{n}\big{[}\frac{\mathbb{I}(A=a_{K+1})}{\hat{\pi}_{0}(a_{1}|X)}\Big{(}\prod_{j=1}^{K}\frac{\hat{\pi}_{j}(a_{j}|X,\overline{M}_{j})}{\hat{\pi}_{j}(a_{j+1}|X,\overline{M}_{j})}\Big{)}\big{(}Y-\hat{\mu}_{K}(X,\overline{M}_{K})\big{)}$
$\displaystyle+\sum_{k=1}^{K}\frac{\mathbb{I}(A=a_{k})}{\hat{\pi}_{0}(a_{1}|X)}\Big{(}\prod_{j=1}^{k-1}\frac{\hat{\pi}_{j}(a_{j}|X,\overline{M}_{j})}{\hat{\pi}_{j}(a_{j+1}|X,\overline{M}_{j})}\Big{)}\big{(}\hat{\mu}_{k}(X,\overline{M}_{k})-\hat{\mu}_{k-1}(X,\overline{M}_{k-1})\big{)}$
$\displaystyle+\hat{\mu}_{0}(X)\big{]}.$ (14)

The multiple robustness and semiparametric efficiency of
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ and
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ are given below.

######  Theorem 3.

Let $\eta_{1}=\\{\pi_{0},f_{1},\ldots f_{K},\mu_{K}\\}$ denote the $K+2$
nuisance functions involved in
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$, and
$\eta_{2}=\\{\pi_{0},\ldots\pi_{K},\mu_{0},\ldots\mu_{K}\\}$ denote the
$2(K+1)$ nuisance functions involved in
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$. Suppose that Assumption 3*
(positivity) and suitable regularity conditions for estimating equations
(e.g., Newey and McFadden 1994, p. 2148) hold. In addition, suppose that
$\mu_{K}(x,\overline{m}_{K})$ is bounded over the support of
$(X,\overline{M}_{K})$. Then, when the elements of $\eta_{1}$ and $\eta_{2}$
are estimated via parametric models,

  1. 1.

$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ is consistent and
asymptotically normal (CAN) if $K+1$ of the $K+2$ nuisance functions in
$\eta_{1}$ are correctly specified and their parameter estimates are
$\sqrt{n}$-consistent; it is semiparametric efficient if all of the $K+2$
nuisance functions in $\mathcal{\eta}_{1}$ are correctly specified and their
parameter estimates are $\sqrt{n}$-consistent.

  2. 2.

$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ is CAN if $\exists
k\in\\{0,\ldots K+1\\}$, the first $k$ treatment models
$\pi_{0},\ldots\pi_{k-1}$ and the last $K+1-k$ outcome models
$\mu_{k},\ldots\mu_{K}$ in $\eta_{2}$ are correctly specified and their
parameter estimates are $\sqrt{n}$-consistent; it is semiparametric efficient
if all of the treatment and outcome models in $\eta_{2}$ are correctly
specified and their parameter estimates are $\sqrt{n}$-consistent.

Both $\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ and
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ are $K+2$-robust in the sense
that they are CAN provided that one of $K+2$ sets of nuisance functions is
correctly specified and the corresponding parameter estimates are
$\sqrt{n}$-consistent. Several special cases are worth noting. First, in the
degenerate case where $K=0$, it is clear that both
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ and
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ reduce to the standard doubly
robust estimator for $\mathbb{E}[Y(a)]$ (Scharfstein _et al._ 1999). Second,
when $K=1$, $\hat{\theta}_{01}^{\textup{eif}_{1}}$ coincides with Tchetgen
Tchetgen and Shpitser’s (2012) triply robust estimator for
$\mathbb{E}[Y\big{(}1,M(0)\big{)}]$. Finally, when $K=2$,
$\hat{\theta}_{010}^{\textup{eif}_{1}}$ is identical to Miles _et al._ ’s
(2020) estimator for $\theta_{010}$. For this case, however, Miles _et al._
provide a slightly weaker condition than that implied by Theorem 3 for
$\hat{\theta}_{010}^{\textup{eif}_{1}}$ to be CAN. Specifically, they showed
that $\hat{\theta}_{010}^{\textup{eif}_{1}}$ remains CAN even if both $f_{1}$
and $\mu_{2}$ are misspecified. In Section 4, we show that the conditions for
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ to be CAN can also be relaxed
for several particular types of PSEs, including the natural path-specific
effect (nPSE), of which $\psi_{010}-\psi_{000}$ is a special case. For the
$K=2$ case, Miles _et al._ also noted that the mediator density ratios in
$\hat{\theta}_{010}^{\textup{eif}_{1}}$ can be indirectly estimated through
models for $\pi_{1}$ and $\pi_{2}$. Clearly, this approach will result in
$\hat{\theta}_{010}^{\textup{eif}_{2}}$ if the $\mu_{k}(x,\overline{m}_{k})$
functions are in the meanwhile estimated through regression-imputation. The
$K+2$-robustness of $\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ and
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$, interestingly, resembles the
multiple robustness of the Bang-Robins (2005) estimator for the mean of a
potential outcome with time-varying treatments and time-varying confounders
(Luedtke _et al._ 2017; Molina _et al._ 2017; Rotnitzky _et al._ 2017).

To gain some intuition as to why
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ is $K+2$-robust, consider
cases in which only one nuisance function in $\eta_{1}$ is misspecified. When
only $\pi_{0}$ is misspecified, all terms inside $\mathbb{P}_{n}[\cdot]$ but
$\hat{\mu}_{0}^{\textup{mle}}(X)$ will have a zero mean (asymptotically),
leaving only $\mathbb{P}_{n}[\hat{\mu}_{0}^{\textup{mle}}(X)]$ (i.e., the MLE
estimator (4)), which is consistent because the corresponding nuisance
functions $\\{f_{1},\ldots f_{K},\mu_{K}\\}$ are all correctly specified. When
only $\mu_{K}$ is misspecified, all terms involving
$\hat{\mu}_{K}(X,\overline{M}_{K})$ and
$\hat{\mu}_{k}^{\textup{mle}}(X,\overline{M}_{k})$ ($k=0,1,\ldots K-1$) inside
$\mathbb{P}_{n}[\cdot]$ will have a zero mean (asymptotically), leaving only a
weighted average of $Y$ (i.e., the weighting estimator (7)), which is
consistent because the corresponding nuisance functions
$\\{\pi_{0},f_{1},\ldots f_{K}\\}$ are all correctly specified. Finally, when
only $f_{k^{\prime}}$ is misspecified (for some $k^{\prime}\in[K]$), it can be
shown that all terms involving $\hat{f}_{k^{\prime}}$ and
$\hat{\mu}_{k}^{\textup{mle}}(X,\overline{M}_{k})$ ($\forall k<k^{\prime}$)
inside $\mathbb{P}_{n}[\cdot]$ will have a zero mean (asymptotically), leaving
only a weighted average of
$\hat{\mu}_{k^{\prime}}^{\textup{mle}}(X,\overline{M}_{k^{\prime}})$. The
latter constitutes a “hybrid” estimator similar to those mentioned in the
previous section, and it is consistent in this case because its nuisance
functions $\\{\pi_{0},f_{1},\ldots f_{k^{\prime}-1},f_{k^{\prime}+1},\ldots
f_{K},\mu_{K}\\}$ are all correctly specified.

The $K+2$-robustness of $\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ is
due to a similar logic to that of
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$. Yet, different from
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$,
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ involves estimating $2(K+1)$
nuisance functions, $K+1$ for the conditional probabilities of treatment and
$K+1$ for the conditional means of observed/imputed outcomes. Also, unlike
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$, the treatment models
involved in $\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ are not
variationally independent; neither are the outcome models. For example, when
$M_{K}\perp\\!\\!\\!\perp A|X,\overline{M}_{K-1}$,
$\pi_{K}(A|X,\overline{M}_{K})$ should be identical to
$\pi_{K-1}(A|X,\overline{M}_{K-1})$; similarly, when $M_{K}\perp\\!\\!\\!\perp
Y|X,A,\overline{M}_{K-1}$, $\mu_{K}(X,\overline{M}_{K})$ should be identical
to $\mu_{K-1}(X,\overline{M}_{K-1})$. Thus, in practice, both the treatment
and outcome models should be specified in a mutually compatible way, otherwise
some of the conditions in Theorem 3 may fail by design.

The local efficiency of $\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ and
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ is due to the fact that both
of the EIF-based estimating equations (12) and (14) have a zero derivative
with respect to the nuisance functions at the truth. This property, referred
to as “Neyman orthogonality” by Chernozhukov _et al._ (2018), implies that
first step estimation of the nuisance functions has no first order effect on
the influence functions of $\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$
and $\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$. This property suggests
that the nuisance functions can be estimated using data-adaptive/machine
learning methods or their ensembles. In this case, these estimators will still
be consistent as long as the nuisance functions associated with one of the
$K+2$ conditions in Theorem 3 are consistently estimated. For
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$, an added advantage of
employing data-adaptive methods to estimate the nuisance functions is that, by
exploring a larger space within $\mathcal{P}_{\textup{np}}$, the risk of model
incompatibility is reduced.

When data-adaptive/machine learning methods are used to estimate the nuisance
functions, it is advisable to use sample splitting to render the empirical
process term asymptotically negligible (Zheng and van der Laan 2011;
Chernozhukov _et al._ 2018; Newey and Robins 2018). For example, Chernozhukov
_et al._ (2018) suggest the method of “cross-fitting,” which involves the
following steps: (a) randomly partition the sample $S$ into $J$ folds:
$S_{1},S_{2}\ldots S_{J}$; (b) for each $j$, obtain a fold-specific estimate
of the target parameter using only data from $S_{j}$ (“main sample”), but with
nuisance functions learned from the remainder of the sample (i.e.,
$S\backslash S_{j}$; “auxiliary sample”); (c) average these fold-specific
estimates to form a final estimate of the target parameter.

When cross-fitting is used, $\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$
and $\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ will be semiparametric
efficient if the corresponding nuisance function estimates are all consistent
and converge at sufficiently fast rates. For example, a sufficient (but not
necessary) condition for $\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ and
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ to attain the semiparametric
efficiency bound is when all of the nuisance function estimates converge at
faster-than-$n^{-1/4}$ rates. More precise conditions are given in Theorem 4.

######  Theorem 4.

Let
$\hat{\eta}_{1}=\\{\hat{\pi}_{0},\hat{f}_{1},\ldots\hat{f}_{K},\hat{\mu}_{K}\\}$
and
$\hat{\eta}_{2}=\\{\hat{\pi}_{0},\ldots\hat{\pi}_{K},\hat{\mu}_{0},\ldots\hat{\mu}_{K}\\}$
denote estimates of the nuisance functions involved in
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ and
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$, respectively. Let
$r_{n}(\cdot)$ denote a mapping from a nuisance function estimator to its
$L_{2}(P)$ convergence rate where $P$ represents the true distribution of
$O=(X,A,\overline{M}_{K},Y)$. Suppose that Assumption 3* (positivity) holds
for both the true distribution $P$ and its estimates implied by
$\hat{\eta}_{1}$ and $\hat{\eta}_{2}$, and that all other assumptions required
for Theorem 3 hold. Then, when the nuisance functions are estimated via data-
adaptive methods and cross-fitting,

  1. 1.

$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ is consistent if $K+1$ of the
$K+2$ elements in $\hat{\eta}_{1}$ are consistent in the $L_{2}$-norm; it is
CAN and semiparametric efficient if all elements in $\hat{\eta}_{1}$ are
consistent in the $L_{2}$-norm and $\sum\limits_{u,v\in\hat{\eta}_{1};u\neq
v}r_{n}(u)r_{n}(v)=o(n^{-1/2})$;

  2. 2.

$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ is consistent if $\exists
k\in\\{0,\ldots K+1\\}$,
$\hat{\pi}_{0},\ldots\hat{\pi}_{k-1},\hat{\mu}_{k},\ldots\hat{\mu}_{K}$ are
all consistent in the $L_{2}$-norm; it is CAN and semiparametric efficient if
all elements in $\hat{\eta}_{2}$ are consistent in the $L_{2}$-norm and
$\sum\limits_{j=0}^{K}r_{n}(\hat{\pi}_{j})r_{n}(\hat{\mu}_{j})=o(n^{-1/2})$.

The multiple robustness result for
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ echoes Theorem 3. Moreover,
the first part of Theorem 4 states that
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ is CAN and semiparametric
efficient if all nuisance functions in $\eta_{1}$ are consistently estimated
and, for every two nuisance functions in $\eta_{1}$, the product of their
convergence rates is $o(n^{-1/2})$. Thus
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ is CAN and semiparametric
efficient if all of the $K+2$ nuisance function estimates are consistent and
converge at faster-than-$n^{-1/4}$ rates, but it will also attain
semiparametric efficiency under alternative conditions. For example, when
estimates of the treatment and mediator models $\\{\hat{\pi}_{0}$,
$\hat{f}_{1},\ldots\hat{f}_{K}\\}$ all converge to the truth at a rate of
$n^{-1/3}$ and estimates of the outcome model $\hat{\mu}_{K}$ converge to the
truth at a rate of $n^{-1/5}$, the product of the convergence rates of any two
elements in $\hat{\eta}_{1}$ is either $O(n^{-1/3})O(n^{-1/3})=O(n^{-2/3})$ or
$O(n^{-1/3})O(n^{-1/5})=O(n^{-8/15})$, both faster than $O(n^{-1/2})$.

The second part of Theorem 4 states that
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ is consistent if there exists
a $k$ such that the first $k$ treatment models and the last $K+1-k$ outcome
models in $\eta_{2}$ are consistently estimated, echoing Theorem 3. As with
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$,
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ will be CAN and
semiparametric efficient if all of the required nuisance functions are
consistently estimated and converge at faster-than-$n^{-1/4}$ rates. The rate
condition
$\sum\limits_{j=0}^{K}r_{n}(\hat{\pi}_{j})r_{n}(\hat{\mu}_{j})=o(n^{-1/2})$
appears to be weaker than that for
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ as it involves the sum of
only $K+1$, rather than ${K+2\choose 2}$, product terms. Because the outcome
models are estimated iteratively, the convergence rate of $\hat{\mu}_{k}$ will
in general depend on the convergence rates of
$\\{\hat{\mu}_{k+1},\ldots\hat{\mu}_{K}\\}$. That is, if
$r_{n}(\hat{\mu}_{k+1})=O(n^{\delta})$, $r_{n}(\hat{\mu}_{k})$ is unlikely to
be faster than $O(n^{\delta})$. Nonetheless,
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ will be CAN and
semiparametric efficient under relatively weak conditions — for example, when
estimates of the treatment models all converge to the truth at a rate of
$n^{-1/3}$ and estimates of the outcome models all converge to the truth at a
rate of $n^{-1/5}$, in which case
$\sum\limits_{j=0}^{K}r_{n}(\hat{\pi}_{j})r_{n}(\hat{\mu}_{j})=\sum\limits_{j=0}^{K}O(n^{-8/15})=o(n^{-1/2})$.

For inference on $\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ and
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$, a simple variance estimator
can be constructed from the empirical analog of the EIF, i.e.,
$\mathbb{P}_{n}[\hat{\varphi}_{\overline{a}}^{2}(O)]/n$. However, unlike
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ and
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$, this variance estimator is
not multiply robust — it will be consistent only if the conditions for
semiparametric efficiency in Theorem 3 or Theorem 4 are satisfied. Thus, when
the nuisance functions are estimated using parametric models, the variance
estimator constructed from the empirical EIF may be inconsistent even when the
corresponding estimator for $\theta_{\overline{a}}$ is CAN — for example, when
only $K+1$ of the $K+2$ nuisance functions involved in
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ are correctly specified. In
this case, the nonparametric bootstrap is a convenient approach to more robust
inference. When the nuisance functions are estimated using data-
adaptive/machine learning methods, however, the nonparametric bootstrap is not
theoretically justified, and the EIF-based variance estimator may still be
preferred.

#### 3.3 Multiply Robust Regression-Imputation Estimators

Both of the multiply robust estimators described above involve inverse
probability weights. When the positivity assumption is nearly violated, the
inverse probability weights tend to be highly variable, which may lead to poor
finite sample performance (Kang and Schafer 2007; Petersen _et al._ 2012). A
variety of methods have been proposed to reduce the influence of highly
variable weights on doubly robust and multiply robust estimators in similar
settings (e.g., Robins _et al._ 2007; Tchetgen Tchetgen and Shpitser 2012;
Seaman and Vansteelandt 2018). Among them, a common strategy is to tailor the
estimating equation of the outcome model(s) such that the terms involving
inverse probability weights will equal zero, leaving only a regression-
imputation or “substitution” estimator that typically resides in the parameter
space of the estimand. Below, we briefly describe how this approach can be
adapted to $\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ and
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$.

Let us start with $\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$, which can
be written as

$\displaystyle\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}=$
$\displaystyle\mathbb{P}_{n}\big{[}\hat{w}_{K}(X,A,\overline{M}_{K})\big{(}Y-\hat{\mu}_{K}(X,\overline{M}_{K})\big{)}$
$\displaystyle+\sum_{k=1}^{K}\hat{w}_{k-1}(X,A,\overline{M}_{k-1})\big{(}\hat{\mu}_{k}(X,\overline{M}_{k})-\hat{\mu}_{k-1}(X,\overline{M}_{k-1})\big{)}$
$\displaystyle+\hat{\mu}_{0}(X)\big{]},$ (15)

where $\hat{w}_{k}(A,X,\overline{M}_{k})$ ($0\leq k\leq K$) are estimates of
the corresponding inverse probability weights as displayed in equation (14).
Note that the nuisance functions $\hat{\mu}_{k}(X,\overline{M}_{k})$ ($0\leq
k\leq K$) here are all estimated via the regression-imputation approach. When
the corresponding outcome models are fitted via generalized linear models
(GLM) with canonical links, one can either (a) fit weighted GLMs (with an
intercept term) for $\hat{\mu}_{k}(X,\overline{M}_{k})$ using
$\hat{w}_{k}(A,X,\overline{M}_{k})$ as weights, or (b) add the corresponding
inverse probability weight as an additional covariate in these regressions
(Robins _et al._ 2007). Either way, the score equations for GLMs will ensure
that all terms inside $\mathbb{P}_{n}[\cdot]$ but $\hat{\mu}_{0}(X)$ have a
sample mean of zero, leaving only $\mathbb{P}_{n}[\hat{\mu}_{0}(X)]$, which
will reside in the parameter space of $\theta_{\overline{a}}$ if the latter
equals the range of the GLM specified for $\mu_{0}(x)$.

Alternatively, one can use the method of targeted maximum likelihood
estimation (TMLE; van Der Laan and Rubin 2006; Zheng and van der Laan 2012),
which, by fitting each of the outcome models in two steps, will also ensure a
zero sample mean for all terms inside $\mathbb{P}_{n}[\cdot]$ but
$\hat{\mu}_{0}(X)$. This approach does not require the first-step models to be
GLM and thus can be used with a wider range of outcome models. In our case, it
involves the following steps:

  1. 1.

For $k=K,\ldots 0$

    1. (a)

Using $\hat{\mu}_{k+1}^{\textup{tmle}}(X,\overline{M}_{k+1})$ (or, in the case
$k=K$, the observed outcome $Y$) as the response variable, obtain a first-step
regression-imputation estimate of $\mu_{k}(X,\overline{M}_{k})$;

    2. (b)

Fit a one-parameter GLM for the conditional mean of
$\hat{\mu}_{k+1}^{\textup{tmle}}(X,\overline{M}_{k+1})$ (or, in the case
$k=K$, the observed outcome $Y$), using $g(\hat{\mu}_{k}(X,\overline{M}_{k}))$
as an offset term and $\hat{w}_{k}(A,X,\overline{M}_{k})$ as the only
covariate (without an intercept term), and obtain an updated estimate
$\hat{\mu}_{k}^{\textup{tmle}}(X,\overline{M}_{k})=g^{-1}\big{(}g(\hat{\mu}_{k}(X,\overline{M}_{k}))+\text{$\hat{\beta}$}_{k}\hat{w}_{k}(A,X,\overline{M}_{k})\big{)}$,
where $g(\cdot)$ is the link function for the GLM and $\hat{\beta}_{k}$ is the
estimated coefficient on $\hat{w}_{k}(A,X,\overline{M}_{k})$;

  2. 2.

Obtain the final estimate
$\hat{\theta}_{\overline{a}}^{\textup{tmle}}=\mathbb{P}_{n}[\hat{\mu}_{0}^{\textup{tmle}}(X)]$.

In the one-mediator case, the above estimator is similar to the TMLE estimator
proposed by Zheng and van der Laan (2012) for the NDE, i.e.,
$\psi_{01}-\psi_{00}$. Since Zheng and van der Laan’s estimand is the NDE
instead of the mediation functional, their TMLE procedure involves fitting a
model for the “mediated mean outcome difference” (p. 6), i.e.,
$\mathbb{E}[\mathbb{E}[Y|X,A=1,M]-\mathbb{E}[Y|X,A=0,M]|X,A=0]$, instead of
the conditional mean of the imputed outcome itself, i.e., $\mu_{0}(X)$.

As with the GLM-based adjustments, the TMLE approach also yields a regression-
imputation estimator that resides in the parameter space of
$\theta_{\overline{a}}$ if the latter equals the range of the model specified
for $\mu_{0}(x)$. It should be noted that when data-adaptive methods are used
to obtain first-step estimates of the nuisance functions, sample splitting
should be employed so that steps 1(a) and steps 1(b) are implemented on
different subsamples. In cross-fitting, for example, steps 1(a) should be
implemented in the auxiliary sample ($S\backslash S_{j}$) and steps 1(b)
implemented in the main sample $S_{j}$. The method of TMLE can also be used to
adjust $\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$, in which case the
first step estimates of $\mu_{k}(X,\overline{M}_{k})$ ($0\leq k\leq K-1$) are
based on equation (13), and the weights $\hat{w}_{k}(A,X,\overline{M}_{k})$
($0\leq k\leq K$) reflect the corresponding terms in equation (12).

###  4 Special Cases

We have so far considered $\theta_{\overline{a}}$ for the unconstrained case
where $a_{1},\ldots a_{K+1}$ can each take 0 or 1. In many applications, the
researcher may be interested in particular causal estimands such as the
natural direct effect (NDE), the natural/total indirect effect (NIE/TIE), and
natural path-specific effects (nPSE; Daniel _et al._ 2015). Below, we discuss
how the multiply robust semiparametric estimators of $\theta_{\overline{a}}$
apply to these estimands. In addition, we discuss a set of cumulative path-
specific effects (cPSEs) that together compose the ATE. In Supplementary
Material E, we connect these cPSEs to noncausal decompositions of between-
group disparities that are widely used in the social sciences. For
illustrative purposes, we focus on estimators based on
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$, although similar results
hold for those based on $\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$.
Throughout this section, we maintain Assumptions 1*-3* so that
$\theta_{\overline{a}}=\psi_{\overline{a}}$.

BaselineNDE:$A$$M_{2}$$M_{1}$$Y$$\textup{NIE}_{M_{1}}$:$A$$M_{2}$$M_{1}$$Y$$\textup{TIE}_{M_{1}}$:$A$$M_{2}$$M_{1}$$Y$$\textup{nPSE}_{M_{2}}$:$A$$M_{2}$$M_{1}$$Y$$\textup{cPSE}_{M_{2}}$:$A$$M_{2}$$M_{1}$$Y$

Comparison$A$$M_{2}$$M_{1}$$Y$$A$$M_{2}$$M_{1}$$Y$$A$$M_{2}$$M_{1}$$Y$$A$$M_{2}$$M_{1}$$Y$$A$$M_{2}$$M_{1}$$Y$

Figure 2: Illustrations of NDE, NIE, TIE, nPSE, and cPSE in the case of two
mediators.

Note: $A$ denotes the treatment, $Y$ denotes the outcome of interest, and
$M_{1}$ and $M_{2}$ denote two causally ordered mediators. Solid and dashed
arrows for $A\to M_{1}$, $A\to M_{2}$, and $A\to Y$ denote activated ($A=1$)
and unactivated ($A=0$) paths, respectively. Gray arrows $M_{1}\to M_{2}$,
$M_{1}\to Y$, and $M_{2}\to Y$ signify that the mediators $M_{1}$ and $M_{2}$
are not under direct intervention.

#### 4.1 Natural Direct Effect (NDE)

The NDE measures the effect of switching treatment status from $0$ to $1$ in a
hypothetical world where the mediators $(M_{1},\ldots M_{K})$ were all set to
values they would have “naturally” taken for each unit under treatment status
$A=0$. It is thus given by
$\psi_{\overline{0}_{K},1}-\psi_{\overline{0}_{K+1}}$. The first row of Figure
2 illustrates the baseline and comparison interventions associated with the
NDE for the case of $K=2$, where the black solid and dashed arrows for $A\to
M_{1}$, $A\to M_{2}$, and $A\to Y$ denote activated ($A=1$) and unactivated
($A=0$) paths, respectively. A semiparametric efficient estimator for the NDE
can be constructed as

$\widehat{\textup{NDE}}^{\textup{eif}_{2}}=\hat{\theta}_{\overline{0}_{K},1}^{\textup{eif}_{2}}-\hat{\theta}_{\overline{0}_{K+1}}^{\textup{eif}_{2}}.$
(16)

If we treat $\overline{M}_{K}=(M_{1},\ldots M_{K})$ as a whole,
$\psi_{\overline{0}_{K},1}-\psi_{\overline{0}_{K+1}}$ coincides with the NDE
defined in the single mediator setting. In fact,
$\widehat{\textup{NDE}}^{\textup{eif}_{2}}$ is akin to the semiparametric
estimator of the NDE given in Zheng and van der Laan (2012). By contrast, if
we use $\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ instead of
$\hat{\theta}_{\overline{a}}^{\textup{eif}_{2}}$ in equation (16), we obtain
Tchetgen Tchetgen and Shpitser’s (2012) estimator of the NDE.

Setting $a_{1}=\ldots a_{K+1}=0$ in equation (14), we have

$\hat{\theta}_{\overline{0}_{K+1}}^{\textup{eif}_{2}}=\mathbb{P}_{n}\big{[}\frac{\mathbb{I}(A=0)}{\hat{\pi}_{0}(0|X)}\big{(}Y-\hat{\mu}_{0}(X)\big{)}+\hat{\mu}_{0}(X)\big{]},$
(17)

where $\mu_{0}(X)=\mathbb{E}[Y|X,A=0]$. Not surprisingly,
$\hat{\theta}_{\overline{0}_{K+1}}^{\textup{eif}_{2}}$ is the standard doubly
robust estimator for $\mathbb{E}[Y(0)]$, which is consistent if either
$\hat{\pi}_{0}(0|X)$ or $\hat{\mu}_{0}(X)$ is consistent. Similarly, by
setting $a_{1}=\ldots a_{K}=0$ and $a_{K+1}=1$ in equation (14), we have

$\hat{\theta}_{\overline{0}_{K},1}^{\textup{eif}_{2}}=\mathbb{P}_{n}\big{[}\frac{\mathbb{I}(A=1)}{\hat{\pi}_{0}(0|X)}\frac{\hat{\pi}_{K}(0|X,\overline{M}_{K})}{\hat{\pi}_{K}(1|X,\overline{M}_{K})}\big{(}Y-\hat{\mu}_{K}(X,\overline{M}_{K})\big{)}+\frac{\mathbb{I}(A=0)}{\hat{\pi}_{0}(0|X)}\big{(}\hat{\mu}_{K}(X,\overline{M}_{K})-\hat{\mu}_{0,K}(X)\big{)}+\hat{\mu}_{0,K}(X)\big{]}.$

In contrast to the general case where $\overline{a}_{K}$ is unconstrained,
$\hat{\theta}_{\overline{0}_{K},1}^{\textup{eif}_{2}}$ involves estimating
only four nuisance functions: $\pi_{0}(a|x)$, $\pi_{K}(a|x,\overline{m}_{K})$,
$\mu_{0,K}(x)$, and $\mu_{K}(x,\overline{m}_{K})$, where
$\mu_{K}(x,\overline{m}_{K})=\mathbb{E}[Y|x,A=1,\overline{m}_{K}]$ and
$\mu_{0,K}(x)=\mathbb{E}[\mu_{K}(X,\overline{M}_{K})|x,A=0]$. Hence
$\mu_{0,K}(x)$ can be estimated by fitting a model for the conditional mean of
$\hat{\mu}_{K}(X,\overline{M}_{K})$ given $(X,A)$ and then setting $A=0$ for
all units. It follows from Theorem 3 that
$\hat{\theta}_{\overline{0}_{K},1}^{\textup{eif}_{2}}$ is triply robust in
that it is consistent if one of the following three conditions holds: (a)
$\hat{\pi}_{0}$ and $\hat{\pi}_{K}$ are consistent; (b) $\hat{\pi}_{0}$ and
$\hat{\mu}_{K}$ are consistent; and (c) $\hat{\mu}_{0,K}$ and $\hat{\mu}_{K}$
are consistent. In the meantime, we know that
$\hat{\theta}_{\overline{0}_{K+1}}^{\textup{eif}_{2}}$ is consistent if either
$\hat{\pi}_{0}$ or $\hat{\mu}_{0}$ is consistent. By taking the intersection
of the multiple robustness conditions for
$\hat{\theta}_{\overline{0}_{K},1}^{\textup{eif}_{2}}$ and
$\hat{\theta}_{\overline{0}_{K+1}}^{\textup{eif}_{2}}$, we deduce that
$\widehat{\textup{NDE}}^{\textup{eif}_{2}}$ is also triply robust, as detailed
in Corollary 1.

######  Corollary 1.

Suppose all assumptions required for Theorem 4 hold. When the nuisance
functions are estimated via parametric models,
$\widehat{\textup{NDE}}^{\textup{eif}_{2}}$ is CAN provided that one of the
following three sets of nuisance functions is correctly specified and its
parameter estimates are $\sqrt{n}$-consistent: $\\{\pi_{0},\pi_{K}\\}$,
$\\{\pi_{0},\mu_{K}\\}$, $\\{\mu_{0},\mu_{0,K},\mu_{K}\\}$.
$\widehat{\textup{NDE}}^{\textup{eif}_{2}}$ is semiparametric efficient if all
of the above nuisance functions are correctly specified and their parameter
estimates $\sqrt{n}$-consistent. When the nuisance functions are estimated via
data-adaptive methods and cross-fitting,
$\widehat{\textup{NDE}}^{\textup{eif}_{2}}$ is CAN and semiparametric
efficient if all of the nuisance functions are consistently estimated and
$r_{n}(\hat{\pi}_{0})r_{n}(\hat{\mu}_{0,K})+r_{n}(\hat{\pi}_{K})r_{n}(\hat{\mu}_{K})+r_{n}(\hat{\pi}_{0})r_{n}(\hat{\mu}_{0})=o(n^{-1/2})$.

#### 4.2 Natural and Total Indirect Effects for $M_{1}$

In Section 2.1, we noted that $\psi_{100}-\psi_{000}$ and
$\psi_{111}-\psi_{011}$ correspond to the NIE and TIE for the first mediator
$M_{1}$ (illustrated in the second and third rows of Figure 2). This
correspondence extends naturally to the case of $K$ mediators, where the NIE
and TIE for $M_{1}$ are given by

$\textup{NIE}_{M_{1}}=\psi_{1,\underline{0}_{2}}-\psi_{\overline{0}_{K+1}},\quad\textup{TIE}_{M_{1}}=\psi_{\overline{1}_{K+1}}-\psi_{0,\underline{1}_{2}},$

where $\underline{0}_{2}=(0,\ldots 0)$ and $\underline{1}_{2}=(1,\ldots 1)$
are vectors of length $K$ representing the fact that $a_{2}=\ldots=a_{K+1}=0$
in $\textup{NIE}_{M_{1}}$ and $a_{2}=\ldots=a_{K+1}=1$ in
$\textup{TIE}_{M_{1}}$. Since $\textup{TIE}_{M_{1}}$ can be obtained by
switching the 0s and 1s in $\textup{NIE}_{M_{1}}$ and then flipping the sign,
we focus on $\textup{NIE}_{M_{1}}$ below, noting that analogous results hold
for $\textup{TIE}_{M_{1}}$.

A semiparametric efficient estimator of $\textup{NIE}_{M_{1}}$ can be
constructed as

$\widehat{\textup{NIE}}_{M_{1}}^{\textup{eif}_{2}}=\hat{\theta}_{1,\underline{0}_{2}}^{\textup{eif}_{2}}-\hat{\theta}_{\overline{0}_{K+1}}^{\textup{eif}_{2}}.$

As shown previously, $\hat{\theta}_{\overline{0}_{K+1}}^{\textup{eif}_{2}}$ is
given by the doubly robust estimator (17). Setting $a_{1}=1$ and $a_{2}=\ldots
a_{K+1}=0$ in equation (14), we obtain

$\displaystyle\hat{\theta}_{1,\underline{0}_{2}}^{\textup{eif}_{2}}=$
$\displaystyle\mathbb{P}_{n}\big{[}\frac{\mathbb{I}(A=0)}{\hat{\pi}_{0}(1|X)}\frac{\hat{\pi}_{1}(1|X,M_{1})}{\hat{\pi}_{1}(0|X,M_{1})}\big{(}Y-\hat{\mu}_{1}(X,M_{1})\big{)}+\frac{\mathbb{I}(A=1)}{\hat{\pi}_{0}(1|X)}\big{(}\hat{\mu}_{1}(X,M_{1})-\hat{\mu}_{0,1}(X)\big{)}+\hat{\mu}_{0,1}(X)\big{]}.$

Like $\hat{\theta}_{\overline{0}_{K},1}^{\textup{eif}_{2}}$,
$\hat{\theta}_{1,\underline{0}_{2}}^{\textup{eif}_{2}}$ also involves
estimating four nuisance functions: $\pi_{0}(a|x)$, $\pi_{1}(a|x,m_{1})$,
$\mu_{0,1}(x)$, and $\mu_{1}(x,m_{1})$, where
$\mu_{1}(x,m_{1})=\mathbb{E}[Y|x,A=0,m_{1}]$ and
$\mu_{0,1}(x)=\mathbb{E}[\mu_{1}(X,M_{1})|x,A=1]$. It follows from Theorem 3
that $\hat{\theta}_{1,\underline{0}_{2}}^{\textup{eif}_{2}}$ is triply robust
in that it is consistent if one of the following three conditions holds: (a)
$\hat{\pi}_{0}$ and $\hat{\pi}_{1}$ are consistent; (b) $\hat{\pi}_{0}$ and
$\hat{\mu}_{1}$ are consistent; and (c) $\hat{\mu}_{0,1}$ and $\hat{\mu}_{1}$
are consistent. By taking the intersection of the multiple robustness
conditions for $\hat{\theta}_{1,\underline{0}_{2}}^{\textup{eif}_{2}}$ and
$\hat{\theta}_{\overline{0}_{K+1}}^{\textup{eif}_{2}}$, we deduce that
$\widehat{\textup{NIE}}_{M_{1}}^{\textup{eif}_{2}}$ is also triply robust, as
detailed in Corollary 2.

######  Corollary 2.

Suppose all assumptions required for Theorem 4 hold. When the nuisance
functions are estimated via parametric models,
$\widehat{\textup{NIE}}_{M_{1}}^{\textup{eif}_{2}}$ is CAN provided that one
of the following three sets of nuisance functions is correctly specified and
its parameter estimates are $\sqrt{n}$-consistent: $\\{\pi_{0},\pi_{1}\\}$,
$\\{\pi_{0},\mu_{1}\\}$, $\\{\mu_{0},\mu_{0,1},\mu_{1}\\}$.
$\widehat{\textup{NIE}}_{M_{1}}^{\textup{eif}_{2}}$ is semiparametric
efficient if all of the above nuisance functions are correctly specified and
their parameter estimates $\sqrt{n}$-consistent. When the nuisance functions
are estimated via data-adaptive methods and cross-fitting,
$\widehat{\textup{NIE}}_{M_{1}}^{\textup{eif}_{2}}$ is semiparametric
efficient if all of the nuisance functions are consistently estimated and
$r_{n}(\hat{\pi}_{0})r_{n}(\hat{\mu}_{0,1})+r_{n}(\hat{\pi}_{1})r_{n}(\hat{\mu}_{1})+r_{n}(\hat{\pi}_{0})r_{n}(\hat{\mu}_{0})=o(n^{-1/2})$.

#### 4.3 Natural Path-Specific Effects (nPSE) for $M_{k}$ ($k\geq 2$)

In the same spirit of the NIE for $M_{1}$, the natural path-specific effect
(nPSE; Daniel _et al._ 2015) for mediator $M_{k}$ ($k\geq 2$) is defined as

$\textup{nPSE}_{M_{k}}=\psi_{\overline{0}_{k-1},1,\underline{0}_{k+1}}-\psi_{\overline{0}_{K+1}}.$

It can be interpreted as the effect of activating the path $A\to
M_{k}\rightsquigarrow Y$ while all other causal paths are “switched off,” as
shown in the fourth row of Figure 2. A semiparametric efficient estimator of
$\textup{nPSE}_{M_{k}}$ can be constructed as

$\widehat{\textup{nPSE}}_{M_{k}}^{\textup{eif}_{2}}=\hat{\theta}_{\overline{0}_{k-1},1,\underline{0}_{k+1}}^{\textup{eif}_{2}}-\hat{\theta}_{\overline{0}_{K+1}}^{\textup{eif}_{2}}.$

If, instead, we use $\hat{\theta}_{\overline{a}}^{\textup{eif}_{1}}$ in the
above equation, the resulting estimator
$\widehat{\textup{nPSE}}_{M_{k}}^{\textup{eif}_{1}}$ can be seen as Miles _et
al._ ’s (2020) estimator of $\theta_{010}-\theta_{000}$ applied to
$\widetilde{M}_{1}=(M_{1},M_{2},\ldots M_{k-1})$ and
$\widetilde{M}_{2}=M_{k}$.

Again, $\hat{\theta}_{\overline{0}_{K+1}}^{\textup{eif}_{2}}$ is given by the
doubly robust estimator (17). Setting $a_{1}=\ldots a_{k-1}=a_{k+1}=\ldots
a_{K+1}=0$ and $a_{k}=1$ in equation (14), we obtain

$\displaystyle\hat{\theta}_{\overline{0}_{k-1},1,\underline{0}_{k+1}}^{\textup{eif}_{2}}=$
$\displaystyle\mathbb{P}_{n}\big{[}\frac{\mathbb{I}(A=0)}{\hat{\pi}_{0}(0|X)}\frac{\hat{\pi}_{k-1}(0|X,\overline{M}_{k-1})}{\hat{\pi}_{k-1}(1|X,\overline{M}_{k-1})}\frac{\hat{\pi}_{k}(1|X,\overline{M}_{k})}{\hat{\pi}_{k}(0|X,\overline{M}_{k})}\big{(}Y-\hat{\mu}_{k}(X,\overline{M}_{k})\big{)}.$
$\displaystyle+\frac{\mathbb{I}(A=1)}{\hat{\pi}_{0}(0|X)}\frac{\hat{\pi}_{k-1}(0|X,\overline{M}_{k-1})}{\hat{\pi}_{k-1}(1|X,\overline{M}_{k-1})}\big{(}\hat{\mu}_{k}(X,\overline{M}_{k})-\hat{\mu}_{k-1,k}(X,\overline{M}_{k-1})\big{)}$
$\displaystyle+\frac{\mathbb{I}(A=0)}{\hat{\pi}_{0}(0|X)}\big{(}\hat{\mu}_{k-1,k}(X,\overline{M}_{k-1})-\hat{\mu}_{0,k-1,k}(X)\big{)}+\hat{\mu}_{0,k-1,k}(X)\big{]}.$

We can see that
$\hat{\theta}_{\overline{0}_{k-1},1,\underline{0}_{k+1}}^{\textup{eif}_{2}}$involves
estimating six nuisance functions: $\pi_{0}(a|x)$,
$\pi_{k-1}(a|x,\overline{m}_{k-1})$, $\pi_{k}(a|x,\overline{m}_{k})$,
$\mu_{0,k-1,k}(x)$, $\mu_{k-1,k}(x,\overline{m}_{k-1})$, and
$\mu_{k}(x,\overline{m}_{k})$, where
$\mu_{k}(X,\overline{M}_{k})=\mathbb{E}[Y|X,A=0,\overline{M}_{k}]$,
$\mu_{k-1,k}(X,\overline{M}_{k-1})=\mathbb{E}[\mu_{k}(X,\overline{M}_{k})|X,A=1,\overline{M}_{k-1}]$,
and $\mu_{0,k-1,k}(X)=\mathbb{E}[\mu_{k-1,k}(X,\overline{M}_{k-1})|X,A=0]$.
Hence $\mu_{k-1,k}(x)$ can be estimated by fitting a model for the conditional
mean of $\hat{\mu}_{k}(X,\overline{M}_{k})$ given ($X,A,\overline{M}_{k-1})$
and then setting $A=1$ for all units, and $\mu_{0,k-1,k}(x)$ can be estimated
by fitting a model for the conditional mean of
$\hat{\mu}_{k-1,k}(X,\overline{M}_{k-1})$ given ($X,A)$ and then setting $A=0$
for all units. It follows from Theorem 3 that
$\hat{\theta}_{\overline{0}_{k-1},1,\underline{0}_{k+1}}^{\textup{eif}_{2}}$
is quadruply robust in that it is consistent if one of the following four
conditions holds: (a) $\hat{\pi}_{0}$, $\hat{\pi}_{k-1}$, and $\hat{\pi}_{k}$
are consistent; (b) $\hat{\pi}_{0}$, $\hat{\pi}_{k-1}$, and $\hat{\mu}_{k}$
are consistent; (c) $\hat{\pi}_{0}$, $\hat{\mu}_{k-1,k}$, and $\hat{\mu}_{k}$
are consistent; and (d) $\hat{\mu}_{0,k-1,k}$, $\hat{\mu}_{k-1,k}$, and
$\hat{\mu}_{k}$ are consistent. By taking the intersection of the multiple
robustness conditions for
$\hat{\theta}_{1,\underline{0}_{2}}^{\textup{eif}_{2}}$ and
$\hat{\theta}_{\overline{0}_{K+1}}^{\textup{eif}_{2}}$, we deduce that
$\widehat{\textup{nPSE}}_{M_{k}}^{\textup{eif}_{2}}$ is also quadruply robust,
as detailed in Corollary 3.

######  Corollary 3.

Suppose all assumptions required for Theorem 4 hold. When the nuisance
functions are estimated via parametric models,
$\widehat{\textup{nPSE}}_{M_{k}}^{\textup{eif}_{2}}$ is CAN provided that one
of the following four sets of nuisance functions is correctly specified and
its parameter estimates are $\sqrt{n}$-consistent:
$\\{\pi_{0},\pi_{k-1},\pi_{k}\\}$, $\\{\pi_{0},\pi_{k-1},\mu_{k}\\}$,
$\\{\pi_{0},\mu_{k-1,k},\mu_{k}\\}$,
$\\{\mu_{0},\mu_{0,k-1,k},\mu_{k-1,k},\mu_{k}\\}$.
$\widehat{\textup{nPSE}}_{M_{k}}^{\textup{eif}_{2}}$ is semiparametric
efficient if all of the above nuisance functions are correctly specified and
their parameter estimates $\sqrt{n}$-consistent. When the nuisance functions
are estimated via data-adaptive methods and cross-fitting,
$\widehat{\textup{nPSE}}_{M_{k}}^{\textup{eif}_{2}}$ is semiparametric
efficient if all of the nuisance functions are consistently estimated and
$r_{n}(\hat{\pi}_{0})r_{n}(\hat{\mu}_{0,k-1,k})+r_{n}(\hat{\pi}_{k-1})r_{n}(\hat{\mu}_{k-1,k})+r_{n}(\hat{\pi}_{k})r_{n}(\hat{\mu}_{k})+r_{n}(\hat{\pi}_{0})r_{n}(\hat{\mu}_{0})=o(n^{-1/2})$.

#### 4.4 Cumulative Path-Specific Effects (cPSE) for $M_{k}$ ($k\geq 2$)

The NDE, NIE, and nPSE are all defined as the effect of activating one causal
path while keeping all other causal paths “switched off.” By contrast, in
equation (3), the ATE is decomposed into $K+1$ components, each of which
reflects the cumulative contribution of a specific mediator to the ATE.
Specifically, the component
$\psi_{\overline{0}_{K},1}-\psi_{\overline{0}_{K+1}}$equals the NDE, the
component $\psi_{\underline{1}_{1}}-\psi_{0,\underline{1}_{2}}$ equals
$\textup{TIE}_{M_{1}}$, and the component
$\psi_{\overline{0}_{k-1},\underline{1}_{k}}-\psi_{\overline{0}_{k},\underline{1}_{k+1}}$
gauges the additional contribution of the causal path $A\to
M_{k}\rightsquigarrow Y$ after the causal paths $A\to M_{k+1}\rightsquigarrow
Y,\ldots A\to M_{K}\rightsquigarrow Y,A\to Y$ are “switched on.” Such a
decomposition will be useful in applications where the investigator aims to
partition the ATE into its path-specific components.

We define the cumulative path-specific effect (cPSE) for mediator $M_{k}$
($k\geq 2$) as

$\textup{cPSE}_{M_{k}}=\psi_{\overline{0}_{k-1},\underline{1}_{k}}-\psi_{\overline{0}_{k},\underline{1}_{k+1}}.$

The last row of Figure 2 gives the baseline and comparison interventions
associated with $\textup{cPSE}_{M_{2}}$ in the case of $K=2$. A semiparametric
efficient estimator for $\textup{cPSE}_{M_{k}}$ can be constructed as

$\widehat{\textup{cPSE}}_{M_{k}}^{\textup{eif}_{2}}=\hat{\theta}_{\overline{0}_{k-1},\underline{1}_{k}}^{\textup{eif}_{2}}-\hat{\theta}_{\overline{0}_{k},\underline{1}_{k+1}}^{\textup{eif}_{2}}.$

Setting $a_{1}=\ldots a_{k}=0$ and $a_{k+1}=\ldots=a_{K+1}=1$ in equation
(14), we obtain

$\displaystyle\hat{\theta}_{\overline{0}_{k},\underline{1}_{k+1}}^{\textup{eif}_{2}}$
$\displaystyle=\mathbb{P}_{n}\big{[}\frac{\mathbb{I}(A=1)}{\hat{\pi}_{0}(0|X)}\frac{\hat{\pi}_{k}(0|X,\overline{M}_{k})}{\hat{\pi}_{k}(1|X,\overline{M}_{k})}\big{(}Y-\hat{\mu}_{k}(X,\overline{M}_{k})\big{)}+\frac{\mathbb{I}(A=0)}{\hat{\pi}_{0}(0|X)}\big{(}\hat{\mu}_{k}(X,\overline{M}_{k})-\hat{\mu}_{0,k}(X)\big{)}+\hat{\mu}_{0,k}(X)\big{]},$
(18)

where $\mu_{k}(X,\overline{M}_{k})=\mathbb{E}[Y|X,A=1,\overline{M}_{k}]$ and
$\mu_{0,k}(X)=\mathbb{E}[\mu_{k}(X,\overline{M}_{k})|X,A=0]$. It follows from
Theorem 3 that
$\hat{\theta}_{\overline{0}_{k},\underline{1}_{k+1}}^{\textup{eif}_{2}}$ is
triply robust in that it is consistent if one of the following three
conditions holds: (a) $\hat{\pi}_{0}$ and $\hat{\pi}_{k}$ are consistent; (b)
$\hat{\pi}_{0}$ and $\hat{\mu}_{k}$ are consistent; and (c) $\hat{\mu}_{0,k}$
and $\hat{\mu}_{k}$ are consistent. By replacing $k$ with $k-1$ in equation
(18), we obtain a similar expression for
$\hat{\theta}_{\overline{0}_{k-1},\underline{1}_{k}}^{\textup{eif}_{2}}$,
which is also triply robust in that it is consistent if one of the following
three conditions holds: (a) $\hat{\pi}_{0}$ and $\hat{\pi}_{k-1}$ are
consistent; (b) $\hat{\pi}_{0}$ and $\hat{\mu}_{k-1}$ are consistent; and (c)
$\hat{\mu}_{0,k-1}$ and $\hat{\mu}_{k-1}$ are consistent. As a result,
$\widehat{\textup{cPSE}}_{M_{k}}^{\textup{eif}_{2}}$ involves fitting seven
working models — for $\pi_{0}(a|x)$, $\pi_{k-1}(a|x,\overline{m}_{k-1})$,
$\pi_{k}(a|x,\overline{m}_{k})$, $\mu_{k-1}(x,\overline{m}_{k-1})$,
$\mu_{0,k-1}(x)$, $\mu_{k}(x,\overline{m}_{k})$, and $\mu_{0,k}(x)$. By taking
the intersection of the multiple robustness conditions for
$\hat{\theta}_{\overline{1}_{k},\underline{0}_{k+1}}^{\textup{eif}_{2}}$and
$\hat{\theta}_{\overline{1}_{k-1},\underline{0}_{k}}^{\textup{eif}_{2}}$, we
deduce that $\widehat{\textup{cPSE}}_{M_{k}}^{\textup{eif}_{2}}$ is quintuply
robust in that it is consistent if one of five sets of nuisance functions is
correctly specified and consistently estimated, as detailed in Corollary 4.

######  Corollary 4.

Suppose all assumptions required for Theorem 4 hold. When the nuisance
functions are estimated via parametric models,
$\widehat{\textup{cPSE}}_{M_{k}}^{\textup{eif}_{2}}$ is CAN provided that one
of the following five sets of nuisance functions is correctly specified and
its parameter estimates are $\sqrt{n}$-consistent:
$\\{\pi_{0},\pi_{k-1},\pi_{k}\\}$;$\\{\pi_{0},\pi_{k-1},\mu_{k}\\}$;$\\{\pi_{0},\mu_{k-1},\pi_{k}\\}$;$\\{\pi_{0},\mu_{k-1},\mu_{k}\\}$;$\\{\mu_{0,k-1},\mu_{0,k},\mu_{k-1},\mu_{k}\\}$.
$\widehat{\textup{cPSE}}_{M_{k}}^{\textup{eif}_{2}}$ is semiparametric
efficient if all of the above nuisance functions are correctly specified and
their parameter estimates $\sqrt{n}$-consistent. When the nuisance functions
are estimated via data-adaptive methods and cross-fitting,
$\widehat{\textup{cPSE}}_{M_{k}}^{\textup{eif}_{2}}$ is semiparametric
efficient if all of the nuisance functions are consistently estimated and
$r_{n}(\hat{\pi}_{0})r_{n}(\hat{\mu}_{0,k-1})+r_{n}(\hat{\pi}_{0})r_{n}(\hat{\mu}_{0,k})+r_{n}(\hat{\pi}_{k-1})r_{n}(\hat{\mu}_{k-1})+r_{n}(\hat{\pi}_{k})r_{n}(\hat{\mu}_{k})=o(n^{-1/2})$.

###  5 A Simulation Study

In this section, we conduct a simulation study to demonstrate the robustness
of various estimators under different forms of model misspecification.
Specifically, we consider a binary treatment $A$, a continuous outcome $Y$,
two causally ordered mediators $M_{1}$ and $M_{2}$, and four pretreatment
covariates $X_{1},X_{2},X_{3},X_{4}$ generated from the following model:

$\displaystyle(U_{1},U_{2},U_{3},U_{XY})$ $\displaystyle\sim N(0,I_{4}),$
$\displaystyle X_{j}$ $\displaystyle\sim
N((U_{1},U_{2},U_{3},U_{XY})\beta_{X_{j}},1),\quad j=1,2,3,4,$ $\displaystyle
A$
$\displaystyle\sim\textup{Bernoulli}\big{(}\textup{logit}^{-1}[(1,X_{1},X_{2},X_{3},X_{4})\beta_{A}]\big{)},$
$\displaystyle M_{1}$ $\displaystyle\sim
N\big{(}(1,X_{1},X_{2},X_{3},X_{4},A)\beta_{M_{1}},1\big{)},$ $\displaystyle
M_{2}$ $\displaystyle\sim
N\big{(}(1,X_{1},X_{2},X_{3},X_{4},A,M_{1})\beta_{M_{2}},1\big{)},$
$\displaystyle Y$ $\displaystyle\sim
N\big{(}(1,U_{XY},X_{1},X_{2},X_{3},X_{4},A,M_{1},M_{2})\beta_{Y},1\big{)}.$

The coefficients $\beta_{X_{j}}(\text{1$\leq j\leq 4$})$, $\beta_{A}$,
$\beta_{M_{1}}$, $\beta_{M_{2}},$ $\beta_{Y}$ are produced from a set of
uniform distributions (see Supplementary Material F for more details). Given
the coefficients, we generate 1,000 Monte Carlo samples of size 2,000. Note
that in the above model, the unobserved variable $U_{XY}$ confounds the
$X$-$Y$ relationship but does not pose an identification threat for
$\psi_{\overline{a}}$ and the associated PSEs (i.e., Assumption 2 still
holds).

Without loss of generality, we focus on the estimand $\textup{cPSE}_{M_{2}}$,
which we estimate by $\hat{\theta}_{011}-\hat{\theta}_{001}$. To highlight the
general results stated in Theorem 3, we use only estimators for the generic
$\theta_{\overline{a}}$ (i.e., those described in Section 3). First, we
consider the weighting estimator $\hat{\theta}_{\overline{a}}^{\textup{w-a}}$,
the regression-imputation estimator
$\hat{\theta}_{\overline{a}}^{\textup{ri}}$, and the hybrid estimators
$\hat{\theta}_{\overline{a}}^{\textup{ri-w-w}}$ and
$\hat{\theta}_{\overline{a}}^{\textup{ri-ri-w}}$, where the mediator density
ratio involved in $\hat{\theta}_{\overline{a}}^{\textup{ri-w-w}}$ is estimated
via the corresponding odds ratio of the treatment variable. We then consider
four EIF-based estimators
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{par,eif}_{2}}$,
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{par}_{2},\textup{eif}_{2}}$,
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{np,eif}_{2}}$, and
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{tmle,eif}_{2}}$. For
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{par,eif}_{2}}$ and
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{par}_{2},\textup{eif}_{2}}$,
the nuisance functions are estimated via GLMs.
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{par}_{2},\textup{eif}_{2}}$
differs from $\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{par,eif}_{2}}$
in that the outcome models $\mu_{2}(x,m_{1},m_{2})$, $\mu_{1}(x,m_{1})$, and
$\mu_{0}(x)$ are fitted using a set of weighted GLMs such that in equation
(15), all terms inside $\mathbb{P}_{n}[\cdot]$ but $\hat{\mu}_{0}(X)$ have a
zero sample mean, yielding a regression-imputation estimator that may perform
better in finite samples.

All of the above estimators are constructed using estimates of six nuisance
functions: $\pi_{0}(a|x)$, $\pi_{1}(a|x,m_{1})$, $\pi_{2}(a|x,m_{1},m_{2})$,
$\mu_{0}(x)$, $\mu_{1}(x,m_{1})$, and $\mu_{2}(x,m_{1},m_{2})$. To demonstrate
the consequences of model misspecification and the multiple robustness of
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{par,eif}_{2}}$ and
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{par}_{2},\textup{eif}_{2}}$,
we generate a set of “false covariates”
$Z=\big{(}X_{1},e^{X_{2}/2},(X_{3}/X_{1})^{1/3},X_{4}/(e^{X_{1}/2}+1)\big{)}$
and use them to fit a misspecified GLM for each of the nuisance functions
(with only the main effects of $Z_{1},Z_{2},Z_{3},Z_{4}$). We evaluate each of
the parametric estimators under five different cases: (a) only $\pi_{0}$,
$\pi_{1}$, $\pi_{2}$ are correctly specified; (b) only $\pi_{0}$, $\pi_{1}$,
$\mu_{2}$ are correctly specified; (c) only $\pi_{0}$, $\mu_{1}$, $\mu_{2}$
are correctly specified; (d) only $\mu_{0}$, $\mu_{1}$, $\mu_{2}$ are
correctly specified; and (e) all of the six nuisance functions are
misspecified. In theory, $\hat{\theta}_{\overline{a}}^{\textup{w-a}}$ is
consistent in case (a), $\hat{\theta}_{\overline{a}}^{\textup{ri-w-w}}$ is
consistent in case (b), $\hat{\theta}_{\overline{a}}^{\textup{ri-ri-w}}$ is
consistent in case (c), $\hat{\theta}_{\overline{a}}^{\textup{ri}}$ is
consistent in case (d), and
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{par,eif}_{2}}$ and
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{par}_{2},\textup{eif}_{2}}$
are consistent in cases (a)-(d). The corresponding estimators of
$\textup{cPSE}_{M_{2}}$ should follow the same properties.

For the two nonparametric estimators,
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{np,eif}_{2}}$ is based on
estimating equation (14), and
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{tmle,eif}_{2}}$ is based on
the method of TMLE. Like
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{par}_{2},\textup{eif}_{2}}$,
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{tmle,eif}_{2}}$ is a
regression-imputation estimator, which may have better finite-sample
performance than
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{np,eif}_{2}}$. For both
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{np,eif}_{2}}$ and
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{tmle,eif}_{2}}$, the nuisance
functions are estimated via a super learner (van der Laan _et al._ 2007)
composed of Lasso and random forest, where the feature matrix consists of
first-order, second-order, and interaction terms of the false covariates $Z$.
The super learner is more flexible than a misspecified GLM consisting of only
the main effects of $Z$, but it remains agnostic about the true nuisance
functions, which are either logit or linear models that depend on
$X=(Z_{1},2\log(Z_{2}),Z_{1}Z_{3}^{3},(1+e^{Z_{1}/2})Z_{4})$. We obtain
nonparametric estimates of $\textup{cPSE}_{M_{2}}$ using both five-fold cross-
fitting and no cross-fitting.

Figure 3: Sampling distributions of eight different estimators for $n=2,000$.
Cases (a)-(e) are described in the main text. The symbols y and n denote
whether cross-fitting is used to implement the nonparametric estimators (y =
yes, n = no).

Results from the simulation study are shown in Figure 3, where each panel
corresponds to an estimator, and the y axis is recentered at the true value of
$\textup{cPSE}_{M_{2}}$. The shaded box plots highlight cases under which a
given estimator should be consistent, and the box plots with a lighter shade
in the last two panels denote nonparametric estimators obtained without cross-
fitting. From the first four panels, we can see that the weighting,
regression-imputation, and hybrid estimators all behave as expected. They
center around the true value if the requisite nuisance functions are all
correctly specified, and deviate from the truth in most other cases. The next
four panels show the box plots of the EIF-based estimators. As expected, both
of the parametric EIF-based estimators are quadruply robust, as their sampling
distributions roughly concentrate around the true value in all of the four
cases from (a) to (d). Moreover, it is reassuring to see that when all of the
nuisance functions are misspecified (case (e)), the multiply robust estimators
do not show a larger amount of bias than those of the other parametric
estimators. Finally, both of the nonparametric EIF-based estimators perform
reasonably well. When cross-fitting is used, the estimating equation estimator
$\widehat{\textup{cPSE}}_{M_{2}}^{\textup{np},\textup{eif}_{2}}$ appears to
have a smaller bias than the TMLE estimator
$\widehat{\textup{cPSE}}_{M_{2}}^{\textup{tmle},\textup{eif}_{2}}$, but it
occasionally gives rise to extreme estimates. Their 95% Wald confidence
intervals, constructed using the estimated variance
$\hat{\mathbb{E}}\big{[}\big{(}\hat{\varphi}_{011}-\hat{\varphi}_{001}\big{)}^{2}\big{]}/n$,
have close-to-nominal coverage rates — 95.5% for
$\widehat{\textup{cPSE}}_{M_{2}}^{\textup{np},\textup{eif}_{2}}$ and 90.9% for
$\widehat{\textup{cPSE}}_{M_{2}}^{\textup{tmle},\textup{eif}_{2}}$. Without
cross-fitting, the point estimates exhibit similar distributions, but the
coverage rates of the corresponding 95% confidence intervals are somewhat
lower — 87.3% for
$\widehat{\textup{cPSE}}_{\textup{$M_{2}$}}^{\textup{np},\textup{eif}_{2}}$
and 85.8% for
$\widehat{\textup{cPSE}}_{M_{2}}^{\textup{tmle},\textup{eif}_{2}}$.

###  6 An Empirical Application

In this section, we illustrate semiparametric estimation of PSEs by analyzing
the causal pathways through which higher education affects political
participation. Prior research suggests that college attendance has a
substantial positive effect on political participation in the United States
(e.g., Dee 2004; Milligan _et al._ 2004). Yet, the mechanisms underlying this
causal link remain unclear. The effect of college on political participation
may operate through the development of civic and political interest (e.g.,
Hillygus 2005), through an increase in economic status (e.g., Kingston _et
al._ 2003), or through other pathways such as social and occupational networks
(e.g., Rolfe 2012). To examine these direct and indirect effects, we consider
a causal structure akin to the top panel of Figure 1, where $A$ denotes
college attendance, $Y$ denotes political participation, and $M_{1}$ and
$M_{2}$ denote two causally ordered mediators that reflect (a) economic
status, and (b) civic and political interest, respectively.

In this model, economic status is allowed to affect civic and political
interest but not vice versa, which we consider to be a reasonable
approximation to reality. Nonetheless, the conditional independence assumption
(Assumption 2) is still strong in this context, as it rules out unobserved
confounding for any of the pairwise relationships between college attendance,
economic status, civic and political interest, and political participation.
Thus, the following analyses should be viewed as an illustration of the
proposed methodology rather than a definitive assessment of the PSEs of
interest.

We use data from $n=2,969$ individuals in the National Longitudinal Survey of
Youth 1997 (NLSY97) who were age 15-17 in 1997 and had completed high school
by age 20. The treatment $A$ is a binary indicator for whether the individual
had attended a two-year or four-year college by age 20. The outcome $Y$ is a
binary indicator for whether the individual voted in the 2010 general
election. We measure economic status ($M_{1}$) using the respondent’s average
annual earnings from 2006 to 2009. To gauge civic and political interest
($M_{2}$), we use a set of variables that reflect the respondent’s interest in
government and public affairs and involvement in volunteering, donation,
community group activities between 2007 and 2010. The overlap of the periods
in which $M_{1}$ and $M_{2}$ were measured is a limitation of this analysis,
and it makes our earlier assumption that $M_{2}$ does not affect $M_{1}$
essential for identifying the direct and path-specific effects.

To minimize potential bias due to unobserved confounding, we include a rich
set of pre-college individual and contextual characteristics in the vector of
pretreatment covariates $X$. They include gender, race, ethnicity, age at
1997, parental education, parental income, parental assets, presence of a
father figure, co-residence with both biological parents, percentile score on
the Armed Services Vocational Aptitude Battery (ASVAB), high school GPA, an
index of substance use (ranging from 0 to 3), an index of delinquency (ranging
from 0 to 10), whether the respondent had any children by age 18, college
expectation among the respondent’s peers, and a number of school-level
characteristics. Descriptive statistics on these pre-college characteristics
as well as the mediators and the outcome are given in Supplementary Material
G. Some components of $X$, $M_{1}$, and $M_{2}$ contain a small fraction of
missing values. They are imputed via a random-forest-based multiple imputation
procedure (with ten imputed data sets). The standard errors of our parameter
estimates are adjusted using Rubin’s (1987) method.

Under Assumptions 1-3 given in Section 2.1, a set of PSEs reflecting the
causal paths $A\to Y$, $A\to M_{1}\rightsquigarrow Y$, and $A\to M_{2}\to Y$
are identified. For illustrative purposes, we focus on the cumulative PSEs
(cPSEs) defined in Section 4.4:

$\textup{ATE}=\underbrace{\psi_{001}-\psi_{000}}_{A\to
Y}+\underbrace{\psi_{011}-\psi_{001}}_{A\to M_{2}\to
Y}+\underbrace{\psi_{111}-\psi_{011}}_{A\to M_{1}\rightsquigarrow Y}.$ (19)

Here, the first component is the NDE of college attendance, and the second and
third components reflect the amounts of treatment effect that are additionally
mediated by civic/political interest and economic status, respectively. Since
$M_{2}$ is multivariate, it would be difficult to model its conditional
distributions directly. We thus estimate the PSEs using the estimator
$\hat{\theta}_{a_{1},a_{2},a}^{\textup{eif}_{2}}$. Each of the nuisance
functions is estimated using a super learner composed with Lasso and random
forest. For computational reasons, the feature matrix supplied to the super
learner consists of only first-order terms of the corresponding variables. As
in our simulation study, we implement two versions of this EIF-based
estimator, one based on the original estimating equation
($\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{np,eif}_{2}}$), and one
based on the method of TMLE
($\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{tmle,eif}_{2}}$). Five-fold
cross-fitting is used to obtain the final estimates.

Table 1: Estimates of total and path-specific effects of college attendance on
voting.

| Estimating equation ($\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{np,eif}_{2}}$) | TMLE ($\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{tmle,eif}_{2}}$)
---|---|---
Average total effect | 0.152 (0.022) | 0.156 (0.023)
Through economic status ($A\to M_{1}\rightsquigarrow Y$) | 0.007 (0.005) | 0.002 (0.005)
Through civic/political interest ($A\to M_{2}\to Y$) | 0.042 (0.008) | 0.049 (0.008)
Direct effect ($A\to Y$) | 0.103 (0.021) | 0.105 (0.021)

Note: Numbers in parentheses are estimated standard errors, which are
constructed using sample variances of the estimated efficient influence
functions and adjusted for multiple imputation via Rubin’s (1987) method.

The results are shown in Table 1. We can see that the two estimators yield
similar estimates of the total and path-specific effects. By
$\hat{\theta}_{\textup{$\overline{a}$}}^{\textup{np,eif}_{2}}$, for example,
the estimated total effect of college attendance on voting is 0.152, meaning
that, on average, college attendance increases the likelihood of voting in
2010 by about 15 percentage points. The estimated PSE via $M_{2}$ is 0.042,
suggesting that a small fraction of the college effect operates through the
development of civic and political interest. By contrast, the estimated PSE
via economic status is substantively negligible and statistically
insignificant. A large portion of the college effect appears to be “direct,”
i.e., operating neither through increased economic status nor through
increased civic and political interest.

###  7 Concluding Remarks

By considering the general case of $K$($\geq 1$) causally ordered mediators,
this paper offers several new insights into the identification and estimation
of PSEs. First, under the assumptions associated with Pearl’s NPSEM with
mutually independent errors, we have defined a set of PSEs as contrasts
between the expectations of $2^{K+1}$ potential outcomes, which are identified
via what we call the generalized mediation functional (GMF). Second, building
on its efficient influence function, we have developed two $K+2$-robust and
semiparametric efficient estimators for the GMF. By virtue of their multiple
robustness, these estimators are well suited to the use of data-adaptive
methods for estimating their nuisance functions. For such cases, we have
established the rate conditions required of the nuisance functions for
consistency and semiparametric efficiency.

As we have seen, our proposed methodology is general in that the GMF
encompasses a variety of causal estimands such as the NDE, NIE/TIE, nPSE,
cPSE. Nonetheless, it does not accommodate PSEs that are not identified under
Pearl’s NPSEM, some of which may be scientifically important. For example,
social and biomedical scientists are often interested in testing hypotheses
about “serial mediation,” i.e., the degree to which the effect of a treatment
operates through multiple mediators sequentially, such as that reflected in
the causal path $A\to M_{1}\to M_{2}\to Y$ (e.g., Jones _et al._ 2015). Given
that the corresponding PSEs are not nonparametrically identified under Pearl’s
NPSEM, previous research has proposed strategies that involve either
additional assumptions (Albert and Nelson 2011) or alternative estimands (Lin
and VanderWeele 2017). We consider semiparametric estimation and inference for
these alternative approaches a promising direction for future research.

### References

  * Albert (2012) Albert, J. M. (2012) Mediation analysis for nonlinear models with confounding.  _Epidemiology_ , 23, 879.
  * Albert and Nelson (2011) Albert, J. M. and Nelson, S. (2011) Generalized causal mediation analysis.  _Biometrics_ , 67, 1028–1038.
  * Alwin and Hauser (1975) Alwin, D. F. and Hauser, R. M. (1975) The decomposition of effects in path analysis.  _American Sociological Review_ , 37–47.
  * Avin _et al._ (2005) Avin, C., Shpitser, I. and Pearl, J. (2005) Identifiability of path-specific effects.  In _Proceedings of the 19th International Joint Conference on Artificial Intelligence_ , 357–363. Morgan Kaufmann Publishers Inc.
  * Bang and Robins (2005) Bang, H. and Robins, J. M. (2005) Doubly robust estimation in missing data and causal inference models.  _Biometrics_ , 61, 962–973.
  * Baron and Kenny (1986) Baron, R. M. and Kenny, D. A. (1986) The moderator–mediator variable distinction in social psychological research: Conceptual, strategic, and statistical considerations.  _Journal of Personality and Social Psychology_ , 51, 1173\.
  * Chernozhukov _et al._ (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W. and Robins, J. (2018) Double/debiased machine learning for treatment and structural parameters.  _The Econometrics Journal_ , 21, C1–C68.
  * Daniel _et al._ (2015) Daniel, R., De Stavola, B., Cousens, S. and Vansteelandt, S. (2015) Causal mediation analysis with multiple mediators.  _Biometrics_ , 71, 1–14.
  * Dee (2004) Dee, T. S. (2004) Are there civic returns to education?  _Journal of Public Economics_ , 88, 1697–1720.
  * Fortin _et al._ (2011) Fortin, N., Lemieux, T. and Firpo, S. (2011) Decomposition methods in economics.  In _Handbook of Labor Economics_ , vol. 4, 1–102. Elsevier.
  * Goetgeluk _et al._ (2009) Goetgeluk, S., Vansteelandt, S. and Goetghebeur, E. (2009) Estimation of controlled direct effects.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 70, 1049–1066.
  * Hafeman and VanderWeele (2011) Hafeman, D. M. and VanderWeele, T. J. (2011) Alternative assumptions for the identification of direct and indirect effects.  _Epidemiology_ , 753–764.
  * Han (2014) Han, P. (2014) Multiply robust estimation in regression analysis with missing data.  _Journal of the American Statistical Association_ , 109, 1159–1173.
  * Han and Wang (2013) Han, P. and Wang, L. (2013) Estimation with missing data: Beyond double robustness.  _Biometrika_ , 100, 417–430.
  * Hillygus (2005) Hillygus, D. S. (2005) The missing link: Exploring the relationship between higher education and political engagement.  _Political Behavior_ , 27, 25–47.
  * Imai _et al._ (2010) Imai, K., Keele, L., Yamamoto, T. _et al._ (2010) Identification, inference and sensitivity analysis for causal mediation effects.  _Statistical Science_ , 25, 51–71.
  * Jones _et al._ (2015) Jones, C. L., Jensen, J. D., Scherr, C. L., Brown, N. R., Christy, K. and Weaver, J. (2015) The health belief model as an explanatory framework in communication research: Exploring parallel, serial, and moderated mediation.  _Health Communication_ , 30, 566–576.
  * Kang and Schafer (2007) Kang, J. D. and Schafer, J. L. (2007) Demystifying double robustness: A comparison of alternative strategies for estimating a population mean from incomplete data.  _Statistical Science_ , 22, 523–539.
  * Kingston _et al._ (2003) Kingston, P. W., Hubbard, R., Lapp, B., Schroeder, P. and Wilson, J. (2003) Why education matters.  _Sociology of Education_ , 53–70.
  * Lin and VanderWeele (2017) Lin, S.-H. and VanderWeele, T. (2017) Interventional approach for path-specific effects.  _Journal of Causal Inference_ , 5.
  * Luedtke _et al._ (2017) Luedtke, A. R., Sofrygin, O., van der Laan, M. J. and Carone, M. (2017) Sequential double robustness in right-censored longitudinal models.  _arXiv preprint arXiv:1705.02459_.
  * Miles _et al._ (2017) Miles, C. H., Shpitser, I., Kanki, P., Meloni, S. and Tchetgen Tchetgen, E. J. (2017) Quantifying an adherence path-specific effect of antiretroviral therapy in the nigeria pepfar program.  _Journal of the American Statistical Association_ , 112, 1443–1452.
  * Miles _et al._ (2020) Miles, C. H., Shpitser, I., Kanki, P., Meloni, S. and Tchetgen Tchetgen, E. J. (2020) On semiparametric estimation of a path-specific effect in the presence of mediator-outcome confounding.  _Biometrika_ , 107, 159–172.
  * Milligan _et al._ (2004) Milligan, K., Moretti, E. and Oreopoulos, P. (2004) Does education improve citizenship? evidence from the united states and the united kingdom.  _Journal of Public Economics_ , 88, 1667–1695.
  * Molina _et al._ (2017) Molina, J., Rotnitzky, A., Sued, M. and Robins, J. (2017) Multiple robustness in factorized likelihood models.  _Biometrika_ , 104, 561–581.
  * Newey and McFadden (1994) Newey, K. and McFadden, D. (1994) Large sample estimation and hypothesis.  _Handbook of Econometrics, IV, Edited by RF Engle and DL McFadden_ , 2112–2245.
  * Newey and Robins (2018) Newey, W. K. and Robins, J. R. (2018) Cross-fitting and fast remainder rates for semiparametric estimation.  _arXiv preprint arXiv:1801.09138_.
  * Neyman (1923) Neyman, J. S. (1923) On the application of probability theory to agricultural experiments. essay on principles. section 9.  _Annals of Agricultural Sciences_ , 10, 1–51.
  * Pearl (2001) Pearl, J. (2001) Direct and indirect effects.  In _Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence_ , 411–420. Morgan Kaufmann Publishers Inc.
  * Pearl (2009) Pearl, J. (2009) _Causality: Models, Reasoning, and Inference_.  Cambridge University Press.
  * Petersen _et al._ (2012) Petersen, M. L., Porter, K. E., Gruber, S., Wang, Y. and Van Der Laan, M. J. (2012) Diagnosing and responding to violations in the positivity assumption.  _Statistical Methods in Medical Research_ , 21, 31–54.
  * Petersen _et al._ (2006) Petersen, M. L., Sinisi, S. E. and van der Laan, M. J. (2006) Estimation of direct causal effects.  _Epidemiology_ , 17, 276–284.
  * Robins _et al._ (2007) Robins, J., Sued, M., Lei-Gomez, Q. and Rotnitzky, A. (2007) Comment: Performance of double-robust estimators when “inverse probability” weights are highly variable.  _Statistical Science_ , 22, 544–559.
  * Robins (2003) Robins, J. M. (2003) Semantics of causal dag models and the identification of direct and indirect effects.  _Highly Structured Stochastic Systems_ , 70–81.
  * Robins and Greenland (1992) Robins, J. M. and Greenland, S. (1992) Identifiability and exchangeability for direct and indirect effects.  _Epidemiology_ , 3, 143–155.
  * Robins and Richardson (2010) Robins, J. M. and Richardson, T. S. (2010) Alternative graphical causal models and the identification of direct effects.  _Causality and psychopathology: Finding the determinants of disorders and their cures_ , 103–158.
  * Rolfe (2012) Rolfe, M. (2012) _Voter Turnout: A Social Theory of Political Participation_.  Cambridge University Press.
  * Rotnitzky _et al._ (2017) Rotnitzky, A., Robins, J. and Babino, L. (2017) On the multiply robust estimation of the mean of the g-functional.  _arXiv preprint arXiv:1705.08582_.
  * Rubin (1974) Rubin, D. B. (1974) Estimating Causal Effects of Treatments in Randomized and Nonrandomized Studies.  _Journal of Educational Psychology_ , 66, 688–701.
  * Rubin (1987) Rubin, D. B. (1987) _Multiple Imputation for Nonresponse in Surveys_.  Wiley and Sons.
  * Scharfstein _et al._ (1999) Scharfstein, D. O., Rotnitzky, A. and Robins, J. M. (1999) Adjusting for nonignorable drop-out using semiparametric nonresponse models.  _Journal of the American Statistical Association_ , 94, 1096–1120.
  * Seaman and Vansteelandt (2018) Seaman, S. R. and Vansteelandt, S. (2018) Introduction to double robust methods for incomplete data.  _Statistical Science_ , 33, 184.
  * Shpitser (2013) Shpitser, I. (2013) Counterfactual graphical models for longitudinal mediation analysis with unobserved confounding.  _Cognitive Science_ , 37, 1011–1035.
  * Steen _et al._ (2017) Steen, J., Loeys, T., Moerkerke, B. and Vansteelandt, S. (2017) Flexible mediation analysis with multiple mediators.  _American Journal of Epidemiology_ , 186, 184–193.
  * Tchetgen Tchetgen (2013) Tchetgen Tchetgen, E. J. (2013) Inverse odds ratio-weighted estimation for causal mediation analysis.  _Statistics in Medicine_ , 32, 4567–4580.
  * Tchetgen Tchetgen and Shpitser (2012) Tchetgen Tchetgen, E. J. and Shpitser, I. (2012) Semiparametric theory for causal mediation analysis: Efficiency bounds, multiple robustness, and sensitivity analysis.  _Annals of Statistics_ , 40, 1816.
  * van der Laan _et al._ (2007) van der Laan, M. J., Polley, E. C. and Hubbard, A. E. (2007) Super learner.  _Statistical Applications in Genetics and Molecular Biology_ , 6.
  * van Der Laan and Rubin (2006) van Der Laan, M. J. and Rubin, D. (2006) Targeted maximum likelihood learning.  _The International Journal of Biostatistics_ , 2.
  * VanderWeele (2015) VanderWeele, T. (2015) _Explanation in Causal Inference: Methods for Mediation and Interaction_.  New York, NY: Oxford University Press.
  * VanderWeele and Vansteelandt (2014) VanderWeele, T. and Vansteelandt, S. (2014) Mediation analysis with multiple mediators.  _Epidemiologic Methods_ , 2, 95–115.
  * VanderWeele (2009) VanderWeele, T. J. (2009) Concerning the consistency assumption in causal inference.  _Epidemiology_ , 20, 880–883.
  * VanderWeele _et al._ (2014) VanderWeele, T. J., Vansteelandt, S. and Robins, J. M. (2014) Effect decomposition in the presence of an exposure-induced mediator-outcome confounder.  _Epidemiology_ , 25, 300–306.
  * Vansteelandt _et al._ (2012) Vansteelandt, S., Bekaert, M. and Lange, T. (2012) Imputation strategies for the estimation of natural direct and indirect effects.  _Epidemiologic Methods_ , 1, 131–158.
  * Vansteelandt and Daniel (2017) Vansteelandt, S. and Daniel, R. M. (2017) Interventional effects for mediation analysis with multiple mediators.  _Epidemiology (Cambridge, Mass.)_ , 28, 258.
  * Wodtke and Zhou (2020) Wodtke, G. T. and Zhou, X. (2020) Effect decomposition in the presence of treatment-induced confounding: A regression-with-residuals approach.  _Epidemiology_ , 31, 369–375.
  * Zheng and van der Laan (2011) Zheng, W. and van der Laan, M. J. (2011) Cross-validated targeted minimum-loss-based estimation.  In _Targeted Learning_ , 459–474. New York, NY: Springer.
  * Zheng and van der Laan (2012) Zheng, W. and van der Laan, M. J. (2012) Targeted maximum likelihood estimation of natural direct effects.  _The International Journal of Biostatistics_ , 8.
  * Zhou and Yamamoto (2020) Zhou, X. and Yamamoto, T. (2020) Tracing causal paths from experimental and observational data.  _SocArXiv preprint doi:10.31235/osf.io/2rx6p_.

## Supplementary Materials

###  Appendix A Proof of Theorem 1

Assumption 2* implies that for any $k\in\\{2,\ldots K\\}$ and any
$j\in\\{1,\ldots k-1\\}$,

$\displaystyle\big{(}M_{k}(a_{k},\overline{m}_{k-1}),\ldots
M_{K}(a_{K},\overline{m}_{K-1}),Y(a_{K+1},\overline{m}_{K})\big{)}\perp\\!\\!\\!\perp
M_{k-j}(a_{k-j},\overline{m}_{k-j-1}^{*})|X,A,\overline{M}_{k-j-1}$
$\displaystyle\Rightarrow$
$\displaystyle\big{(}M_{k}(a_{k},\overline{m}_{k-1}),\ldots
M_{K}(a_{K},\overline{m}_{K-1}),Y(a_{K+1},\overline{m}_{K})\big{)}\perp\\!\\!\\!\perp
M_{k-j}(a_{k-j},\overline{m}_{k-j-1}^{*})|X,A=a_{k-j},\overline{M}_{k-j-1}=\overline{m}_{k-j-1}^{*}$
$\displaystyle\Rightarrow$
$\displaystyle\big{(}M_{k}(a_{k},\overline{m}_{k-1}),\ldots
M_{K}(a_{K},\overline{m}_{K-1}),Y(a_{K+1},\overline{m}_{K})\big{)}\perp\\!\\!\\!\perp
M_{k-j}|X,A=a_{k-j},\overline{M}_{k-j-1}=\overline{m}_{k-j-1}^{*}$
$\displaystyle\Rightarrow$
$\displaystyle\big{(}M_{k}(a_{k},\overline{m}_{k-1}),\ldots
M_{K}(a_{K},\overline{m}_{K-1}),Y(a_{K+1},\overline{m}_{K})\big{)}\perp\\!\\!\\!\perp
M_{k-j}|X,A,\overline{M}_{k-j-1}.$ (20)

Setting $\overline{m}_{k-1}^{*}=\overline{m}_{k-1}$, Assumption 2* also
implies that for any $k\in[K]$,

$\big{(}M_{k+1}(a_{k+1},\overline{m}_{k}),\ldots
M_{K}(a_{K},\overline{m}_{K-1}),Y(a_{K+1},\overline{m}_{K})\big{)}\perp\\!\\!\\!\perp
M_{k}(a_{k},\overline{m}_{k-1})|X,A,\overline{M}_{k-1}.$ (21)

Now suppose that for some $j\in\\{1,\ldots k-1\\}$,

$\displaystyle\big{(}M_{k+1}(a_{k+1},\overline{m}_{k}),\ldots
M_{K}(a_{K},\overline{m}_{K-1}),Y(a_{K+1},\overline{m}_{K})\big{)}$
$\displaystyle\perp\\!\\!\\!\perp
M_{k}(a_{k},\overline{m}_{k-1})|X,A,\overline{M}_{k-j}.$ (22)

By the contraction rule of conditional independence, the relationships (20)
and (22) imply

$\big{(}M_{k+1}(a_{k+1},\overline{m}_{k}),\ldots
M_{K}(a_{K},\overline{m}_{K-1}),Y(a_{K+1},\overline{m}_{K})\big{)}\perp\\!\\!\\!\perp
M_{k}(a_{k},\overline{m}_{k-1})|X,A,\overline{M}_{k-j-1}.$

Hence, by the initial relationship (21) and mathematical induction, we have

$\big{(}M_{k+1}(a_{k+1},\overline{m}_{k}),\ldots
M_{K}(a_{K},\overline{m}_{K-1}),Y(a_{K+1},\overline{m}_{K})\big{)}\perp\\!\\!\\!\perp
M_{k}(a_{k},\overline{m}_{k-1})|X,A,\quad\forall k\in[K].$ (23)

In the meantime, because $\big{(}M_{k+1}(a_{k+1},\overline{m}_{k}),\ldots
M_{K}(a_{K},\overline{m}_{K-1}),Y(a_{K+1},\overline{m}_{K})\big{)}\perp\\!\\!\\!\perp
A|X$, we have (by the contraction rule)

$\big{(}M_{k+1}(a_{k+1},\overline{m}_{k}),\ldots
M_{K}(a_{K},\overline{m}_{K-1}),Y(a_{K+1},\overline{m}_{K})\big{)}\perp\\!\\!\\!\perp\big{(}A,M_{k}(a_{k},\overline{m}_{k-1})\big{)}|X,\quad\forall
k\in[K].$

Thus the components in $\big{(}A,M_{1}(a_{1}),\ldots
M_{K}(a_{K},\overline{m}_{K-1}),Y(a_{K+1},\overline{m}_{K})\big{)}$ are
mutually independent given $X$. Therefore,

$\displaystyle\psi_{\overline{a}}$
$\displaystyle=\mathbb{E}[Y(a_{K+1},\overline{M}_{K}(\overline{a}_{K}))]$
$\displaystyle=\int_{x}\int_{\overline{m}_{K}}\mathbb{E}[Y(a_{K+1},\overline{m}_{K})|X=x,A=a_{K+1},M_{1}(a_{1})=m_{1},\ldots
M_{K}(a_{K},\overline{m}_{K-1})=m_{K}]$
$\displaystyle\big{(}\prod_{k=1}^{K}dP_{M_{k}(a_{k},\overline{m}_{k-1})|X,A,M_{1}(a_{1}),\ldots
M_{k-1}(a_{k-1},\overline{m}_{k-2})}(m_{k}|x,a_{K+1},\overline{m}_{k-1})\big{)}dP_{X}(x)$
$\displaystyle=\int_{x}\int_{\overline{m}_{K}}\mathbb{E}[Y(a_{K+1},\overline{m}_{K})|X=x,A=a_{K+1}]\big{(}\prod_{k=1}^{K}dP_{M_{k}(a_{k},\overline{m}_{k-1})|X}(m_{k}|x)\big{)}dP_{X}(x)$
$\displaystyle=\int_{x}\int_{\overline{m}_{K}}\mathbb{E}[Y(a_{K+1},\overline{m}_{K})|x,a_{K+1},\overline{M}_{K}=\overline{m}_{K}]\big{(}\prod_{k=1}^{K}dP_{M_{k}(a_{k},\overline{m}_{k-1})|X,A,\overline{M}_{k-1}}(m_{k}|x,a_{k},\overline{m}_{k-1})\big{)}dP_{X}(x)$
$\displaystyle=\int_{x}\int_{\overline{m}_{K}}\mathbb{E}[Y|x,a_{K+1},\overline{m}_{K}]\big{(}\prod_{k=1}^{K}dP_{M_{k}|X,A,\overline{M}_{k-1}}(m_{k}|x,a_{k},\overline{m}_{k-1})\big{)}dP_{X}(x).$

###  Appendix B Hybrid Estimators of $\theta_{a_{1},a_{2},a}$

For notational brevity, let us use the following shorthands:

$\displaystyle\lambda_{0}^{j}(A|X)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\mathbb{I}(A=a_{j})}{p(a_{j}|X)}$
$\displaystyle\lambda_{1}^{j}(M_{1}|X)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}\frac{p(M_{1}|X,a_{1})}{p(M_{1}|X,a_{j})}$
$\displaystyle\lambda_{2}^{j}(M_{2}|X,M_{1})$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}\frac{p(M_{2}|X,a_{2},M_{1})}{p(M_{2}|X,a_{j},M_{1})},$

In addition, define $\lambda_{0}(A|X)=\mathbb{I}(A=a)/p(a|X)$,
$\lambda_{1}(M_{1}|X)=p(M_{1}|X,a_{1})/p(M_{1}|X,a)$ and
$\lambda_{2}(M_{2}|X,M_{1})=p(M_{2}|X,a_{2},M_{1})/p(M_{2}|X,a,M_{1})$. With
the above notation, the iterated conditional means $\mu_{1}(X,M_{1})$,
$\mu_{0}(X)$, and $\theta_{a_{1},a_{2},a}$ can each be written in several
different forms:

$\displaystyle\mu_{1}(X,M_{1})$
$\displaystyle=\begin{cases}\mathbb{E}[\mu_{2}(X,M_{1},M_{2})|X,a_{2},M_{1}]\\\
\mathbb{E}[\lambda_{2}(M_{2}|X,M_{1})Y|X,a,M_{1}]\end{cases}$
$\displaystyle\mu_{0}(X)$
$\displaystyle=\begin{cases}\mathbb{E}[\mu_{1}(X,M_{1})|X,a_{1}]=\begin{cases}\mathbb{E}\big{[}\mathbb{E}[\mu_{2}(X,M_{1},M_{2})|X,a_{2},M_{1}]|X,a_{1}\big{]}\\\
\mathbb{E}\big{[}\mathbb{E}[\lambda_{2}(M_{2}|X,M_{1})Y|X,a,M_{1}]|X,a_{1}\big{]}\end{cases}\\\
\mathbb{E}[\lambda_{1}^{2}(M_{1}|X)\mu_{2}(X,M_{1},M_{2})|X,a_{2}]\\\
\mathbb{E}[\lambda_{1}(M_{1}|X)\lambda_{2}(X,M_{1},M_{2})Y|X,a]\end{cases}$
$\displaystyle\theta_{a_{1},a_{2},a}$
$\displaystyle=\begin{cases}\mathbb{E}[\mu_{0}(X)]=\begin{cases}\mathbb{E}\Big{[}\mathbb{E}\big{[}\mathbb{E}[\mu_{2}(X,M_{1},M_{2})|X,a_{2},M_{1}]|X,a_{1}\big{]}\Big{]}&\textup{(RI-
RI-RI)}\\\
\mathbb{E}\Big{[}\mathbb{E}\big{[}\mathbb{E}[\lambda_{2}(M_{2}|X,M_{1})Y|X,a,M_{1}]|X,a_{1}\big{]}\Big{]}&\textup{(W-RI-
RI)}\\\
\mathbb{E}\Big{[}\mathbb{E}[\lambda_{1}^{2}(M_{1}|X)\mu_{2}(X,M_{1},M_{2})|X,a_{2}]\Big{]}&\textup{(RI-
W-RI)}\\\
\mathbb{E}\Big{[}\mathbb{E}[\lambda_{1}(M_{1}|X)\lambda_{2}(X,M_{1},M_{2})Y|X,a]\Big{]}&\textup{(W-W-
RI)}\end{cases}\\\
\mathbb{E}[\lambda_{0}^{1}(A|X)\mu_{1}(X,M_{1})]=\begin{cases}\mathbb{E}\big{[}\lambda_{0}^{1}(A|X)\mathbb{E}[\mu_{2}(X,M_{1},M_{2})|X,a_{2},M_{1}]\big{]}&\textup{(RI-
RI-W)}\\\
\mathbb{E}\big{[}\lambda_{0}^{1}(A|X)\mathbb{E}[\lambda_{2}(M_{2}|X,M_{1})Y|X,a,M_{1}]\big{]}&\textup{(W-RI-W)}\end{cases}\\\
\mathbb{E}[\lambda_{0}^{2}(A|X)\lambda_{1}^{2}(M_{1}|X)\mu_{2}(X,M_{1},M_{2})]\quad\;\;\textup{(RI-
W-W)}\\\
\mathbb{E}[\lambda_{0}(A|X)\lambda_{1}(M_{1}|X)\lambda_{2}(M_{2}|X,M_{1})Y]\quad\textup{(W-W-W)}\end{cases}$

The first set of equations suggest two different ways of estimating
$\mu_{1}(x,m_{1})$: (a) fit a model for the conditional mean of
$\hat{\mu}_{2}(X,M_{1},M_{2})$ given $X$, $A$, $M_{1}$ and then set $A=a_{2}$
for all units; (b) fit a model for the conditional mean of
$\hat{\lambda}_{2}(M_{2}|X,M_{1})Y$ given $X$, $A$, and $M_{1}$ and then set
$A=a$ for all units. Similarly, the second set of equations suggest four
different ways of estimating $\mu_{0}(x)$, and the last set of equations point
to eight different ways of estimating $\theta_{a_{1},a_{2},a}$. Each of these
eight estimators corresponds to a unique combination of regression-imputation
and weighting.

###  Appendix C Proof of Theorem 2

To show that equation (11) is the EIF of $\theta_{\overline{a}}$ in
$\mathcal{P}_{\textup{np}}$, it suffices to show

$\frac{\partial\theta_{\overline{a}}(t)}{\partial
t}\biggl{|}_{t=0}=\mathbb{E}[\varphi_{\overline{a}}(O)S_{0}(O)],$ (24)

where $S_{0}(O)$ is the score function for any one-dimensional submodel
$P_{t}(O)$ evaluated at $t=0$. We first note that $S_{t}(O)$ can be written as
$S_{t}(O)=S_{t}(X)+S_{t}(A|X)+\sum_{k=1}^{K}S_{t}(M_{k}|X,A,\overline{M}_{k-1})+S_{t}(Y|X,A,\overline{M}_{K})$,
where $S_{t}(u|v)=\partial\log p_{t}(u|v)/\partial t$ and $p_{t}(u|v)$ is the
conditional probability density/mass function of $U$ given $V$. Using equation
(2) and the product rule, the left-hand side of equation (24) can be written
as

$\displaystyle\frac{\partial\theta_{\overline{a}}(t)}{\partial
t}\biggl{|}_{t=0}$ $\displaystyle=\frac{\partial\iiint
ydP_{t}(y|x,a_{K+1},\overline{m}_{K})\big{[}\prod_{k=1}^{K}dP_{t}(m_{k}|x,a_{k},\overline{m}_{k-1})\big{]}dP_{t}(x)}{\partial
t}\biggl{|}_{t=0}$ $\displaystyle=\underbrace{\iiint
yS_{0}(x)dP_{0}(y|x,a_{K+1},\overline{m}_{K})\big{[}\prod_{k=1}^{K}dP_{0}(m_{k}|x,a_{k},\overline{m}_{k-1})\big{]}dP_{0}(x)}_{=:\textup{$\phi_{0}$}}$
$\displaystyle+\sum_{k=1}^{K}\underbrace{\iiint
yS_{0}(m_{k}|x,a_{k},\overline{m}_{k-1})dP_{0}(y|x,a_{K+1},\overline{m}_{K})\big{[}\prod_{k=1}^{K}dP_{0}(m_{k}|x,a_{k},\overline{m}_{k-1})\big{]}dP_{0}(x)}_{=:\phi_{k}}$
$\displaystyle+\underbrace{\iiint
yS_{0}(y|x,a_{K+1},\overline{m}_{K})dP_{0}(y|x,a_{K+1},\overline{m}_{K})\big{[}\prod_{k=1}^{K}dP_{0}(m_{k}|x,a_{k},\overline{m}_{k-1})\big{]}dP_{0}(x)}_{=:\phi_{K+1}}$
$\displaystyle=\sum_{k=0}^{K+1}\phi_{k}$

where the second equality follows from the fact that $\partial
dP_{t}(u|v)/\partial t=S_{t}(u|v)dP_{t}(u|v).$ Below, we verify that
$\phi_{k}=\mathbb{E}[\varphi_{k}(O)S_{0}(O)]$ for all $k\in\\{0,\ldots
K+1\\}$, where $\varphi_{k}(O)$ is defined in Theorem 2. First,

$\displaystyle\mathbb{E}[\varphi_{0}(O)S_{0}(O)]$ $\displaystyle=$
$\displaystyle\mathbb{E}[\big{(}\mu_{0}(X)-\theta_{\overline{a}}\big{)}S_{0}(O)]$
$\displaystyle=$ $\displaystyle\mathbb{E}[\mu_{0}(X)S_{0}(O)]$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\mu_{0}(X)\big{(}S_{0}(X)+S_{0}(A|X)+\sum_{k=1}^{K}S_{0}(M_{k}|X,A,\overline{M}_{k-1})+S_{0}(Y|X,A,\overline{M}_{K})\big{)}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\mu_{0}(X)S_{0}(X)\big{]}+\mathbb{E}\Big{[}\mu_{0}(X)\underbrace{\mathbb{E}\big{[}S_{0}(A|X)|X\big{]}}_{=0}\Big{]}+\sum_{k=1}^{K}\mathbb{E}\Big{[}\mu_{0}(X)\underbrace{\mathbb{E}\big{[}S_{0}(M_{k}|X,A,\overline{M}_{k-1})|X,A,\overline{M}_{k-1}\big{]}}_{=0}\Big{]}$
$\displaystyle+\mathbb{E}\Big{[}\mu_{0}(X)\underbrace{\mathbb{E}\big{[}S_{0}(Y|X,A,\overline{M}_{K})|X,A,\overline{M}_{K}\big{]}}_{=0}\Big{]}$
$\displaystyle=$ $\displaystyle\int\mu_{0}(x)S_{0}(x)dP_{0}(x)$
$\displaystyle=$ $\displaystyle\iiint
yS_{0}(x)dP_{0}(y|x,a_{K+1},\overline{m}_{K})\big{[}\prod_{k=1}^{K}dP_{0}(m_{k}|x,a_{k},\overline{m}_{k-1})\big{]}dP_{0}(x)$
$\displaystyle=$ $\displaystyle\phi_{0}.$

Second, for $k\in[K]$,

$\displaystyle\mathbb{E}[\varphi_{k}(O)S_{0}(O)]$ $\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\varphi_{k}(O)\big{(}S_{0}(X)+S_{0}(A|X)+\sum_{j=1}^{K}S_{0}(M_{j}|X,A,\overline{M}_{j-1})+S_{0}(Y|X,A,\overline{M}_{K})\big{)}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\mathbb{E}\big{[}\varphi_{k}(O)\big{(}S_{0}(X)+S_{0}(A|X)+\sum_{j=1}^{k-1}S_{0}(M_{j}|X,A,\overline{M}_{j-1})\big{)}|X,A,\overline{M}_{k-1}\big{]}\Big{]}+\mathbb{E}\big{[}\varphi_{k}(O)S_{0}(M_{k}|X,A,\overline{M}_{k-1})\big{)}\big{]}$
$\displaystyle+\sum_{j=k+1}^{K}\mathbb{E}\Big{[}\varphi_{k}(O)\underbrace{\mathbb{E}\big{[}S_{0}(M_{j}|X,A,\overline{M}_{j-1})|X,A,\overline{M}_{j-1}\big{]}}_{=0}\Big{]}+\mathbb{E}\Big{[}\varphi_{k}(O)\underbrace{\mathbb{E}\big{[}S_{0}(Y|X,A,\overline{M}_{K})|X,A,\overline{M}_{K}\big{]}}_{=0}\Big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\big{(}S_{0}(X)+S_{0}(A|X)+\sum_{j=1}^{k-1}S_{0}(M_{j}|X,A,\overline{M}_{j-1})\big{)}\underbrace{\mathbb{E}\big{[}\varphi_{k}(O)|X,A,\overline{M}_{k-1}\big{]}}_{=0}\Big{]}+\mathbb{E}\big{[}\varphi_{k}(O)S_{0}(M_{k}|X,A,\overline{M}_{k-1})\big{)}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\varphi_{k}(O)S_{0}(M_{k}|X,A,\overline{M}_{k-1})\big{)}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\mathbb{E}\big{[}\frac{\mathbb{I}(A=a_{k})}{p(a_{k}|X)}\Big{(}\prod_{j=1}^{k-1}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{k},\overline{M}_{j-1})}\Big{)}\big{(}\mu_{k}(X,\overline{M}_{k})-\mu_{k-1}(X,\overline{M}_{k-1})\big{)}S_{0}(M_{k}|X,A,\overline{M}_{k-1})\big{)}\big{|}X,A,\overline{M}_{k-1}\big{]}\Big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\frac{\mathbb{I}(A=a_{k})}{p(a_{k}|X)}\Big{(}\prod_{j=1}^{k-1}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{k},\overline{M}_{j-1})}\Big{)}\mu_{k}(X,\overline{M}_{k})S_{0}(M_{k}|X,A,\overline{M}_{k-1})\big{)}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\mathbb{E}\big{[}\Big{(}\prod_{j=1}^{k-1}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{k},\overline{M}_{j-1})}\Big{)}\mu_{k}(X,\overline{M}_{k})S_{0}(M_{k}|X,A,\overline{M}_{k-1})\big{)}|X,A=a_{k}\big{]}$
$\displaystyle=$ $\displaystyle\iiint
S_{0}(m_{k}|x,a_{k},\overline{m}_{k-1})\Big{(}\int_{y}\int_{\overline{m}_{K}}ydP_{0}(y|x,a_{K+1},\overline{m}_{K})\prod_{j=k+1}^{K}dP_{0}(m_{j}|x,a_{j},\overline{m}_{j-1})\Big{)}$
$\displaystyle\cdot$ $\displaystyle
dP_{0}(m_{k}|x,a_{k},\overline{m}_{k-1})\Big{(}\prod_{j=1}^{k-1}\frac{p(m_{j}|x,a_{j},\overline{m}_{j-1})}{p(m_{j}|x,a_{k},\overline{m}_{j-1})}\Big{)}\Big{(}\prod_{j=1}^{k-1}dP_{0}(m_{j}|x,a_{k},\overline{m}_{j-1})\Big{)}dP_{0}(x)$
$\displaystyle=$ $\displaystyle\iiint
yS_{0}(m_{k}|x,a_{k},\overline{m}_{k-1})dP_{0}(y|x,a_{K+1},\overline{m}_{K})\big{(}\prod_{j=1}^{K}dP_{0}(m_{j}|x,a_{j},\overline{m}_{j-1})\big{)}dP_{0}(x)$
$\displaystyle=$ $\displaystyle\phi_{k},$

where the fourth equality is due to the fact that

$\displaystyle\mathbb{E}\big{[}\varphi_{k}(O)|X,A,\overline{M}_{k-1}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\frac{\mathbb{I}(A=a_{k})}{p(a_{k}|X)}\Big{(}\prod_{j=1}^{k-1}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{k},\overline{M}_{j-1})}\Big{)}\big{(}\mu_{k}(X,\overline{M}_{k})-\mu_{k-1}(X,\overline{M}_{k-1})\big{)}|X,A,\overline{M}_{k-1}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\Big{(}\prod_{j=1}^{k-1}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{k},\overline{M}_{j-1})}\Big{)}\big{(}\mu_{k}(X,\overline{M}_{k})-\mu_{k-1}(X,\overline{M}_{k-1})\big{)}|X,A=a_{k},\overline{M}_{k-1}\big{]}$
$\displaystyle=$
$\displaystyle\Big{(}\prod_{j=1}^{k-1}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{k},\overline{M}_{j-1})}\Big{)}\underbrace{\mathbb{E}\big{[}\mu_{k}(X,\overline{M}_{k})-\mu_{k-1}(X,\overline{M}_{k-1})|X,A=a_{k},\overline{M}_{k-1}\big{]}}_{=0}$
$\displaystyle=$ $\displaystyle 0.$

Finally,

$\displaystyle\mathbb{E}[\varphi_{K+1}(O)S_{0}(O)]$ $\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\varphi_{K+1}(O)\big{(}S_{0}(X)+S_{0}(A|X)+\sum_{j=1}^{K}S_{0}(M_{j}|X,A,\overline{M}_{j-1})+S_{0}(Y|X,A,\overline{M}_{K})\big{)}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\mathbb{E}\big{[}\varphi_{K+1}(O)\big{(}S_{0}(X)+S_{0}(A|X)+\sum_{j=1}^{K}S_{0}(M_{j}|X,A,\overline{M}_{j-1})\big{)}|X,A,\overline{M}_{K}\big{]}\Big{]}+\mathbb{E}\big{[}\varphi_{K+1}(O)S_{0}(Y|X,A,\overline{M}_{K})\big{)}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\big{(}S_{0}(X)+S_{0}(A|X)+\sum_{j=1}^{K}S_{0}(M_{j}|X,A,\overline{M}_{j-1})\big{)}\underbrace{\mathbb{E}\big{[}\varphi_{K+1}(O)|X,A,\overline{M}_{K}\big{]}}_{=0}\Big{]}+\mathbb{E}\big{[}\varphi_{K+1}(O)S_{0}(Y|X,A,\overline{M}_{K})\big{)}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\varphi_{K+1}(O)S_{0}(Y|X,A,\overline{M}_{K})\big{)}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\mathbb{E}\big{[}\frac{\mathbb{I}(A=a_{K+1})}{p(a_{K+1}|X)}\Big{(}\prod_{j=1}^{K}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{K+1},\overline{M}_{j-1})}\Big{)}\big{(}Y-\mu_{K}(X,\overline{M}_{K})\big{)}S_{0}(Y|X,A,\overline{M}_{K})\big{)}\big{|}X,A,\overline{M}_{K}\big{]}\Big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\mathbb{E}\big{[}\frac{\mathbb{I}(A=a_{K+1})}{p(a_{K+1}|X)}\Big{(}\prod_{j=1}^{K}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{K+1},\overline{M}_{j-1})}\Big{)}YS_{0}(Y|X,A,\overline{M}_{K})\big{)}\big{|}X,A,\overline{M}_{K}\big{]}\Big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\mathbb{E}\big{[}\frac{\mathbb{I}(A=a_{K+1})}{p(a_{K+1}|X)}\Big{(}\prod_{j=1}^{K}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{K+1},\overline{M}_{j-1})}\Big{)}YS_{0}(Y|X,A,\overline{M}_{K})\big{)}\big{|}X,A\big{]}\Big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\Big{[}\mathbb{E}\big{[}\Big{(}\prod_{j=1}^{K}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{K+1},\overline{M}_{j-1})}\Big{)}YS_{0}(Y|X,A,\overline{M}_{K})\big{)}\big{|}X,A=a_{K+1}\big{]}\Big{]}$
$\displaystyle=$ $\displaystyle\iiint
yS_{0}(y|x,a_{K+1},\overline{m}_{K})dP_{0}(y|x,a_{K+1},\overline{m}_{K})\Big{(}\prod_{j=1}^{K}\frac{p(m_{j}|x,a_{j},\overline{m}_{j-1})}{p(m_{j}|x,a_{K+1},\overline{m}_{j-1})}\Big{)}\big{(}\prod_{j=1}^{K}dP_{0}(m_{j}|x,a_{K+1},\overline{m}_{j-1})\big{)}dP_{0}(x)$
$\displaystyle=$ $\displaystyle\iiint
yS_{0}(y|x,a_{K+1},\overline{m}_{K})dP_{0}(y|x,a_{K+1},\overline{m}_{K})\big{(}\prod_{j=1}^{K}dP_{0}(m_{j}|x,a_{j},\overline{m}_{j-1})\big{)}dP_{0}(x)$
$\displaystyle=$ $\displaystyle\phi_{K+1},$

where the third equality is due to the fact that

$\displaystyle\mathbb{E}\big{[}\varphi_{K+1}(O)|X,A,\overline{M}_{K}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\frac{\mathbb{I}(A=a_{K+1})}{p(a_{K+1}|X)}\Big{(}\prod_{j=1}^{K}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{K+1},\overline{M}_{j-1})}\Big{)}\big{(}Y-\mu_{K}(X,\overline{M}_{K})\big{)}|X,A,\overline{M}_{K}\big{]}$
$\displaystyle=$
$\displaystyle\frac{p(a_{K+1}|X,\overline{M}_{K})}{p(a_{K+1}|X)}\Big{(}\prod_{j=1}^{K}\frac{p(M_{j}|X,a_{j},\overline{M}_{j-1})}{p(M_{j}|X,a_{K+1},\overline{M}_{j-1})}\Big{)}\underbrace{\mathbb{E}\big{[}Y-\mu_{K}(X,\overline{M}_{K})|X,A=a_{K+1},\overline{M}_{K}\big{]}}_{=0}$
$\displaystyle=$ $\displaystyle 0.$