# Causal mediation analysis: From simple to more robust strategies for
estimation of marginal natural (in)direct effects

Trang Quynh Nguyen, Elizabeth L. Ogburn, Elizabeth B. Sarker, Noah Greifer,
Ian Schmid, Ina M. Koning, and Elizabeth A. Stuart Affiliations: TQN
<EMAIL_ADDRESS>is Assistant Scientist, Department of Mental Health
(DMH), Johns Hopkins Bloomberg School of Public Health (JHSPH), USA. ELO is
Associate Professor, Department of Biostatistics (DB), JHSPH. EBS is Doctoral
Candidate, DB, JHSPH. NG is Postdoctoral Researcher, DMH, JHSPH. IS is
Doctoral Candidate, DMH, JHSPH. IMK is Assistant Professor, Social and
Behavioral Sciences, University of Utrecht, The Netherlands. EAS is Professor,
DMH, DB and Department of Health Policy and Management,
JHSPH.Acknowledgements: This work is funded by NIMH grant R01MH115487 (PI
Stuart). The authors appreciate helpful feedback from Drs. Guanglei Hong and
Fan Yang on an earlier draft. TQN thanks Drs. Ilya Shpitser and Eric Tchetgen
Tchetgen for insightful discussions on robust estimation and their seminal
article Tchetgen Tchetgen and Shpitser [2012].

###### Abstract

This paper aims to provide practitioners of causal mediation analysis with a
better understanding of estimation options. We take as inputs two familiar
strategies (weighting and model-based prediction) and a simple way of
combining them (weighted models), and show how we can generate a range of
estimators with different modeling requirements and robustness properties. The
primary goal is to help build intuitive appreciation for robust estimation
that is conducive to sound practice. A second goal is to provide a “menu” of
estimators that practitioners can choose from for the estimation of marginal
natural (in)direct effects. The estimators generated from this exercise
include some that coincide or are similar to existing estimators and others
that have not previously appeared in the literature. We note several different
ways to estimate the weights for cross-world weighting based on three
expressions of the weighting function, including one that is novel; and show
how to check the resulting covariate and mediator balance. We use a random
continuous weights bootstrap to obtain confidence intervals, and also derive
general asymptotic (sandwich) variance formulas for the estimators. The
estimators are illustrated using data from an adolescent alcohol use
prevention study.

##  1 Introduction

Causal mediation methodology is complex. There are different types of causal
contrasts: controlled direct effect, natural (in)direct effects [Robins and
Greenland, 1992, Pearl, 2001], interventional (in)direct [Didelez et al.,
2006, VanderWeele et al., 2014] and other interventional effects [Nguyen et
al., 2020b], etc., each with their own set of identification assumptions
[Nguyen et al., 2020a]. The literature on effect estimation is vast, with a
wide variety of estimation methods based on regression [e.g., VanderWeele and
Vansteelandt, 2010, Valeri and VanderWeele, 2013, Muthén and Asparouhov,
2015], weighting [e.g., Hong, 2010, Hong et al., 2015, Tchetgen Tchetgen,
2013, Huber, 2014, VanderWeele and Vansteelandt, 2013], simulation [e.g., Imai
et al., 2010b, a, Vansteelandt et al., 2012, VanderWeele and Vansteelandt,
2013], or some combination of these strategies. To add to the complexity of
the picture, some methods estimate marginal effects [Imai et al., 2010a, Hong,
2010, Hong et al., 2015, Lange et al., 2012, VanderWeele and Vansteelandt,
2013, Tchetgen Tchetgen and Shpitser, 2012], others estimate effects
conditional on covariates [VanderWeele and Vansteelandt, 2010, Valeri and
VanderWeele, 2013, Steen et al., 2017, Tchetgen Tchetgen and Shpitser, 2014,
Tchetgen Tchetgen, 2013]. Most of these methods are parametric and require all
the models used to be correctly specified. Some methods have built in
robustness to model misspecification; these are often presented in highly
technical papers [e.g., Tchetgen Tchetgen and Shpitser, 2012, Zheng and van
der Laan, 2012].

This paper explicates a range of estimation options for causal mediation,
focusing on options with some robustness properties. We take a concrete
approach of using as inputs two familiar strategies (weighting and regression)
and a simple way of combining them, and show how to generate a range of
estimators with different modeling requirements and robustness properties. The
primary goal is to help build intuitive appreciation for robust estimation
that is conducive to sound practice (without requiring prior understanding of
these methods). A secondary goal is to provide a “menu” of estimators that
practitioners can choose from (depending on which modeling components they
feel comfortable with given their specific problem).

The paper focuses on natural (in)direct effects. These decompose the total
causal effect and (when identified) provide insight about effect mechanisms.
This is a direct match to researchers’ common motivation for conducting
mediation analysis – a wish to understand what part of a causal effect is
indirect (operating through a specific intermediate variable) and what part is
direct (not through that variable). The kind of reasoning used to build
estimators here is not specific to these effects; they can be applied to
effects of other types as well. We consider marginal natural (in)direct
effects, which, defined on the additive scale, correspond to the total effect
being the average treatment effect (a popular effect in causal inference).
Adaptation to average effects on the treated or on the controls is trivial.
This paper does not address the estimation of conditional effects as functions
of covariates, which entails a different set of estimation strategies that
should be tackled separately.

Our construction of estimators is a bottom-up exercise and is relatively
exhaustive, so not all the estimators generated have appeared in the
literature. We connect to work that employs, or is related to, the strategies
and estimators discussed in this paper, and comment on the differences (some
quite subtle) between some of these estimators.

To make the paper accessible to a broad audience including not only
statisticians but also methodologically savvy substantive researchers, all the
proofs (about robustness properties and large-sample variance) are placed in
the Technical Appendix. To make things concrete and facilitate application,
the Computing Appendix explains step-by-step implementation on real data, with
adaptable R code.

##  2 Preliminaries

###  2.1 Effect definitions

We consider the setting with a binary exposure $A$, followed in time by a
mediator variable $M$ (which may be multivariate), followed in time by an
outcome $Y$. Our target estimands are marginal natural (in)direct effects,
which decompose the marginal total effect. We focus on effects defined on the
additive scale, but will comment on multiplicative scales.

These effects are defined using potential outcomes [Rubin, 1974, Holland,
1986]. The marginal additive total effect is formally defined as
$\mathrm{TE}:=\text{E}[Y_{1}]-\text{E}[Y_{0}]$, the difference between the
population mean of $Y_{1}$ (potential outcome if exposed to the active
treatment) and that of $Y_{0}$ (potential outcome if exposed to the comparison
condition). Definition of natural (in)direct effects [Pearl, 2001]
additionally employs a nested potential outcome type, $Y_{aM_{a^{\prime}}}$
(for a hypothetical condition where exposure had been set to $a$ and mediator
had been set to its potential value under condition $a^{\prime}$) where $a$
and $a^{\prime}$ can be either 0 or 1. We assume that $Y_{a}=Y_{aM_{a}}$, thus
$\mathrm{TE}=\text{E}[Y_{1M_{1}}]-\text{E}[Y_{0M_{0}}]$. Using a third
potential outcome with mismatched $a$ and $a^{\prime}$, either $Y_{1M_{0}}$
(had exposure been set to the active treatment but mediator set to its
potential value under control) or $Y_{0M_{1}}$ (the other way around), TE is
decomposed in two ways, giving rise to two pairs of natural (in)direct
effects:

$\text{TE}=\underbrace{\text{E}[Y_{1M_{1}}]-\text{E}[Y_{1M_{0}}]}_{\textstyle\text{NIE}_{1}}\leavevmode\nobreak\
+\leavevmode\nobreak\
\underbrace{\text{E}[Y_{1M_{0}}]-\text{E}[Y_{0M_{0}}]}_{\textstyle\text{NDE}_{0}},$
$\text{TE}=\underbrace{\text{E}[Y_{1M_{1}}]-\text{E}[Y_{0M_{1}}]}_{\textstyle\text{NDE}_{1}}\leavevmode\nobreak\
+\leavevmode\nobreak\
\underbrace{\text{E}[Y_{0M_{1}}]-\text{E}[Y_{0M_{0}}]}_{\textstyle\text{NIE}_{0}}.$

The marginal TE defined above is equal to the mean of the causal effect on the
individual, $Y_{1}−Y_{0}$, and thus is usually known as the average treatment
effect in the non-mediation literature. The natural (in)direct effects defined
above can also be seen as averages of effects on individuals. This
marginal/average connection holds for effects defined on the additive (but not
multiplicative) scale.

Note that each decomposition here includes only one indirect effect. In a
situation where $M$ is a set of more than one mediator (as in our data
example), this is the effect mediated by all the mediators combined.
Alternatively, one may be interested in path-specific effects involving
different mediators or subsets of mediators; that is a more complex problem
that is outside the scope of this paper.

For conciseness, the rest of the paper addresses one of the two effect pairs:
the $\text{NDE}_{0}$ and $\text{NIE}_{1}$ (also called the pure direct effect
and total indirect effect [Robins and Greenland, 1992]). The other effect pair
mirrors this one in all content covered here.

###  2.2 Assumptions for effect identification

As the current focus is estimation, we simply assume that the effects of
interest are identified, noting that this is a matter for careful judgment in
applications. By “identified” we mean that the effects, which are functions of
potential outcomes, can be equated (under certain assumptions) to some
functions of the observed data distribution which we then attempt to estimate.
Below are the assumptions we make for ($\text{NDE}_{0}$, $\text{NIE}_{1}$)
identification; for more detailed explication, see VanderWeele and
Vansteelandt [2009], Imai et al. [2010b], Pearl [2012] or Nguyen et al.
[2020a].

The first assumption is that there is consistency between observed and
potential outcomes or mediator values, and between potential outcomes of
several types. Specifically,

$\displaystyle Y=Y_{a}\text{ if }A=a,$ $\displaystyle Y=Y_{1m}\text{ if
}A=1,M=m,$ $\displaystyle M=M_{0}\text{ if }A=0,$ $\displaystyle
Y_{a}=Y_{aM_{a}},$ $\displaystyle Y_{1M_{0}}=Y_{1m}\text{ if }A=0,M=m,$

for $a$ being either 0 or 1, and $m$ being any mediator value. Essentially we
have invoked this assumption in defining the effects above.

Figure 1: Unconfoundedness holds when no mediator-outcome confounder is
influenced by exposure and a set of observed pre-exposure covariates $C$
captures all confounding.$C$$A$$M$$Y$

The second assumption may be called ignorability, exchangeability or
unconfoundedness. This assumption requires that there is a set of observed
pre-exposure covariates $C$ (where “pre-exposure” means either preceding
exposure in time or simply not being influenced by exposure) that provides
several conditional independence relationships. Specifically,

$\displaystyle A\perp\\!\\!\\!\perp Y_{a},Y_{1m},M_{0}\mid C,$ $\displaystyle
M\perp\\!\\!\\!\perp Y_{1m}\mid C,A=1,$ $\displaystyle
M_{0}\perp\\!\\!\\!\perp Y_{1m}\mid C,$

for $a$ being either 0 or 1, and $m$ being any value in the distribution of
the mediator given covariates $C$ in the unexposed. The first two of the three
elements above fit with the usual notion of ignorability, where once we
condition on some variables the observed exposure (or mediator value) does not
carry any information about certain potential variables. The last element is
different in that it requires conditional independence of two potential
variables in two different worlds (thus commonly known as the cross-world
independence assumption).

In practice the usual way to deal with the unconfoundedness assumption is to
ask (i) whether there is any mediator-outcome confounder (observed or not)
that is influenced by exposure (this violates cross-world independence); and
if not, (ii) whether there is a set of pre-exposure covariates $C$ (all of
which observed) that captures all exposure-mediator, exposure-outcome and
mediator-outcome confounders. If either the answer to (i) is yes or the answer
to (ii) is no, then the unconfoundedness assumption does not hold. Note though
that while we could use substantive knowledge to judge the plausibility of
these assumptions, these assumptions are not testable using data.

Since identification involves conditioning on covariates $C$, what is also
required is that for all covariate levels there are positive chances of
observing relevant potential mediator/outcomes. This is the third assumption,
termed positivity. Specifically,

$\displaystyle\text{P}(A=a\mid C)>0,$ $\displaystyle\text{P}(M=m\mid
C,A=1)>0,$

for $a$ being either 0 or 1, and $m$ being any value in the distribution of
the mediator given covariates $C$ in the unexposed. The first element implies
positive chances of observing $Y_{1},Y_{0},M_{0}$; both combined imply
positive chances of observing $Y_{1m}$.

###  2.3 A heuristic view of identification that clarifies the estimation task

Identification of $(\text{NDE}_{0},\text{NIE}_{1})$ amounts to identifying the
means of the three potential outcomes $Y_{1},Y_{0}$ and $Y_{1M_{0}}$. Under
the assumptions above, the identification results [Pearl, 2001, VanderWeele
and Vansteelandt, 2009, Imai et al., 2010b, Pearl, 2012] of these three means
are

$\displaystyle\text{E}[Y_{1}]=\text{E}_{C}\\{\text{E}[Y\mid C,A=1]\\},$
$\displaystyle\text{E}[Y_{0}]=\text{E}_{C}\\{\text{E}[Y\mid C,A=0]\\},$
$\displaystyle\text{E}[Y_{1M_{0}}]=\text{E}_{C}(\text{E}_{M\mid
C,A=0}\\{\text{E}[Y\mid C,M,A=1]\\}),$

where the right-hand sides are functions of the observed data distribution.

To make these results more intuitive to readers who may find them unfamiliar,
we offer a heuristic visualization in Fig. 2. This figure has three columns.
The left column shows the data that we have: the full sample, which is
comprised of the treated subsample and the control subsample. As we are
interested in the $\text{NDE}_{0}$ and $\text{NIE}_{1}$ that together contrast
the means of $Y_{1},Y_{0},Y_{1M_{0}}$, what we would ideally like to have
instead is shown in the middle column: three full samples that all resemble
the actual full sample pre-exposure, but are then set to three conditions: the
treated (1) condition, the control (0) condition, and the cross-world
condition (where exposure is set to 1 but mediator is set to $M_{0}$); this
would allow us to average the outcome in the three samples to estimate the
three potential outcome means. Unfortunately, we do not observe these three
full samples. To remedy the situation, we invoke the assumptions above, which
give us the additional information in the right column: in such a treated
(control) full sample, the outcome distribution given $C$ would be the same as
that in the observed treated (control) units; and in such a cross-world full
sample, the mediator distribution given $C$ would be the same as that in the
observed control units, while the outcome distribution given $(C,M)$ would be
the same as that in the observed treated units.

Figure 2: Heuristic visualization of the identification result

This sheds light on the estimation puzzle we need to solve. If we take the
obvious approach of estimating the three potential outcome means, the task of
estimating $\text{E}[Y_{1}]$ (or $\text{E}[Y_{0}]$) would be a puzzle of
obtaining the outcome mean for a hypothetical full sample with the
distribution of $C$ from the actual full sample and the outcome distribution
given $C$ from the actual treated (control) units. The task of estimating
$\text{E}[Y_{1M_{0}}]$ would be another puzzle of obtaining the outcome mean
for a hypothetical sample with the distribution of $C$ from the actual full
sample, the mediator distribution given $C$ from the control units, and the
outcome distribution given $(C,M)$ from the treated units. This is what is
conveyed in the identification results stated above.

###  2.4 Preview of approaches and strategies for effect estimation

We will take two approaches to effect estimation: estimating the three
potential outcome means and deriving the effects from them; and estimating
effects directly. Within the first approach (Section 3), we will use the two
simple strategies of weighting and regression-based prediction, and combine
them to form a variety of estimators that range from non-robust (i.e.,
requiring all models used to be correctly specified) to more robust (i.e.,
allowing some model misspecification). These estimators are different
solutions to the puzzle posed above. Within the second approach (Section 4),
we rely on two separate ideas: emulating a randomized trial and estimating
effects as if analyzing trial data; and regression if (a proxy of) the
individual effect is available for some individuals. These two sections
effectively provide a menu of estimators with varying properties of reliance
on, and robustness to misspecification of, different models.

With respect to interval estimation (Section 5), we use a bootstrap procedure
to obtain confidence intervals that applies to all the estimators discussed.
We also derive general formulas for the asymptotic variance of the estimators.

###  2.5 Data example

We illustrate the estimators using data from The Prevention of Alcohol Use in
Students (PAS) trial in the Netherlands, which were previously used to
illustrate a method for estimating the effect on a binary outcome mediated by
multiple mediators [Nguyen et al., 2016]. In the PAS trial, middle schools
were randomized to one of four conditions: student intervention (promoting
healthy attitudes and strengthening refusal skills), parent intervention
(encouraging parental rule setting), student and parent combined intervention,
and control condition (regular biology curriculum covering effects of
alcohol). The combined intervention was effective in reducing drinking onset
[Koning et al., 2009, 2011] and drinking frequency [Koning et al., 2009], and
Koning et al. [2010] found that student attitudes towards alcohol, perceived
self-control in situations involving alcohol, and student-reported parental
rules about alcohol mediated the effect of the combined intervention and onset
of weekly drinking. In the current paper we consider the effect of the
combined intervention relative to the control condition on weekly drinking at
22 months, with the same mediators measured at six months.

It is important to note that in assembling the illustrative dataset, we make
simplifications that would not be appropriate for a full analysis: ignoring
the clustered nature of the data and leaving out observations with missing
data (mostly due to absence on data collection days, but some also due to item
nonresponse). Estimation outputs here should therefore not be interpreted as
substantive results. This use of data is simply to illustrate the estimation
strategies for a simple setting.

Table 1: Summary of baseline covariates in the analytic sample from the PAS
study

|  |  | Treated |  |  | Control |  |  | Total
---|---|---|---|---|---|---|---|---|---
|  |  | (n=533) |  |  | (n=653) |  |  | (n=1186)
Age |  |  |  |  |  |  |  |  |  |  |  |
11 |  |  | 26 | (4.9%) |  |  | 31 | (4.7%) |  |  | 57 | (4.8%)
12 |  |  | 388 | (12.8%) |  |  | 494 | (75.7%) |  |  | 882 | (74.4%)
13 |  |  | 119 | (22.3%) |  |  | 128 | (19.6%) |  |  | 247 | (20.8%)
Sex |  |  |  |  |  |  |  |  |  |  |  |
female |  |  | 221 | (41.5%) |  |  | 310 | (47.5%) |  |  | 531 | (44.8%)
male |  |  | 312 | (58.5%) |  |  | 343 | (52.5%) |  |  | 655 | (55.2%)
Religion |  |  |  |  |  |  |  |  |  |  |  |
Catholic |  |  | 40 | (7.5%) |  |  | 232 | (35.5%) |  |  | 272 | (22.9%)
Protestant/other Christian |  |  | 67 | (12.6%) |  |  | 84 | (12.9%) |  |  | 151 | (12.7%)
Islam |  |  | 27 | (5.1%) |  |  | 19 | (2.9%) |  |  | 46 | (3.9%)
Not religiously socialized |  |  | 374 | (70.2%) |  |  | 301 | (46.1%) |  |  | 675 | (56.9%)
Other |  |  | 25 | (4.7%) |  |  | 17 | (2.6%) |  |  | 42 | (3.5%)
Education tract |  |  |  |  |  |  |  |  |  |  |  |
vocational |  |  | 176 | (33.0%) |  |  | 365 | (55.9%) |  |  | 541 | (45.6%)
academic |  |  | 357 | (67.0%) |  |  | 288 | (44.1%) |  |  | 645 | (54.4%)
Baseline weekly drinking |  |  |  |  |  |  |  |  |  |  |  |
yes |  |  | 62 | (11.6%) |  |  | 105 | (16.1%) |  |  | 167 | (14.1%)
no |  |  | 440 | (82.6%) |  |  | 486 | (74.4%) |  |  | 926 | (78.1%)
no answer |  |  | 31 | (5.8%) |  |  | 62 | (9.5%) |  |  | 93 | (7.8%)
Baseline attitude |  |  |  |  |  |  |  |  |  |  |  |
negative re. alcohol use |  |  | 376 | (70.5%) |  |  | 416 | (63.7%) |  |  | 792 | (66.8%)
less negative |  |  | 157 | (29.5%) |  |  | 237 | (36.3%) |  |  | 394 | (33.2%)
Baseline parental rule |  |  |  |  |  |  |  |  |  |  |  |
strict |  |  | 402 | (75.4%) |  |  | 430 | (65.8%) |  |  | 832 | (70.2%)
not strict |  |  | 131 | (24.6%) |  |  | 223 | (34.2%) |  |  | 354 | (29.8%)
Baseline self control |  |  |  |  |  |  |  |  |  |  |  |
mean (SD) |  |  | 3.65 | (0.51) |  |  | 3.60 | (0.50) |  |  | 3.62 | (0.51)
median [min, max] |  |  | 3.69 | [1.92,4.77] |  | 3.62 | [2.00,4.92] |  | 3.62 | [1.92,4.92]

In this example $A$ is binary variable treat indicating whether a student is
in the treatment (i.e., combined intervention) or control condition. $Y$ is
binary variable drink indicating whether the student engages in weekly
drinking at 22 months. $M$ consists of three mediators measured at six months:
attitudes towards alcohol consumption (binary variable att indicating
attitudes against consumption), self-control in situations involving alcohol
(continuous variable sfc), and parental rules regarding alcohol (binary
variable rul indicating strict rules). Baseline covariates $C$ include
demographic variables age, sex, religion, education track (academic or
vocational); baseline measures of the mediators (att0, sfc0, rul0); and
baseline measure of the outcome (drink0).

The illustrative sample is of size 1186, including 533 students in the
treatment condition and 653 in the control condition. There is some covariate
imbalance in this analysis sample. The unconfoundedness assumption means that
the listed baseline covariates are sufficient to remove exposure-mediator,
exposure-outcome and mediator-outcome confounding. Table 1 summarizes the
baseline covariates in the full sample and in the treated and control
subsamples.

##  3 First approach: estimating the potential outcome means

Here the job is to estimate the means of $Y_{0},Y_{1}$ (which we refer to as
the regular potential outcomes) and the mean of $Y_{1M_{0}}$ (the cross-world
potential outcome). We address these two tasks separately in sections 3.1 and
3.2 before combining them to estimate natural (in)direct effects in section
3.3.

We label estimators with double labels: a front label for what is being
estimated (“reg” for regular and “cross” for cross-world potential outcomes)
and a back label for estimation strategy.

Proof of robustness properties of the estimators discussed is provided in the
Technical Appendix.

###  3.1 Estimating the regular potential outcome means $\text{E}[Y_{0}]$ and
$\text{E}[Y_{1}]$

We present three ways (visualized in Fig. 3) to solve the estimation puzzle
outlined earlier. The first two each use a simple strategy, weighting or
model-based prediction. The third combines both.

####  3.1.1 reg|psYobs: observed outcome in pseudo treated/control samples

With this method we directly use the observed $Y_{1}$ values in treated units
and the observed $Y_{0}$ values in control units. But instead of simply
averaging them over the treated and control subsample (which generally results
in bias), we average them over the so-called pseudo treated sample and pseudo
control sample. (In the estimator label “ps” stands for pseudo sample, and
“Yobs” for observed outcome.) These pseudo samples are obtained by weighting
treated units and control units to mimic the full sample covariate
distribution; they stand in for two of the three full samples we wish we had –
shown in Fig. 2. The weights are $\omega_{1}(C)=\mfrac{1}{\text{P}(A=1\mid
C)}$ for treated units and $\omega_{0}(C)=\mfrac{1}{\text{P}(A=0\mid C)}$ for
control units. This estimator is often known as the inverse probability
weighting (IPW) estimator.

The $\omega_{1}(C)$ and $\omega_{0}(C)$ weights are commonly estimated via
propensity score [Rosenbaum and Rubin, 1983] modeling and checking covariate
balance using inverse propensity score weights. An alternative approach is to
estimate the weights directly, finding weights that reduce the difference
between the full and pseudo samples’ covariate distributions. For example,
several methods [e.g., Imai and Ratkovic, 2014, Hainmueller, 2012] directly
target balance on covariate moments specified by the user, and one method
[Huling and Mak, 2020] minimizes a measure of distance between multivariate
distributions called energy distance [Székely and Rizzo, 2013].

Figure 3: How three reg| estimators estimate $\text{E}[Y_{1}]$

The balance desired here is three-way covariate balance, meaning balance
between each pseudo sample and the full sample and balance between the two
pseudo samples (see illustration in our data application). The propensity
score methods literature (with the average treatment effect as the central
estimand) typically focuses on treated-control covariate balance, leaving the
target covariate distribution implicit. As the current problem involves two
effects and three potential outcomes, we explicitly anchor on the full sample
covariate distribution and emphasize multi-way balance.

This estimator is nonrobust in that it is only consistent if the weights are
consistently estimated. It relies completely on the pseudo samples mimicking
well the full sample covariate distribution, but in practice there is often
some residual covariate imbalance after weighting. This estimator also often
suffers from high influence of observations with large weights.

####  3.1.2 reg|fuYpred(ss): outcome prediction on full sample based on models
fit to treated and control subsamples

Recall from Fig. 2 that if we were to observe a full treated (control) sample,
it would have the same outcome distribution given $C$ as in the actual treated
(control) units – under our assumption of unconfoundedness given $C$. This
means we can learn this distribution from the treated (control) subsample and
apply it to the full sample. We thus fit a model regressing $Y$ on $C$ in the
treated subsample, use this model to predict $Y_{1}$ for every individual in
the full sample, and average the predicted values over the full sample to
estimate $\text{E}[Y_{1}]$. For $\text{E}[Y_{0}]$, we fit the model to the
control subsample and use it to predict $Y_{0}$.111A variant of this strategy
is to combine observed $Y_{1}$ and $Y_{0}$ values (in treated and control
units, respectively) with predicted $Y_{1}$ andd $Y_{0}$ values (for control
and treated units, respectively). (Again, the label tells what the estimator
does, predicting outcomes (Ypred) on the full sample (fu) based on models fit
to subsamples (ss).)

The two models here are models for $\text{E}[Y∣C,A=1]$ and for
$\text{E}[Y∣C,A=0]$. Below we often refer to them collectively as
$\text{E}[Y∣C,a]$ (keeping implicit that $a=0,1$), but this is just an
abbreviation, as the key is that these are models that allow predicting two
different variables, $Y_{1}$ and $Y_{0}$. These models can be fit separately.
Compared to a joint model regressing $Y$ on $C,A$, the separate fitting of the
two models has the advantage that it allows tailoring the models to the
subsamples while avoiding the risk of (conscious or unconscious) fishing for a
desired treatment effect estimate.

This estimator is also nonrobust. If the outcome models are misspecified,
predictions may be poor, leading to estimation bias. The problem may be
exacerbated by extrapolation if the covariate distributions of the subsamples
(to which the models are fit) differ substantially from that of the full
sample (on which outcomes are predicted).

There is a class of estimators that are doubly robust, i.e., they combine the
outcome models and the inverse probability weights (or the propensity scores),
and are consistent if one of these two components (but not necessarily both)
is correct. Many such estimators exist [see e.g., Kang and Schafer, 2007,
Robins et al., 2007]. We consider one that builds directly on the two simple
estimators above [Robins et al., 2007, Vansteelandt and Keiding, 2011] via a
strategy that readily extends to the other estimation tasks later in the
paper.

A practical motivation for this estimator is that in most applications, we do
not know the true model that generated the data, all models are misspecified
and at best can be thought of as approximations of the true model. Other than
making the model more flexible to reduce misspecification, one way to improve
the approximation is to fit the model on (approximately) the same covariate
space on which we will use the model for prediction, which is a guard against
extrapolation. This idea characterizes the third estimator.

####  3.1.3 reg|fuYpred(ps): outcome prediction on full sample based on models
fit to pseudo treated and control samples

Like reg|fuYpred(ss), this estimator relies on model-based prediction of
$Y_{1}$ and $Y_{0}$ on the full sample. The difference is that the outcome
models used for prediction are fit to the pseudo samples (hence the (ps)
label), not to the subsamples. This is an improvement because if the weights
are (approximately) correct, the pseudo samples (approximately) share the
covariate distribution of the full sample, thus avoiding extrapolation. This
means even if the outcome models are just approximations, the approximation is
for the right covariate space. This estimator is also an improvement over
reg|psYobs in that it is protected against bias due to the pseudo samples not
perfectly mimicking the full sample, as outcome prediction is on the full
sample.222Related to the nonrobust variant in footnote 1, the corresponding
robust variant here would predict $Y_{1}$ for control units based on an
outcome model fit to a weighted treated subsample that mimics the control
subsample (using odds weights $\mfrac{\text{P}(A=0\mid C)}{\text{P}(A=1\mid
C)}$), and predict $Y_{0}$ for treated units based on an outcome model fit to
a weighted control subsample that mimics the treated subsample (using odds
weights $\mfrac{\text{P}(A=1\mid C)}{\text{P}(A=0\mid C)}$).

To make this estimator doubly robust, the outcome models are required to
satisfy a mild condition: in the sample to which the model is fit (here the
pseudo sample), the average model-predicted outcome equals the average
observed outcome [Robins et al., 2007, equation 8]. This condition, which we
label mean recovery, is satisfied by generalized linear models with canonical
link and an intercept (e.g., the usual linear regression, logistic regression,
Poisson regression). Note that this is not the only choice. For example, an
estimating equations approach may accommodate other link functions. Or if the
outcome model is fit by machine learning to improve prediction, this condition
may be achieved using targeted maximum likelihood estimation (TMLE) [van der
Laan and Rose, 2011]. These topics are outside the scope of the current paper.

Looking forward, we require this same condition for all estimators that
involve fitting outcome models to pseudo samples. It is this condition that
confers the estimators their more robust properties.

###  3.2 Estimating the cross-world potential outcome mean
$\text{E}[Y_{1M_{0}}]$

Compared to the regular potential outcome means, there are a larger number of
strategies for estimating the cross-world potential outcome mean. This is
because the latter requires combining different pieces of information from the
full sample and from the treated and control conditions – recall that the
cross-world sample in Fig. 2 has the $C$ distribution of the full sample, the
$M$ given $C$ distribution of the controls, and the $Y$ given $C,M$
distribution of the treated – and there are different ways those pieces could
be obtained and combined.

In this section we present nine estimators. The first estimator is based
purely on weighting to create a pseudo cross-world sample; it requires that
the weights are consistently estimated. The other estimators come as four
pairs. In each pair, the first estimator is nonrobust, while the second
estimator is more robust as it replaces the subsample(s) (used in the first
estimator) with corresponding pseudo sample(s), and thus allows certain
modeling elements to be misspecified if the weights that form the pseudo
sample(s) are consistent. To supplement the explanations here, Fig. 4
visualizes the gist of these estimators, and Table 2 lists their
implementation steps.

####  3.2.1 cross|pxYobs: observed outcome in pseudo cross-world sample (first
row in Fig. 4)

This estimator is in the same spirit as reg|psYobs as it is purely based on
weighting. Here the treated units are weighted to form a pseudo cross-world
sample that mimics the $C$ distribution in the full sample and mimics the $M$
given $C$ distribution in control units. (This stands in for the hypothetical
full cross-world sample that we wish we had.) We simply average the observed
outcome in this pseudo cross-world sample to estimate $\text{E}[Y_{1M_{0}}]$.
Needless to say, cross|pxYobs is only consistent if the weights are
consistently estimated.

A quick note on labels: “px” is the abbreviation for pseudo cross-world sample
(“x” symbolizing the crossing of worlds). The labels in this section use “s1”,
“s0”, “p1”, “p0”, “px” to refer specifically to the treated subsample, the
control subsample, the pseudo treated sample, the pseudo control sample, and
the pseudo cross-world sample. There is also an “sx” that we will explain
separately.

Denote the weights that form the pseudo cross-world sample out of treated
units by $\omega_{\text{x}}(C,M)$. These weights have several equivalent
expressions that point to several ways they may be estimated. We will now
pause to consider them in some detail – before continuing our survey of
different solutions to the puzzle of estimating the cross-world potential
outcome mean.

##### Three expressions for the cross-world weights.

The first expression of $\omega_{\text{x}}(C,M)$ builds on the inverse
probability weights $\omega_{1}(C)$, which weight treated units to the full
sample with respect to the covariate distribution, in a sense doing half of
the job. Such weighting does not change the mediator given covariates
distribution (which is still the distribution of $M_{1}$ given $C$). To morph
this distribution to mimic the $M_{0}$ given $C$ distribution, we use density
ratio weighting (or probability ratio if the mediator is discrete), with the
weighting function $\mfrac{\text{P}(M\mid C,A=0)}{\text{P}(M\mid C,A=1)}$,
where the numerator and denominator are the densities (or probabilities) of
the observed mediator value $M$ conditional on $C$ and on $A=0$ and $A=1$,
respectively. This weighting scheme was proposed by Hong [2010] [see also Hong
et al., 2015]. The result of the combined weighting is a pseudo sample that
mimics the $C$ distribution of the full sample and the $M$ given $C$
distribution of control units, but retains its original $Y$ given $(C,M)$
distribution (which is that of treated units). This is the desired pseudo
cross-world sample. Thus we have

$\displaystyle\omega_{\text{x}}(C,M)=\frac{1}{\text{P}(A=1\mid
C)}\frac{\text{P}(M\mid C,A=0)}{\text{P}(M\mid C,A=1)}.$ (1)

Another expression is due to the fact that by Bayes’ rule the ratio of
mediator densities above is equal to the ratio of two odds of exposure,
$\mfrac{\text{P}(A=0\mid C,M)}{\text{P}(A=1\mid
C,M)}\Big{/}\mfrac{\text{P}(A=0\mid C)}{\text{P}(A=1\mid C)}$ [noted in Zheng
and van der Laan, 2012]. The resulting expression,

$\displaystyle\omega_{\text{x}}(C,M)=\frac{\text{P}(A=0\mid
C,M)}{\text{P}(A=1\mid C,M)}\frac{1}{\text{P}(A=0\mid C)},$ (2)

is the product of two terms: an odds weight and an inverse probability weight.
(This expression appears in an identification result in Huber, 2014.)333A side
note: The odds weight in (2) also appears as part of the weight formula in
Tchetgen Tchetgen’s (2013) inverse odds ratio weighting method for estimation
of the conditional natural direct effect. Both there and here, the role of
weighting by this odds is to construct a weighted outcome distribution that
reflects the distribution of the cross-world potential outcome (conditional on
covariates). This formula provides another interpretation of the weighting: it
could be thought of as first morphing the treated subsample to mimic the joint
distribution of $(C,M)$ in the control subsample (this is what odds weighting
does), and then morphing the $C$ distribution (which now reflects the
distribution under control) to mimic that in the full sample (this is what
inverse probability weighting does).

In addition, we found a novel third expression (see derivation in the
Technical Appendix). This expression is best viewed in its version for
stabilized weights. The $\omega_{\text{x}}(C,M)$ weights in treated units have
mean equal to $\mfrac{1}{\text{P}(A=1)}$; stabilized weights are simply
$\omega_{\text{x}}(C,M)$ scaled down to mean 1 by multiplying with
$\text{P}(A=1)$. The third expression is

$\displaystyle\omega_{\text{x}}^{\text{stabilized}}(C,M)=\frac{\text{P}(C,M\mid
A=0)\frac{\text{P}(A=0)}{\text{P}(A=0\mid C)}}{\text{P}(C,M\mid A=1)},$ (3)

which could be seen as the ratio of two densities of $(C,M)$: the denominator
is the density in the treated subsample, and the weighted density in the
numerator turns out to be the density in the pseudo control sample. That is,
the weighting morphs the treated subsample such that it mimics the joint
$(C,M)$ distribution in the pseudo control sample. This makes sense, as the
pseudo control sample has the $C$ distribution of the full sample and the $M$
given $C$ distribution of control units.

For $\omega_{\text{x}}(C,M)$ estimation, the first two of the three
expressions above can be used directly as formulas. In addition to the
propensity score model, we fit either two mediator density models,
$\text{P}(M\mid C,A=0)$ and $\text{P}(M\mid C,A=1)$ (if using the first
formula), or a model for exposure given covariates and mediators,
$\text{P}(A\mid C,M)$ (if using the second formula); and plug the estimated
elements into the formula. The insight from the third expression reveals
another method: stacking the treated subsample with the pseudo control sample,
fitting a model for $A$ given $C,M$ to the stacked data, and computing
$\omega_{\text{x}}(C,M)$ as the model-predicted odds of being in the pseudo
control sample rather than the treated subsample. Of these three ways for
weights estimation based on fitting probability models, the first one has the
appeal that the models are variationally independent,444This weights
estimation method uses models for $\text{P}(A\mid C)$, $\text{P}(M\mid C,A=1)$
and $\text{P}(M\mid C,A=0)$, which correspond to a factorization of the
likelihood. These model components thus do not put constraints on one another.
but its disadvantage is that density estimation is generally a harder problem
than mean estimation (more on this later). The other two methods require
fitting fewer models (only two) and they are conditional mean models.

The third expression of $\omega_{\text{x}}(C,M)$ also provides a simple way to
use tools that estimate weights directly (moments balancing or distributional
distance minimizing): seek weights that morph the treated subsample to mimic
the pseudo control sample with respect to the joint $(C,M)$
distribution.555The appeal of direct weights estimation tools is that they
directly seek to mimic a target distribution, rather than indirectly mimic
that distribution via targeting certain probabilities. This requires data that
reflect that target distribution. This method makes use of the fact that the
pseudo-control sample reflects the target $(C,M)$ distribution.

For balance checking, the balance desired includes two components: covariate
balance between the pseudo cross-world sample and the full sample as well as
the other two pseudo samples; and covariate-and-mediator balance between the
pseudo cross-world and pseudo control samples. The checking of these balance
components is demonstrated in the data application.

Figure 4: How cross| estimators estimate $\text{E}[Y_{1M_{0}}]$

####  3.2.2 cross|p0Ypred pair: outcome prediction (given covariates and
mediators) on pseudo control sample based on model fit to treated subsample
(s1) or pseudo cross-world sample (px) (rows 2-3 of Fig. 4)

These estimators are anchored on the pseudo control sample, where the mediator
is $M_{0}$, and (if the weights are correct) the $C$ distribution mimics that
of the full sample. This give us two of the three required pieces of
information. To complete the puzzle, we replace the outcome in this pseudo
control sample (which is $Y_{0}$) with predicted $Y_{1M_{0}}$ given the units’
$(C,M_{0})$ values, where the prediction is based on a model for
$\text{E}[Y\mid C,M,A=1]$. We then average these predicted $Y_{1M_{0}}$ values
over the pseudo control sample to estimate $\text{E}[Y_{1M_{0}}]$.

With the nonrobust estimator cross|p0Ypred(s1), the $\text{E}[Y\mid C,M,A=1]$
model is estimated by regressing $Y$ on $(C,M)$ in the treated subsample. For
this estimator to be consistent, the $\omega_{0}(C)$ weights have to be
consistently estimated and this outcome model has to be correctly specified.

The more robust estimator cross|p0Ypred(px) fits the outcome model to the
pseudo cross-world sample instead of the treated subsample, and this model is
required to satisfy the mean recovery condition mentioned in section 3.1. Like
its nonrobust sibling, this estimator is not consistent if the $\omega_{0}(C)$
weights are not consistent. But if they are, then cross|p0Ypred(px) is
consistent if either the outcome model is correctly specified or the
$\omega_{\text{x}}(C,M)$ weights are consistent. That is, cross|p0Ypred(px)
has two chances to be correct, while cross|p0Ypred(s1) has only one. Note
though that if either estimator is used, it is crucial to obtain good balance
with the pseudo control sample.

Relating to the literature, cross|p0Ypred(s1) coincides with part of the
weighting-based estimator in VanderWeele and Vansteelandt [2013]. Albert
[2012] uses this strategy – using a model fit in one treatment condition to
predict the cross-world outcome on units in the other condition and weighting
those units to standardize the covariate distribution – for a general target
population.

While these two estimators rely on outcome prediction, they are weighted
estimators (as predictions are averaged over weighted units), and weighted
estimators tend to have large variance due to the variability of the estimated
weights. An alternative strategy is to find a way to predict $Y_{1M_{0}}$ on
the full sample instead of on the pseudo control sample. The next three pairs
of estimators are different versions of this.

####  3.2.3 cross|fuYpred pair: outcome prediction (given covariates) on full
sample based on model fit to pseudo cross-world subsample (sx) or pseudo
cross-world sample (px) (rows 4-5 of Fig. 4)

In the full sample, we cannot predict $Y_{1M_{0}}$ based on observed
covariates and mediators (because the observed $M$ a mixture of both $M_{0}$
and $M_{1}$). Instead, this estimator relies on $Y_{1M_{0}}$ prediction based
on covariates only. To do this, we need a model that informs of the mean of
the cross-world potential outcome given covariates, $\text{E}[Y_{1M_{0}}\mid
C]$.

The trick employed by the simpler estimator of the pair (cross|fuYpred(sx)) is
to weight the treated subsample to mimic the control subsample with respect to
the joint distribution of $(C,M)$. This weighted subsample shares features of
both subsamples; importantly, the $M$ given $C$ distribution is like in the
controls and the $Y$ given $(C,M)$ distribution is like in the treated. We
call it a pseudo cross-world subsample (abbreviated by “sx”). In this pseudo
subsample, we fit a model regressing outcome on covariates, which estimates
$\text{E}[Y_{1M_{0}}\mid C]$. (To simplify language, we loosely call this the
$\text{E}[Y_{1M_{0}}\mid C]$ model.) Based on this model, we predict
$Y_{1M_{0}}$ in the full sample and average the predicted values to estimate
$\text{E}[Y_{1M_{0}}]$.

To morph the treated subsample to mimic the control subsample’s joint $(C,M)$
distribution, we use odds weights
$\omega_{\text{sx}}(C,M)=\mfrac{\text{P}(A=0\mid C,M)}{\text{P}(A=1\mid
C,M)}$, which can be estimated based on a $\text{P}(A\mid C,M)$ model.666An
alternative is to mimic the $M$ given $C$ distribution using density ratio
weights $\mfrac{\text{P}(M\mid C,A=0)}{\text{P}(M\mid C,A=1)}$; we do not
recommend this because there is no simple way to check balance on conditional
distributions.

cross|fuYpred(sx) is nonrobust. For it to be consistent, the
$\omega_{\text{sx}}(C,M)$ weights have to be consistently estimated and the
$\text{E}[Y_{1M_{0}}\mid C]$ model has to be correctly specified.

Table 2: Implementation steps of reg| and cross| estimators. nonR, MR and R =
nonrobust, more robust and robust estimators. These are combined to form
effect pair estimators in Table LABEL:tab:regcrossEstimators.

reg|psYobs 1. Estimate $\omega_{1}(C)$ and $\omega_{0}(C)$ weights* 2. Average
the observed outcome in the pseudo control and pseudo treated samples
---
reg|fuYpred pair: nonrobust (ss) and robust (ps) 1. For the R estimator,
estimate $\omega_{1}(C)$ and $\omega_{0}(C)$ weights* 2. Model $Y$ given $C$
in the treated and control subsamples (nonR) or pseudo treated and control
samples (R) 3. Based on these models, predict $Y_{1}$ and $Y_{0}$ given $C$ in
full sample 4. Average predicted $Y_{1}$ and $Y_{0}$ in full sample
cross|pxYobs 1. Estimate $\omega_{\text{x}}(C,M)$ weights* 2. Average the
observed outcome in the pseudo cross-world sample
cross|p0Ypred pair: nonrobust (s1) and more robust (px) 1. Estimate
$\omega_{0}(C)$ weights* 2. For the MR estimator, also estimate
$\omega_{\text{x}}(C,M)$ weights* 3. Model $Y$ given $C,M$ in the treated
subsample (nonR) or pseudo cross-world sample** (MR) 4. Based on model,
predict $Y_{1M_{0}}$ given $C,M$ in control units 5. Average predicted
$Y_{1M_{0}}$ in pseudo control sample
cross|fuYpred pair: nonrobust (sx) and more robust (px) 1. Estimate
$\omega_{\text{x}}^{\text{sub}}(C,M)$ weights (nonR) or
$\omega_{\text{x}}(C,M)$ weights* (MR) 2. Model $Y$ given $C$ in the pseudo
cross-world subsample (nonR) or pseudo cross-world sample** (MR) 3. Based on
model, predict $Y_{1M_{0}}$ given $C$ in full sample 4. Average predicted
$Y_{1M_{0}}$ in full sample
cross|fuY2pred pair: nonrobust (s1s0) and robust (pxp0) 1. For the R
estimator, estimate $\omega_{0}(C)$ and $\omega_{\text{x}}(C,M)$ weights* 2.
Model $Y$ given $C,M$ in the treated subsample (nonR) or pseudo cross-world
sample** (R) 3. Based on model, predict $Y_{1M_{0}}$ given $C,M$ in control
units 4. Model predicted $Y_{1M_{0}}$ given $C$ in the control subsample
(nonR) or pseudo control sample** (R) 5. Based on model, predict $Y_{1M_{0}}$
given $C$ in full sample 6. Average predicted $Y_{1M_{0}}$ in full sample
cross|fuMsimYpred pair: nonrobust (s0s1) and more robust (p0px) 1. For the MR
estimator, estimate $\omega_{0}(C)$ and $\omega_{\text{x}}(C,M)$ weights* 2.
Model the density of $M$ given $C$ in the control subsample (nonR) or pseudo
control sample (MR) 3. Model $Y$ given $C,M$ in the treated subsample (nonR)
or pseudo cross-world sample** (MR) 4. Do many times in full sample: i. Based
on first model, simulate $M_{0}$ given $C$ ii. Based on second model, predict
$Y_{1M_{0}}$ given the combination of $C$ and predicted $M_{0}$ 5. Average all
predicted $Y_{1M_{0}}$ values in full sample
* Weights estimation includes balance checking.
* This model for the MR/R estimator is required to satisfy the mean recovery condition.

With the more robust cross|fuYpred(px) estimator, the outcome given covariates
model is fit to the pseudo cross-world sample (px) instead. We require this
model to satisfy the mean recovery condition. This estimator has two chances
to be correct: (1) if the weights $\omega_{\text{x}}(C,M)$ are consistently
estimated, cross|fuYpred(px) is consistent even if the
$\text{E}[Y_{1M_{0}}\mid C]$ model is misspecified; and (2) if the outcome
model is correctly specified and only the mediator-related part of
$\omega_{\text{x}}(C,M)$ is correct, cross|fuYpred(px) is consistent.

To clarify the second case, the mediator-related part of the weights is the
term that controls the mediator distribution in the pseudo cross-world sample.
It varies by the weights estimation method: in the first and second methods it
is $\mfrac{\text{P}(M\mid C,A=0)}{\text{P}(M\mid C,A=1)}$ and
$\mfrac{\text{P}(A=0\mid C,M)}{\text{P}(A=1\mid C,M)}$, respectively; in the
third method it is the odds of being in the pseudo control sample rather than
the treated subsample (where the pseudo control sample may be incorrectly
weighted). When only the mediator-related part of the weights is correct, the
weighting gets the $C$ distribution wrong but gets the $M$ given $C$
distribution right, and it is the latter that ensures that the treated units’
$Y$ given $C$ distribution is appropriately morphed into the target
$Y_{1M_{0}}$ given $C$ distribution. As the outcome regression model
conditions on $C$, it is (if correctly specified) not affected by the
incorrectly weighted $C$ distribution.

It should be pointed out that while these two estimators are based on
$Y_{1M_{0}}$ prediction for the full sample, they are essentially also
weighted estimators because they rely completely on weighting to obtain data
that reflect the distribution of $Y_{1M_{0}}$ given $C$. It is crucial to
achieve good balance with the pseudo cross-world sample/subsample.

####  3.2.4 cross|fuY2pred pair: outcome prediction on full sample based on
double outcome model fit – to treated and control subsamples (s1s0) or to
pseudo cross-world and pseudo control samples (pxp0) (rows 6-7 of Fig. 4)

These estimators also rely on predicting $Y_{1M_{0}}$ in the full sample based
on an estimated $\text{E}[Y_{1M_{0}}\mid C]$ model. Here this model is
estimated in two steps: first fitting a $\text{E}[Y\mid C,M,A=1]$ model and
using it to predict $Y_{1M_{0}}$ in control units, then regressing the
predicted $Y_{1M_{0}}$ on covariates to estimate the $\text{E}[Y_{1M_{0}}\mid
C]$ model. (We loosely refer to these two models with one building directly on
the other as a double outcome model fit.)

With the nonrobust estimator cross|fuY2pred(s1s0), these two models are fit to
the treated and control subsamples, respectively. For this estimator to be
consistent, both models have to be correctly specified.

The robust estimator cross|fuY2pred(pxp0) fits the outcome models to the
pseudo cross-world and pseudo control samples instead; and these models are
required to satisfy the mean recovery condition. This estimator has three
chances to be consistent: (1) both of the outcome models are correctly
specified; (2) both the $\omega_{\text{x}}(C,M)$ and $\omega_{0}(C)$ weights
(that form pseudo cross-world and pseudo control samples) are consistently
estimated; or (3) the $\omega_{\text{x}}(C,M)$ weights are consistently
estimated and the $\text{E}[Y\mid C,A=1,M]$ model is correctly specified.

A side technical note: These robustness properties of cross|fuY2pred(pxp0) are
similar to those of the robust estimator in Tchetgen Tchetgen and Shpitser
[2012] that solves the efficient influence function (EIF) based estimating
equation (equation 4 in that paper). cross|fuY2pred(pxp0) is also based on the
EIF, but its steps are intuitively meaningful without requiring knowledge of
the EIF.

Notice that we have called cross|fuY2pred(pxp0) robust rather than more
robust. That reflects the fact that this estimator does not require any single
estimation component to be correct. The estimators dubbed more robust, on the
other hand, are more robust than their nonrobust siblings, but all depend on a
certain estimation component being correct.

####  3.2.5 cross|fuMsimYpred pair: mediator simulation and outcome prediction
on full sample based on models fit to control and treated subsamples (s0s1) or
pseudo control and cross-world samples (p0px) (rows 8-9 of Fig. 4)

These estimators involve fitting a model for the conditional mediator density
$\text{P}(M\mid C,A=0)$ and a model for the conditional outcome mean
$\text{E}[Y\mid C,M,A=1]$. Having learned these models, we put the observed
mediator and outcome aside. For all units in the full sample, we simulate
$M_{0}$ based on the first model, and with the simulated $M_{0}$ and observed
$C$, predict $Y_{1M_{0}}$ based on the second model. We do this multiple
times, resulting in multiple sets of predicted $Y_{1M_{0}}$ values. We average
these predicted values to estimate $\text{E}[Y_{1M_{0}}]$.777A variant is to
use the observed mediator in control units and only simulate $M_{0}$ for
treated units.

With the nonrobust estimator cross|fuMsimYpred(s0s1), the mediator model is
fit to the control subsample and the outcome model is fit to the treated
subsample. For this estimator to be consistent, both models have to be
correctly specified.

This mediator simulation strategy was used in Imai et al.’s (2010a) method for
natural (in)direct effects estimation, implemented in the R package mediation
[Tingley et al., 2014]. Their method differs from the current estimator in
that it uses this strategy to estimate not only the cross-world potential
outcome means but all potential outcome means (thus relying on more models),
and its implementation uses models for $\text{P}(M\mid C,A)$ and
$\text{E}[Y\mid C,A,M]$ fit to the full sample rather than exposure-condition-
specific models.

The more robust cross|fuMsimYpred(p0px) instead fits the mediator model to the
pseudo control sample and the outcome model to the pseudo cross-world sample,
with the outcome model safisfying the mean recovery condition. Like its
nonrobust sibling, this estimator is inconsistent if the mediator density is
misspecified. Assuming correct specification of this model, this estimator has
two chances to be correct: either the $\omega_{\text{x}}(C,M)$ weights are
consistently estimated or the outcome model is correctly specified.888In the
more robust version of the variant mentioned in foot note 7, the mediator
model is fit to a weighted control subsample that mimics the treated
subsample’s $C$ distribution – using odds weights $\mfrac{\text{P}(A=1\mid
C)}{\text{P}(A=0\mid C)}$.

Relating to the existing literature, a version of cross|fuMsimYpred(p0px)
where the $\omega_{\text{x}}(C,M)$ weights are estimated based on mediator
density models is an implementation of the estimator in section 5 in Tchetgen
Tchetgen and Shpitser [2012]. This estimator integrates the conditional
outcome mean function $\text{E}[Y\mid C,M,A=1]$ over the conditional mediator
density $\text{P}(M\mid C,A=0)$; the simulation-prediction-averaging procedure
here is a numerical evaluation of that integral.

Note that these estimators’ dependence on correct specification of the
mediator density model is an important drawback, as density estimation is a
harder problem than mean estimation, an issue also raised by Albert [2012].
(Estimators that do not involve simulating $M_{0}$, or more generally,
integrating over an estimated conditional density of $M_{0}$, avoid this
problem.) For example, with a continuous variable, if the conditional mean is
of interest, a common model choice is the linear model, which assumes a
specific functional form for the mean but makes no other assumption. If the
conditional density is of interest, one might still use the linear model but
has to make additional distributional assumptions (e.g., the error is normally
distributed or follows some other distribution) which are likely incorrect. In
the special case where the mediator is a single binary variable, the
distribution is fully described by the probability so the model reduces to a
conditional mean model.

In a setting with a multivariate mediator (like in our data example), we need
to model the joint distribution of the mediators given covariates in control
units. To do this, we factor the joint into conditional
densities/probabilities. Let $M=(M^{a},M^{b},M^{c})$ where $M^{a},M^{b},M^{c}$
are three mediators.

$\text{P}(M\mid C,A=0)=\text{P}(M^{a}\mid C,A=0)\,\text{P}(M^{b}\mid
C,A=0,M^{b})\,\text{P}(M^{c}\mid C,A=0,M^{a},M^{b}).$

In the control subsample (or the pseudo control sample), we fit three models
for the three mediators. All three models condition on $C$, the second model
conditions additionally on $M^{a}$, and the third model conditions
additionally on both $M^{a}$ and $M^{b}$. The order of variables in the
factorization can be chosen for modeling convenience (see our data application
for an example). Simulation follows the order of the fitted models.

To sum up, in this section we have presented a purely weighting-based
estimator (cross|pxYobs) followed by four pairs of estimators: the
cross|p0Ypred pair relies on outcome prediction on the pseudo control sample,
the cross|fuYpred pair relies on outcome prediction on the full sample based
on a model fit to the pseudo cross-world (sub)sample, the cross|fuY2pred pair
relies on outcome prediction on the full sample based on double outcome model
fit, and the cross|fuMsimYpred pair relies on mediator simulation before
outcome prediction. In each pair, the second estimator is more robust than the
first as it does not require that all estimation components are correct. Among
the more robust estimators, cross|fuY2pred(p0px) is the most robust as it does
not require any specific estimation component to be correctly
specified/consistent; on the other hand cross|p0Ypred(px) and
cross|fuYpred(px) are not robust to inconsistent weights, and
cross|fuMsimYpred(p0px) is not robust to misspecification of the conditional
mediator density model.

###  3.3 Combining the three potential outcome means to estimate the effects

Effects are estimated by the differences of outcome mean estimates from reg|
and cross| estimators, using the definitions
$\mathrm{NDE}_{0}=\text{E}[Y_{1M_{0}}]-\text{E}[Y_{0}]$ and
$\mathrm{NIE}_{1}=\text{E}[Y_{1}]-\text{E}[Y_{1M_{0}}]$. Table
LABEL:tab:regcrossEstimators combines reg| and cross| estimators based on
similar ideas into estimators of this effect pair. The table lists the reg|
and cross| being combined, but from here on we will also label each
combination by combining the labels of the elements, e.g., POs|psYobs-pxYobs
is the combination of reg|psYobs and cross|pxYobs.

The table includes the purely weighting-based estimator (which is nonrobust,
marked with ~~R~~) and five pairs of estimators, each pair containing a
nonrobust estimator (also marked with ~~R~~) and a more robust sibling (marked
with MR for more robust or R for robust). For each of these MR/R estimators,
the table lays out the different subsets of estimation components that if
consistent ensure the estimator’s consistency. The four MR estimators, while
more robust than their nonrobust siblings, depend on specific modeling
components being correct: one is not robust to inconsistent $\omega_{1}(C)$
and $\omega_{0}(C)$ weights, one is not robust to inconsistent $\omega_{0}(C)$
weights, one is not robust to inconsistent estimation of the mediator-related
part of $\omega_{\text{x}}(C,M)$ weights, and one is not robust to
misspecification of the mediator density model $\text{P}(M\mid C,A=0)$. The
one estimator marked R, POs|fuYpred(ps)-fuY2pred(pxp0), does not depend on a
specific modeling component being correct.

Tying our results to the literature, several of the ~~R~~ estimators are
connected to known estimators. The purely weighting-based POs|psYobs-pxYobs
has a version that is Hong’s (2010) original mediator-probability-ratio
weighting estimator, which also appears in Lange et al. [2012].
POs|psYobs-p0Ypred(s1) is the weighting-based estimator among the estimators
in VanderWeele and Vansteelandt [2013]; the cross| part of this estimator uses
the same strategy used in Albert [2012]. POs|fuYpred(ss)-fuMsimYpred(s0s1) is
similar to Imai et al. [2010b]’s quasi-Bayesian estimator (which is also
nonrobust); the difference is that Imai et al.’s estimator uses the mediator
simulation strategy to estimate all potential outcome means, not just the
cross-world ones. We have not seen the other ~~R~~ estimators in our reading
of the literature.

Among the MR estimators, the cross| component of
POs|fuYpred(ps)-fuMsimYpred(p0px) contains a version that coincides with the
estimator in section 5 of Tchetgen Tchetgen and Shpitser [2012]. We have not
seen the other MR estimators in our reading of the literature. Interestingly,
it seems the robust estimator POs|fuYpred(ps)-fuY2pred(pxp0) also has not
appeared in the literature. It is however a close sibling to an estimator of
the natural direct effect which we will cover in section 4.1.1 (the latter may
be thought of as an extension of the former), and that other estimator is
closely related to an existing estimator.

Like many estimation methods for mediation analysis, quite a few of the
estimators here rely on cross-world weighting using $\omega_{\text{x}}(C,M)$
weights. We note that the cross-world weights have several equivalent
expressions including two previously known [Hong, 2010, Huber, 2014] and a
novel one, expanding options for the estimation of these weights. This is
relevant to all estimators that involve weighting to construct the pseudo
cross-world sample, i.e., all the more robust estimators plus the pure
weighting estimator.

###  3.4 Data example application

####  3.4.1 Weights estimation, weights checking and balance checking

For the purpose of illustrating the different effect estimators, we choose to
estimate weights via parametric probability models, for this approach is
familiar to many applied researchers. $\omega_{1}(C)$ and $\omega_{0}(C)$
weights are estimated based on a propensity score model.
$\omega_{\text{x}}(C,M)$ weights are estimated using the second method based
on the propensity score model and a model for exposure given covariates and
mediators. (We opt not to use the first method which would require fitting six
models for the three mediators.) We estimate both of the exposure models by
logistic regression, with the following model formulas, which include some
interactions and spline terms.999The ns() function implements cubic splines
with 4 degrees of freedom on the continuous variable sfc0. These are the
result of several rounds of model fitting and balance checking.

$\text{P}(A\mid C)$: | treat $\sim$ sex + age + edu + religion + drink0 + att0 * rul0 * ns(sfc0,4)
---|---
$\text{P}(A\mid C,M)$: | treat $\sim$ sex + age + edu + religion + drink0 + att0 * rul0 * att * rul + ns(sf0,4) * ns(sfc,4)

Figure 5: Data example application: Distributions of weights for the pseudo
control, treated and cross-world samples. For comparability, stabilized
weights are shown. ess = effective sample size.
Figure 6: Data example application: Balance checking

The distributions of the $\omega_{0}(C)$, $\omega_{1}(C)$ and
$\omega_{\text{x}}(C,M)$ weights are shown in Fig 5. None of the weights are
extreme. Balance checks are presented in Fig 6. The top part of this figure
shows covariate balance of the pseudo treated and control samples (relative to
the full sample and to each other). Covariate balance (in terms of categorical
variable prevalence and continuous variable density) noticeably improves as a
result of weighting: compared to the treated and control subsamples, the
pseudo treated and control samples mimic the full sample much better with
regards to covariate distribution. The bottom part of Fig 6 show covariate-
and-mediator balance of the pseudo cross-world sample (relative to the full
sample and the other pseudo samples). Covariate balance is slightly worse for
the pseudo cross-world sample than for the other two pseudo samples. Mediator
balance is substantially improved after weighting. In this example, the
weighting achieves quite good balance. In cases where substantial imbalance
remains, that might be a concern for estimators that rely heavily on
weighting.

For the cross|fuYpred(sx) estimator, we estimate the $\omega_{\text{sx}}(C,M)$
weights using the same $\text{P}(A\mid C,M)$ model above. Weights distribution
and balance are not shown.

Figure 7: Data example application: Effect estimates from all estimators
(point estimates and 95% percentile bootstrap intervals) – for super
population inference. We leave out labels for where models are fit, and
instead mark the more robust and robust estimators with .MR and .R.

####  3.4.2 Outcome and mediator models

The various estimators call for outcome-given-covariates-and-mediators and
outcome-given-covariates models fit to subsamples and pseudo samples. As the
outcome is binary, we use logistic regression with the following model
formulas

$\text{E}[Y\mid C,A=a]$: |  drink $\sim$ sex + age + edu + religion + drink0 + att0 + rul0 + ns(sfc0,3) (for $a=1,0$)
---|---
$\text{E}[Y_{1M_{0}}\mid C]$: | drink $\sim$ sex + age + edu + religion + drink0 + att0 + rul0 + ns(sfc0,3)
$\text{E}[Y\mid C,M,A=1]$: | drink $\sim$ sex + age + edu + religion + drink0 + att0 + rul0 + ns(sfc0,3) + att + rul + ns(sf0,3)

The cross|fuMsimYpred estimators require mediator density modeling in the
control subsample or pseudo control sample. We fit logit models for the two
binary mediators att and rul, and a linear model for the continuous mediator
sfc, with the following formulas

$\text{P}(M^{a}\mid C,A=0)$: | att $\sim$ age + sex + edu + religion + att0 + rul0 + ns(sfc0,3) + drink0
---|---
$\text{P}(M^{b}\mid C,M^{a},A=0)$: | rul $\sim$ age + sex + edu + religion + att0 + rul0 + ns(sfc0,3) + drink0 + att
$\text{P}(M^{c}\mid C,M^{a},M^{b},A=0)$: | sfc $\sim$ age + sex + edu + religion + att0 + rul0 + ns(sfc0,3) + drink0 + att + rul

In simulating $M^{c}_{0}$ (sfc), we assume errors are normally distributed.
(Of the three mediators, we choose to model sfc last for convenience. This
avoids having to specify models with sfc as another continuous predictor,
which would be more complicated.)

####  3.4.3 Results

Effect estimates from the eleven POs| estimators (as well as from the direct
estimators of effects to be discussed in Section 4) are presented in Fig. 7.
The confidence intervals are from bootstrapping (see details in Section 5).
The estimates are quite similar. As expected, weighting estimators have larger
variance; and within each pair of nonrobust/(more) robust estimators, the
(more) robust one tends to have larger variance than the nonrobust one.
Overall, it appears that the effect of the intervention on drinking consists
of a small part mediated by the mediators being considered (alcohol-related
attitudes, parental rules, and self-control), and a larger direct effect.

##  4 Second approach: effects estimation

For direct estimation of effects, we consider two different approaches. The
first approach is based on the view of the marginal additive effect as the
mean of the individual additive effect; thus it focuses on obtaining (for a
set of units) a reasonable proxy for the unobserved variable that is the
individual effect. The second approach treats pseudo samples as if they were
arms of a randomized trial and taps into methods developed to estimate the
average treatment effect in randomized trials.

###  4.1 Marginal additive effect as the mean of individual effects

The idea is to find a proxy for the individual effect that has the same mean
given covariates as the effect itself. For control units, one proxy for the
individual $N\\!D\\!E_{0}=Y_{1M_{0}}-Y_{0}$ is the difference between
predicted $Y_{1M_{0}}$ (based on an appropriate model) and observed $Y_{0}$.
This simple idea leads to the pair of estimators below. Note that this
approach is asymmetric, as no similar proxy for the individual
$N\\!I\\!E_{1}=Y_{1}-Y_{1M_{0}}$ is available.

####  4.1.1 NDE|fuNDEpred pair: NDE prediction on full sample based on a
regression of the difference between predicted $Y_{1M_{0}}$ and observed
$Y_{0}$, where the model used to predict $Y_{1M_{0}}$ and the difference
regression are fit to subsamples (s1s0) or pseudo samples (pxp0)

These estimators are a direct extension of the cross|fuY2pred pair, which
relies on predicting $Y_{1M_{0}}$ on control units and regressing it on
covariates. The current estimators also predict $Y_{1M_{0}}$ on control units
but instead regress the difference between predicted $Y_{1M_{0}}$ and observed
$Y_{0}$ on covariates.

These estimators involve fitting a model for $\text{E}[Y\mid C,M,A=1]$ in the
treated subsample (nonrobust version) or the pseudo cross-world sample (robust
version), and using it to predict $Y_{1M_{0}}$ in control units. In these
units, the difference between this predicted $Y_{1M_{0}}$ and the observed
outcome serves as the proxy for the unit’s $N\\!D\\!E_{0}$. The next step is
to regress this difference on covariates in the control subsample (nonrobust)
or pseudo control sample (robust); we call this the
$\text{E}[N\\!D\\!E_{0}\mid C]$ model. This model is then used to predict the
individual $N\\!D\\!E_{0}$ for all units in the full sample, and these
predicted individual effects are averaged to estimate the average
$\text{NDE}_{0}$. The conditional mean and effect models used for the robust
estimator are required to satisfy the mean recovery condition.

For the nonrobust NDE|fuNDEpred(s1s0) to be consistent, both the outcome model
$\text{E}[Y\mid C,M,A=1]$ and effect model $\text{E}[N\\!D\\!E_{0}\mid C]$
need to be correctly specified. The robust NDE|fuNDEpred(pxp0), on the other
hand, has three chances to be correct: (i) if both the outcome and effect
models are correctly specified; or (ii) if the $\omega_{0}(C)$ and
$\omega_{\text{x}}(C,M)$ weights that form the two pseudo samples are
consistent; or (iii) if the $\omega_{\text{x}}(C,M)$ weights and the
$\text{E}[Y\mid C,A=1,M]$ model are consistent.

Table 3: Implementation steps of direct estimators of effects. These are
combined to form several estimators of the effect pair in Table
LABEL:tab:effectEstimators. nonR, MR and R = nonrobust, more robust and
robust.

NDE|fuNDEpred pair: nonrobust (s1s0) and robust (pxp0) 1. For the R estimator,
estimate $\omega_{0}(C)$ and $\omega_{\text{x}}(C,M)$ weights* 2. Model $Y$
given $C,M$ in the treated subsample (nonR) or pseudo cross-world sample** (R)
3. Based on model, predict $Y_{1M_{0}}$ given $C,M$ in control units 4. In
control units compute proxy of the individual $\text{NDE}_{0}$ as predicted
$Y_{1M_{0}}$ minus $Y$ 5. Model the $\text{NDE}_{0}$ proxy given $C$ in the
control subsample (nonR) or pseudo control sample** (R) 6. Based on model,
predict $\text{NDE}_{0}$ given $C$ in full sample 7. Average predicted
$\text{NDE}_{0}$ in full sample
---
NDE&NIE|pxsCadj 1. Estimate $\omega_{1}(C)$, $\omega_{0}(C)$ and
$\omega_{\text{x}}(C,M)$ weights* (to mimic 3-arm trial) 2. Covariate adjust
using a working model (see the bottom row of the table)
NIE|psYpredCadj pair: nonrobust (s1) and more robust (px) 1. Estimate
$\omega_{1}(C)$ and $\omega_{0}(C)$ weights* (to mimic 2-arm trial) 2. Model
$Y$ given $C,M$ in the treated subsample (nonR) or pseudo cross-world sample**
(MR) 3. Based on model, predict $Y_{1M_{0}}$ in pseudo control sample and
replace observed outcome 4. Covariate adjust using a working model (see the
bottom row of the table)
TE|psCadj 1. Estimate $\omega_{1}(C)$ and $\omega_{0}(C)$ weights* (to mimic
2-arm trial) 2. Covariate adjust using a working model (see the bottom row of
the table)
Covariate adjustment steps: • stack the trial-mimicking weighted data • fit a
working model*** for the outcome given arm and covariates • if linear model,
take the coefficient(s) of the variable(s) indicating arm as the effect
estimate(s) • if nonlinear model, use fitted model to predict potential
outcomes on the full sample as and take the differences of potential outcome
means.
* Weights estimation includes balance checking.
* This model for the MR/R estimator is required to satisfy the mean recovery condition.
** The working model (not assumed to be correct) is required to satisfy the
group mean recovery condition.

An aside: Instead of using the observed outcome in the construction of the
proxy for $N\\!D\\!E_{0}$ in control units, a variant of these estimators
replaces it with a predicted value of this outcome based on a $\text{E}[Y\mid
C,M,A=0]$ model. The robust version of this variant (where all models are fit
to relevant pseudo samples) is closely related to the TMLE estimator in Zheng
and van der Laan [2012]. (Both solve a common set of estimating equations.
While we restrict our attention to parametric models here, the TMLE estimator
can fit models using nonparametric (machine learning) algorithms and satisfy
the mean recovery condition.) This alternative estimator also has three
chances to be correct under similar conditions to those listed above, except
that condition (i) additionally requires correct specification of the
$\text{E}[Y\mid C,M,A=0]$ model.

###  4.2 Treating pseudo samples as mimicking a randomized trial

The estimators presented so far could be seen as arising from (and extending)
the tradition of estimating the mean of a variable in the presence of missing
data. We now consider a different approach, prompted by a practice commonly
seen in applied research of first (i) preprocessing data to balance covariates
and then (ii) fitting a model regressing outcome on treatment and covariates.
The approach we take assumes that the preprocessing that forms the pseudo
samples is correct and treats the pseudo samples as mimicking a randomized
trial. We use regression to reduce estimation variance, but do not assume the
regression model is correct [in the spirit of Wang et al., 2019, Steingrimsson
et al., 2017].

Note that there are other approaches and goals that result in estimators with
the same two steps. Keeping the same marginal effect estimands, a doubly
robust approach (similar to what we did earlier when estimating outcome means)
would result in estimators that are consistent for the effects if either the
preprocessing components are correct or the outcome regression is correct for
all the potential outcomes being contrasted. These estimators turn out to be
equivalent to several of the POs| estimators from section 3, so we do not
discuss this approach separately. In a different scenario where the goal is to
estimate effects conditional on covariates, there are estimators with the same
two steps [see e.g., Steen et al., 2017] that are consistent if the
preprocessing components are correct AND the outcome regression is correct for
the conditional effects. The approach we are considering (for marginal effects
estimation) refrains from counting on the outcome regression model to be
correct for either potential outcomes or effects (this is mindful of the fact
that applied researchers may use quite simple outcome regression models).
Instead, it counts on successful construction of the pseudo samples, which
require that the weights are correct. This admittedly is a tall order, but we
think there is more hope for this, say, through the common practice of balance
checking.

This approach treats the pseudo samples as mimicking a randomized trial, and
borrows methods for analyzing trial data. Let us for a moment put aside the
mediation analysis and consider the task of estimating the average additive
treatment effect (ATE) with a two-arm randomized trial. The simplest estimator
of the ATE is the unadjusted (difference in means) estimator. It is unbiased
but has large variance. Adjusting for covariates helps improve estimation
precision. For a continuous outcome, a simple covariate adjusted estimator is
the ANCOVA estimator – the coefficient of treatment in a linear model [Wang et
al., 2019] regressing the outcome on treatment and covariates (with no
treatment-covariate interaction terms). More generally, covariate adjustment
may be done by fitting a generalized linear model regressing the outcome on
exposure and covariates and using this model to predict both potential
outcomes for every unit (in both arms), and averaging their difference over
all the units [Steingrimsson et al., 2017]. In each case, the regression model
is a working model, meaning it is not assumed to be correct, but is used
simply as a tool to leverage covariates to soak up variance. The estimator’s
consistency is due to the treatment being randomized, and due to the working
model satisfying the group mean recovery condition, i.e., mean recovery within
each treatment arm. Other covariate adjustment methods for randomized trials
[Colantuoni and Rosenblum, 2015] may also be used, but we will restrict our
attention to these simple methods. Borrowing these methods lead to several
estimators below for our current problem.

####  4.2.1 NDE&NIE|psxCadj: estimating both NDE and NIE with pseudo samples
mimicking a three-way randomized trial using covariate adjustment

Assuming the $\omega_{1}(C),\omega_{0}(C)$ and $\omega_{\text{x}}(C,M)$
weights are correct, we treat the pseudo treated, control and cross-world
samples as mimicking a three-arm randomized trial, with the three outcome
distributions we wish to contrast. The simplest estimator of the pair
$(\text{NDE}_{0},\text{NIE}_{1})$ is the difference-in-means estimator
(labeled NDE&NIE|psxYobs), which exactly coincides with the POs|psYobs-pxYobs
estimator from section 3.3. Applying the working model covariate adjustment
method above to this pseudo trial obtains estimators we label NDE&NIE|psxCadj.

Consider first covariate adjustment by linear regression. With three arms
there are two versions. One version is to estimate the two effects separately:
for each effect, we stack the two relevant pseudo samples, add a variable
indicating arm (1 = cross-world vs. 0 = control for $\text{NDE}_{0}$
estimation, and 1 = treated vs. 0 = cross-world for $\text{NIE}_{1}$
estimation), fit a main effects linear model regressing the outcome on arm and
covariates, and take the coefficient of the arm variable as the effect
estimate. The other version is to stack all three pseudo samples and use a
single working model regressing outcome on arm and covariates. Both versions
are consistent (assuming the weights are correct). If the covariate parts of
the two models in the first version are similar, the second version tends to
have smaller variance; if they are not similar, the second version tends to
have larger variance.

If a working generalized linear (instead of a linear) model is used, there are
also the two options for joint or separate estimation of the two effects. In
addition, the procedure is slightly modified in that after fitting the model,
the outcome prediction and averaging is done on the original full sample
instead of on the pseudo samples.

For clarity, two comments are warranted. First, while this last version of the
NDE&NIE|psxCadj estimator with outcome prediction looks similar to the
POs|fuYpred(ps)-fuYpred(px) estimator in that both involve model-based
prediction, this resemblance is superficial. The latter has the robustness
property that even with inconsistent $\omega_{1}(C)$ and $\omega_{0}(C)$
weights, the estimator is still consistent if the conditional outcome mean
models are correctly specified and the mediator-related component of
$\omega_{\text{x}}(C,M)$ is consistent. NDE&NIE|psxCadj, on the other hand,
relies on consistent estimation of all $\omega_{1}(C),\omega_{0}(C)$ and
$\omega_{\text{x}}(C,M)$ weights, as we do not count on the working regression
model to be correct for the conditional outcome means.

Second, NDE&NIE|psxCadj bears a resemblance to certain estimators of
conditional natural (in)direct effects – e.g., those described in Lange et al.
(2012, eq. 5) and implemented in Steen et al.’s (2017) R package medflex.
These other estimators use the same weighting scheme that forms the pseudo
samples here, and rely on a generalized linear model to estimate conditional
natural (in)direct effects given covariates (where the link function
(identity/log/logit) matches the scale on which the conditional effects are
defined). The key difference, though, is that these estimators assume that
such a model is correctly specified for the conditional effects. For the
NDE&NIE|psxCadj estimator of the marginal additive effects, the model we use
is simply a working model.

####  4.2.2 NIE|psYpredCadj pair: estimating NIE by covariate adjustment, with
pseudo treated and control samples (the latter with predicted $Y_{1M_{0}}$)
mimicking a two-way randomized trial

These estimators of the NIE rely on the pseudo treated and pseudo control
samples, but replace the outcome in control units with the predicted
$Y_{1M_{0}}$ based on a $\text{E}[Y\mid C,A=1,M]$ model.101010Variants of
these estimators may replace the observed outcome in the pseudo treated sample
with predictions of that outcome given covariates (or covariates and
mediators). This is unnecessary given the covariate adjustment that will be
accomplished via the working model.

Again, we take the difference-in-means estimator as the starting point,
labeling it NIE|psYpred(s1) or NIE|psYpred(px) depending on whether the
$\text{E}[Y\mid C,A=1,M]$ model is fit to the treated subsample or to the
pseudo cross-world sample. (These are exactly the NIE part of the estimators
POs|psYobs-p0Ypred(s1) and POs|psYobs-p0Ypred(px) from section 3.3.) Applying
the covariate adjustment method above obtains the estimators
NIE|psYpred(s1)Cadj and NIE|psYpred(px)Cadj with smaller variance.

These covariate adjusted estimators assume that the $\omega_{1}(C)$ and
$\omega_{0}(C)$ weights are consistently estimated. In addition,
NIE|psYpred(s1)Cadj also assumes that model for $\text{E}[Y\mid C,M,A=1]$ is
correctly specified; its sibling NIE|psYpred(px)Cadj is more robust as it
requires either this model or the $\omega_{\text{x}}(C,M)$ weights to be
correct (but not both).

Similar to NDE&NIE|psxCadj, the current estimators have some apparent
resemblance to the NIE part of the POs|fuYpred(ps)-p0Ypred(px), and to the
imputation-based estimators of conditional effects described in Vansteelandt
et al. [2012] and implemented in the medflex package. Again, a key difference
is that the regression model we use is a working model that we do not assume
to be correct either for conditional outcome means or for conditional effects.

Note the asymmetry of the current strategy. We do not estimate the NDE by
treating the pseudo control sample with predicted $Y_{1M_{0}}$ and the pseudo
control sample with observed $Y_{0}$ as a two-arm trial. That is because these
predicted $Y_{1M_{0}}$ and observed $Y_{0}$ are on the same units. We could
instead use them to estimate the NDE as the mean of the individual effects,
which we did in section 4.1.

####  4.2.3 TE|psCadj: estimating TE with pseudo samples mimicking a two-way
randomized trial using covariate adjustment

We could use the same strategy to estimate TE, using the pseudo treated and
control samples. This estimator assumes the $\omega_{1}(C)$ and
$\omega_{0}(C)$ weights are consistent.

###  4.3 Combining to form estimators of the pair of effects

Apart from NDE&NIE|psxCadj (which is an estimator of both effects), other
estimators are for single effects and need to be combined. We combine each of
the two NDE estimators with an NIE estimator or a TE estimator. In addition to
the TE estimator TE|psCadj above, we add options by repurposing two reg|
estimators (fuYpred(ss) and fuYpred(ps)) as TE estimators.

Table LABEL:tab:effectEstimators thus includes four nonrobust estimators, one
estimator that is more robust but requires consistent $\omega_{1}(C)$ and
$\omega_{0}(C)$ weights, and one estimator that does not depend on any
specific estimation component being consistent. This last one, which combines
NDE|fuNDEpred(pxp0) and TE|fuYpred(ps), is very similar to the robust
estimator POs|fuYpred(ps)-fuY2pred(pxp0) from the outcome mean estimation
approach.

###  4.4 Data example application

We use the same weights estimated from section 3.4. For $Y_{1M_{0}}$
prediction (which is required for NDE|fuNDEpred estimators and NIE|psYpredCadj
estimators), we use the same $\text{E}[Y\mid C,M,A=1]$ model from section 3.4.

Model for $N\\!D\\!E_{0}$ given covariates. The NDE|fuNDEpred estimators
involve regressing the proxy of the individual $N\\!D\\!E_{0}$ (predicted
$Y_{1M_{0}}$ minus observed $Y_{0}$) on covariates. As the difference between
the two (predicted and observed) binary outcomes is bounded in the $-1$ to 1
interval, we use a simple transformation by adding 1 then dividing by 2 to map
to the 0 to 1 interval, and fit the regression model using the transformed
difference using logit link. The predictions are transformed back to the
original scale by multiplying by 2 then subtracting 1, before being averaged
to estimate $\text{NDE}_{0}$. (This is equivalent to the use of the tanh link
to estimate a function bounded in $[-1,1]$ [Wang and Tchetgen Tchetgen, 2018];
the transformation trick allows us to fit the model using standard software
with logit link.) The formula we use for this model is: transformed.diff
$\sim$ age + sex + edu + religion + att0 + rul0 + ns(sfc0,3) + drink0.

Working model for covariate adjustment. For all covariate adjustment
estimators, we use a simple logit model as the working model (for estimation
of all three effects, $\text{NDE}_{0},\text{NIE}_{1}$ and TE). The formula we
use for this working model is: drink $\sim$ arm + age + sex + edu + religion +
att0 + rul0 + sfc0 + drink0.

Results. Estimates from the six direct estimators of effects are presented in
Fig. 7 alongside effects from the POs| estimators.

##  5 Confidence interval estimation

Sections 3 and 4 cover point estimators of the NDE and NIE. We now turn our
attention to obtaining confidence intervals. We consider two options: using
asymptotic variance and bootstrapping.

We derive (see Technical Appendix) general formulas for the asymptotic
(sandwich) variance estimators for each of the estimators when all estimation
components (weights, conditional mediator density, conditional mean
outcome/effect) are based on (semi)parametric models – by applying Stefanski
and Boos’s (2002) calculus of M-estimators.

We use bootstrapping as a generic tool to obtain confidence intervals for all
the estimators. As our data example includes quite a few categorical
covariate/mediator variables, a challenge when using the simple resampling
bootstrap [Efron, 1979] is that some bootstrap samples do not cover all values
(and combinations of values) of those variables, resulting in predictors being
dropped from models and predictions being distorted. To avoid this problem we
use a continuous weights bootstrap [Xu et al., 2020] instead. With both
procedures, the making of a bootstrap sample can be seen as weighting the
observations of the original sample by a set of random weights that are
identically distributed: the resampling bootstrap uses integer weights drawn
from a uniform Multinomial distribution; the continuous weights bootstrap
draws weights from a continuous distribution. We use the version in Xu et al.
based on the uniform Dirichlet distribution, where the weights sum to sample
size $n$, have mean 1 and variance $(n-1)/(n+1)$. (For comparison, resampling
weights sum to $n$, have mean 1 and variance $(n-1)/n$.) Continuous weights
bootstrap samples retain all observations, thus they do not lose data
patterns. The confidence intervals for the data application in Fig. 7 are
obtained using this bootstrap procedure.

##  6 Discussion

In this paper we have shown how a range of estimators may be constructed based
on two strategies that are familiar to many who are involved in statistical
analyses (weighting and model-based prediction) and a simple way of combining
them (weighted models); this is the paper’s primary goal. The key ideas of
this exercise, which are not specific to natural effects but apply generally,
are (i) to use these tools flexibly to put together the different pieces of
the estimation puzzle, where the puzzle is defined by the identification
result of the effects of interest; and (ii) to induce robustness on pieces of
the puzzle by using weighted models. One thing to note is that we can correct
bias due to misspecified conditional mean models if we have correct weights,
but we cannot correct bias due to misspecified conditional density models in
the same way.

The secondary goal of the paper is to provide, for the specific estimand
considered (marginal natural (in)direct effects) a menu of estimation options.
The menu we have created is summarized in Tables LABEL:tab:regcrossEstimators
and LABEL:tab:effectEstimators. There are different ways the menu can be used.
As nonrobust estimators are simpler, one approach is to consider whether to
pick one of those estimators. Among the nonrobust options we do not recommend
the weighting-based estimators, as they are inefficient. Otherwise, we
recommend considering the set of estimation components (weights, mediator
density and outcome/effect mean models) required by each estimator and
deciding which set one is most feasible to implement well without severe
misspecification. The disadvantage of the simpler nonrobust estimators, of
course, is the lack of robustness. Another approach is to consider the more
robust estimators, looking at the components that must be correctly specified
(last column of the tables) as a way to rule out estimators that depend on
hard-to-model components. Again, we note that density estimation is generally
more challenging than mean estimation. Several of these estimators depend on
correct weights. While this cannot be determined, balance checking is a useful
tool. If the weights do not achieve near perfect balance, it is advisable not
to use these estimators whose validity hangs on the weights. And last but very
important, the menu includes two triply robust estimators, which do not depend
on any specific modeling component being correct. Note that for each of these,
we still need at least one of the three sets of estimation components to be
correct, and preferably more than one.

## References

  * Albert [2012] Jeffrey M Albert.  Distribution-free mediation analysis for nonlinear models with confounding.  _Epidemiology_ , 23(6):879–88, 2012.  doi: 10.1097/EDE.0b013e31826c2bb9.
  * Colantuoni and Rosenblum [2015] Elizabeth Colantuoni and Michael Rosenblum.  Leveraging prognostic baseline variables to gain precision in randomized trials.  _Statistics in Medicine_ , 34(18):2602–2617, 2015\.  doi: 10.1002/sim.6507.
  * Didelez et al. [2006] Vanessa Didelez, A Philip Dawid, and Sara Geneletti.  Direct and indirect effects of sequential treatments.  In _Proceedings of the 22nd Conference on Uncertainty in Artificial Intelligence_ , pages 138–146. AUAI Press, 2006.
  * Efron [1979] Bradley Efron.  Bootstrap methods: Another look at the jackknife.  _The Annals of Statistics_ , 11(1):1–26, 1979\.
  * Hainmueller [2012] Jens Hainmueller.  Entropy balancing for causal effects: A multivariate reweighting method to produce balanced samples in observational studies.  _Political Analysis_ , 20(1):25–46, 2012.  doi: 10.1093/pan/mpr025.
  * Holland [1986] Paul W Holland.  Statistics and causal inference.  _Journal of the American Statistical Association_ , 81(396):945, 1986.  doi: 10.2307/2289064.
  * Hong [2010] Guanglei Hong.  Ratio of mediator probability weighting for estimating natural direct and indirect effects.  In _Proceedings of the American Statistical Association, Biometrics Section_ , pages 2401–2415, 2010.
  * Hong et al. [2015] Guanglei Hong, Jonah Deutsch, and Heather D. Hill.  Ratio-of-mediator-probability weighting for causal mediation analysis in the presence of treatment-by-mediator interaction.  _Journal of Educational and Behavioral Statistics_ , 40(3):307–340, 2015.  doi: 10.3102/1076998615583902.
  * Huber [2014] Martin Huber.  Identifying causal mechanisms (primarily) based on inverse probability weighting.  _Journal of Applied Econometrics_ , 29(6):920–943, 2014.  ISSN 08837252.  doi: 10.1002/jae.2341.
  * Huling and Mak [2020] Jared D. Huling and Simon Mak.  Energy balancing of covariate distributions.  _arXiv_ , pages 1–68, 2020.
  * Imai and Ratkovic [2014] Kosuke Imai and Marc Ratkovic.  Covariate balancing propensity score.  _Journal of the Royal Statistical Society. Series B: Statistical Methodology_ , 76(1):243–263, 2014.
  * Imai et al. [2010a] Kosuke Imai, Luke Keele, and Dustin Tingley.  A general approach to causal mediation analysis.  _Psychological Methods_ , 15(4):309–34, 2010a.  doi: 10.1037/a0020761.
  * Imai et al. [2010b] Kosuke Imai, Luke Keele, and Teppei Yamamoto.  Identification, inference and sensitivity analysis for causal mediation effects.  _Statistical Science_ , 25(1):51–71, 2010b.  doi: 10.1214/10-STS321.
  * Kang and Schafer [2007] Joseph D. Y. Kang and Joseph L. Schafer.  Demystifying Double Robustness: A Comparison of Alternative Strategies for Estimating a Population Mean from Incomplete Data.  _Statistical Science_ , 22(4):523–539, 2007.  ISSN 0883-4237.  doi: 10.1214/07-sts227rej.
  * Koning et al. [2009] Ina M Koning, Wilma a M Vollebergh, Filip Smit, Jacqueline E E Verdurmen, Regina J J M Van Den Eijnden, Tom F M Ter Bogt, Håkan Stattin, and Rutger C M E Engels.  Preventing heavy alcohol use in adolescents (PAS): cluster randomized trial of a parent and student intervention offered separately and simultaneously.  _Addiction_ , 104(10):1669–78, 2009.  doi: 10.1111/j.1360-0443.2009.02677.x.
  * Koning et al. [2010] Ina M Koning, Regina J J M van den Eijnden, Rutger C M E Engels, Jacqueline E E Verdurmen, and Wilma a M Vollebergh.  Why target early adolescents and parents in alcohol prevention? The mediating effects of self-control, rules and attitudes about alcohol use.  _Addiction_ , 106(3):538–46, 2010.  doi: 10.1111/j.1360-0443.2010.03198.x.
  * Koning et al. [2011] Ina M. Koning, Regina J. Van Den Eijnden, Jacqueline E. Verdurmen, Rutger C. Engels, and Wilma a. Vollebergh.  Long-term effects of a parent and student intervention on alcohol use in adolescents: A cluster randomized controlled trial.  _American Journal of Preventive Medicine_ , 40:541–547, 2011.  doi: 10.1016/j.amepre.2010.12.030.
  * Lange et al. [2012] Theis Lange, Stijn Vansteelandt, and Maarten Bekaert.  A simple unified approach for estimating natural direct and indirect effects.  _American Journal of Epidemiology_ , 176(3):190–195, 2012.  doi: 10.1093/aje/kwr525.
  * Muthén and Asparouhov [2015] Bengt O Muthén and Tihomir Asparouhov.  Causal effects in mediation modeling: An introduction with applications to latent variables.  _Structural Equation Modeling_ , 22(1):12–23, 2015.  ISSN 15328007.  doi: 10.1080/10705511.2014.935843.
  * Nguyen et al. [2016] Trang Quynh Nguyen, Yenny Webb-Vargas, Ina M. Koning, and Elizabeth A. Stuart.  Causal mediation analysis with a binary outcome and multiple continuous or ordinal mediators: Simulations and application to an alcohol intervention.  _Structural Equation Modeling: A Multidisciplinary Journal_ , 23(3):368–383, 2016.  doi: 10.1080/10705511.2015.1062730.
  * Nguyen et al. [2020a] Trang Quynh Nguyen, Ian Schmid, Elizabeth L. Ogburn, and Elizabeth A. Stuart.  Clarifying causal mediation analysis for the applied researcher: Effect identification via three assumptions and five potential outcomes.  2020a.  URL http://arxiv.org/abs/2011.09537.
  * Nguyen et al. [2020b] Trang Quynh Nguyen, Ian Schmid, and Elizabeth A. Stuart.  Clarifying causal mediation analysis for the applied researcher: Defining effects based on what we want to learn.  _Psychological Methods_ , 2020b.  doi: 10.1037/met0000299.  URL http://arxiv.org/abs/1904.08515.
  * Pearl [2001] Judea Pearl.  Direct and indirect effects.  _Proceedings of the Seventeenth Conference on Uncertainty and Artificial Intelligence_ , pages 411–420, 2001.
  * Pearl [2012] Judea Pearl.  The causal mediation formula–a guide to the assessment of pathways and mechanisms.  _Prevention Science_ , 13(4):426–36, 2012.  ISSN 1573-6695.  doi: 10.1007/s11121-011-0270-1.
  * Robins et al. [2007] James Robins, Mariela Sued, Quanhong Lei-Gomez, and Andrea Rotnitzky.  Comment: Performance of double-robust estimators when “inverse probability” weights are highly variable.  _Statistical Science_ , 22(4):544–559, 2007.  doi: 10.1214/07-STS227D.
  * Robins and Greenland [1992] James M Robins and Sander Greenland.  Identifiability and exchangeability for direct and indirect effects.  _Epidemiology_ , 3(2):143–155, 1992.
  * Rosenbaum and Rubin [1983] Paul R Rosenbaum and Donald B Rubin.  The Central Role of the Propensity Score in Observational Studies for Causal Effects.  _Biometrika_ , 70(1):41, 1983.  doi: 10.2307/2335942.
  * Rubin [1974] Donald B Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of Educational Psychology_ , 66(5):688–701, 1974.  doi: 10.1037/h0037350.
  * Steen et al. [2017] Johan Steen, Tom Loeys, Beatrijs Moerkerke, and Stijn Vansteelandt.  Medflex: An R package for flexible mediation analysis using natural effect models.  _Journal of Statistical Software_ , 76(11), 2017.  doi: 10.18637/jss.v076.i11.
  * Stefanski and Boos [2002] Leonard A. Stefanski and Dennis D. Boos.  The calculus of M-estimation.  _The American Statistician_ , 56(1):29–38, 2002\.  doi: 10.1198/000313002753631330.
  * Steingrimsson et al. [2017] Jon Arni Steingrimsson, Daniel F. Hanley, and Michael Rosenblum.  Improving precision by adjusting for prognostic baseline variables in randomized trials with binary outcomes, without regression model assumptions.  _Contemporary Clinical Trials_ , 54:18–24, 2017.  doi: 10.1016/j.cct.2016.12.026.
  * Székely and Rizzo [2013] Gábor J. Székely and Maria L. Rizzo.  Energy statistics: A class of statistics based on distances.  _Journal of Statistical Planning and Inference_ , 143(8):1249–1272, 2013.  doi: 10.1016/j.jspi.2013.03.018.
  * Tchetgen Tchetgen [2013] Eric J. Tchetgen Tchetgen.  Inverse odds ratio-weighted estimation for causal mediation analysis.  _Statistics in Medicine_ , 32(26):4567–4580, 2013\.  ISSN 02776715.  doi: 10.1002/sim.5864.
  * Tchetgen Tchetgen and Shpitser [2012] Eric J. Tchetgen Tchetgen and Ilya Shpitser.  Semiparametric theory for causal mediation analysis: Efficiency bounds, multiple robustness and sensitivity analysis.  _The Annals of Statistics_ , 40(3):1816–1845, 2012.  ISSN 0090-5364.  doi: 10.1214/12-AOS990.
  * Tchetgen Tchetgen and Shpitser [2014] Eric J. Tchetgen Tchetgen and Ilya Shpitser.  Estimation of a semiparametric natural direct effect model incorporating baseline covariates.  _Biometrika_ , 101(4):849–864, 2014.  ISSN 14643510.  doi: 10.1093/biomet/asu044.
  * Tingley et al. [2014] Dustin Tingley, Teppei Yamamoto, Kentaro Hirose, Luke Keele, and Kosuke Imai.  mediation: R package for causal mediation analysis.  _Journal of Statistical Software_ , 59(5):1–38, 2014.  ISSN 15487660.  doi: 10.18637/jss.v059.i05.
  * Valeri and VanderWeele [2013] Linda Valeri and Tyler J VanderWeele.  Mediation analysis allowing for exposure–mediator interactions and causal interpretation: Theoretical assumptions and implementation with SAS and SPSS macros.  _Psychological methods_ , 18(2):137–150, 2013\.  ISSN 1939-1463.  doi: 10.1037/a0031034.
  * van der Laan and Rose [2011] Mark J. van der Laan and Sherri Rose.  _Targeted Learning: Causal Inference for Observational and Experimental Data_.  Springer Series in Statistics. Springer New York, 2011.  ISBN 9781441997821.  URL https://books.google.com/books?id=RGnSX5aCAgQC.
  * VanderWeele and Vansteelandt [2009] Tyler J VanderWeele and Stijn Vansteelandt.  Conceptual issues concerning mediation, interventions and composition.  _Statistics and its Interface_ , 2:457–468, 2009.
  * VanderWeele and Vansteelandt [2010] Tyler J VanderWeele and Stijn Vansteelandt.  Odds ratios for mediation analysis for a dichotomous outcome.  _American Journal of Epidemiology_ , 172(12):1339–1348, 2010.  doi: 10.1093/aje/kwq332.
  * VanderWeele and Vansteelandt [2013] Tyler J VanderWeele and Stijn Vansteelandt.  Mediation analysis with multiple mediators.  _Epidemiologic Methods_ , 2(1):95–115, 2013.  doi: 10.1515/em-2012-0010.
  * VanderWeele et al. [2014] Tyler J VanderWeele, Stijn Vansteelandt, and James M Robins.  Effect decomposition in the presence of an exposure-induced mediator-outcome confounder.  _Epidemiology_ , 25(2):300–6, 2014.  ISSN 1531-5487.  doi: 10.1097/EDE.0000000000000034.
  * Vansteelandt and Keiding [2011] Stijn Vansteelandt and Niels Keiding.  Invited commentary: G-computation–Lost in translation?  _American Journal of Epidemiology_ , 173(7):739–742, 2011.  doi: 10.1093/aje/kwq474.
  * Vansteelandt et al. [2012] Stijn Vansteelandt, Maarten Bekaert, and Theis Lange.  Imputation strategies for the estimation of natural direct and indirect effects.  _Epidemiologic Methods_ , 1(1):7, 2012.  doi: 10.1515/2161-962X.1014.
  * Wang et al. [2019] Bingkai Wang, Elizabeth L. Ogburn, and Michael Rosenblum.  Analysis of Covariance in Randomized Trials: More Precision, Less Conditional Bias, and Valid Confidence Intervals, Without Model Assumptions.  _Biometrics_ , 75(4):1391–1400, 2019.  doi: 10.1111/biom.13062.
  * Wang and Tchetgen Tchetgen [2018] Linbo Wang and Eric Tchetgen Tchetgen.  Bounded, efficient and multiply robust estimation of average treatment effects using instrumental variables.  _Journal of the Royal Statistical Society. Series B: Statistical Methodology_ , 80(3):531–550, 2018.  doi: 10.1111/rssb.12262.
  * Xu et al. [2020] Li Xu, Chris Gotwalt, Yili Hong, Caleb B. King, and William Q. Meeker.  Applications of the fractional-random-weight bootstrap.  _American Statistician_ , 1305:1–32, 2020.  doi: 10.1080/00031305.2020.1731599.
  * Zheng and van der Laan [2012] Wenjing Zheng and Mark J van der Laan.  Targeted maximum likelihood estimation of natural direct effects.  _The International Journal of Biostatistics_ , 8(1), 2012\.  doi: 10.2202/1557-4679.1361.

See pages - of TechAppendix20210510.pdf