# A QCD motivated model for soft processes

A. Kormilitzin  E. Levin

###### Abstract

In this talk we give a brief description of a QCD motivated model for both
hard and soft interactions at high energies. In this model the long distance
behaviour of the scattering amplitude is determined by the dipole scattering
amplitude in the saturation domain. All phenomenological parameters for
dipole-proton interaction were fitted from the deep inelastic scattering data
and the soft processes are described with only one new parameter, related to
the wave function of hadron. It turns out that we do not need to introduce the
so called soft Pomeron that has been used in high energy phenomenology for
four decades.

###### Keywords:

High density QCD, saturation, Pomeron structure, diffraction

###### :

13.8.5.-t, 13.85.Hd,11.55.-m,11.55.Bq

##  1 Main ideas and theoretical input

This talk is the brief description of our attempt to find a self-consistent
theoretical approach to soft (long distance) interaction at high energy (see
our paper of Ref.KOLE  for full presentation). Our main idea is that there is
no other dimemsionful scales except the saturation scale for high energy
interaction in the entire kinematic region. This idea looks strange since the
traditional approach to soft high energy interaction is based on the
phenomenological soft Pomeron, whose typical dimensionful scale is given by
the slope of the Pomeron trajectory ($\alpha_{P}\,\approx 0.25\,GeV^{-2}$). In
simple words, we would like to replace the soft Pomeron exchange in the
scattering amplitude by the QCD scattering amplitude in the saturation region.
$\alpha_{P}=0$ for this amplitude and the saturation momentum ( $Q_{s}$) is
the only dimensionful scale for it. Since this amplitude governs both short
distances and long distances processes, all needed parameters in a such
amplitude can be found from the DIS on the contrary to soft phenomenology in
which all parameters are extracted from the soft scattering. The idea is not
new and it is scattered in a number of reviews and talks (see references in
Ref.KOLE ). However the first practical realization was suggested in Ref.BATAV
with an encouraging result: it is possible to describe the data on total cross
section in the framework of such ideas. In this talk as well as in our paper
KOLE  we wish to check these ideas against the full set of experimental data
in the region of accessible energies and even give some predictions for the
LHC energy range. Our approach is based on two theoretical inputs.

Factorization of short and long distances: Using the dispersion relation we
can prove (see Ref.KOLE  that at least at large impact parameters the
scattering amplitude for the BFKL Pomeron can be written in the form

$A_{BFKL}(Q^{2},x;b)\,=\,\stackrel{{\scriptstyle\mbox{ \small short
distances}}}{{A_{BFKL}(Q^{2},x;t=0)}}\,\times\,{\stackrel{{\scriptstyle\mbox{\small
long distances}}}{{S(b)}}}$ (1)

The scattering amplitude of dipole (x,y) near to the saturation boundary  has
the formLMP :

$N(Y-Y_{0};x,y)\,=\,1-e^{-\int\,d^{2}x^{\prime}\,d^{2}y^{\prime}\,\,P(Y-Y_{0};x,y;x^{\prime},y^{\prime})\gamma(Y_{0};x^{\prime},y^{\prime})}$
(2)

where $\gamma(Y_{0};x^{\prime},y^{\prime})$ is the scattering amplitude of
dipole $(x^{\prime},y^{\prime})$ with the target at low energy and
$P(Y-Y_{0};x,y;x^{\prime},y^{\prime})$ is the Green’s function of the BFKL
Pomeron.

##  2 The model

Using these two theory inputs we built our model assuming

  * •

$\gamma_{mod}(Y_{0};x^{\prime},y^{\prime})\,=\,\,\,\frac{\pi^{2}\,\alpha_{S}(\mu^{2}_{0})}{3}\,r^{2}\,x_{0}G^{DGLAP}(x_{0},\mu^{2}_{0})\,S(b)$;

  * •

$S(b)\,=\,\frac{2}{\pi\,R^{2}}(\frac{\sqrt{8}\,b}{R})\,K_{1}(\frac{\sqrt{8}\,b}{R})$\-
the electromagnetic proton form factor;

  * •

Instead of the BFKL Pomeron in Eq.2 we use the solution to the DGLAP equation;

  * •

The scale of hardness in the DGLAP evolution ia chosen in the form:
$\mu^{2}\,=\,\mu^{2}_{0}\,\,+\,\,\frac{C}{r^{2}}$.

Finally,

$N(Y-Y_{0};x,y)\,=\,1-\exp\left(-\stackrel{{\scriptstyle\mbox{ \small short
distances}}}{{\frac{\pi^{2}\,\alpha_{S}(\mu^{2})}{3}\,r^{2}\,xG^{DGLAP}(x,\mu)}}\,\times\,{\stackrel{{\scriptstyle\mbox{
\small long distances}}}{{S(b)}}}\right)$ (3)

All parameters in Eq.3 we found from fitting the DIS data with
$\chi^{2}$/d.o.f. =1.02. The formulae for calculation of the physical
observables are simple and can be found in Ref. KOLE  and the structure of
them is seen from the expression for the total cross section for meson proton
scattering: $\sigma_{tot}\,=\,\int d^{2}bd^{2}r|\Psi(r)|^{2}\,N(Y-Y_{0},r)$
where $\Psi$ is the wave function of a meson which consists of one colourless
dipole with size $r$, and $N$ is dipole scattering amplitude. Therefore, we
need to specify the wave hadronic wave function in terms of t dipoles. As have
been mentioned, we assume that mesons consist of one dipole, for baryon we
have two possibilities for the dipole content: two dipoles or three dipoles.
We used both in our comparison with the experiment. For one dipole in the
hadron we used the simplest Gaussian parameterization, namely
$\Psi(r)\,=\frac{1}{\pi S}\,\exp(-\frac{r^{2}_{i}}{S})$ where $S$ was chosen
from comparison with the experiment. This is the only parameter that we found
from the data on soft scattering. The second input that we need to calculate
the scattering amplitude is to define the energy variable $x$ (or $Y-Y_{0}$).
We did this in the following way, using our the only dimensionful scale
(saturation momentum $Q_{s}$):

$x_{soft}\,=\,\frac{Q^{2}+Q^{2}_{s}(x_{soft})}{s}\,\,\stackrel{{\scriptstyle
Q^{2}\to
0}}{{\rightarrow}}\,\,\frac{Q^{2}_{s}(x_{soft})}{s};\,\,\,\,x_{soft}\,=\,\frac{Q^{2}+Q^{2}_{s}(x_{soft})}{s}\,\,\stackrel{{\scriptstyle
Q^{2}\gg Q^{2}_{s}}}{{\rightarrow}}\,\,x_{Bjorken};$ (4)

##  3 Comparison with the data

Using the model we obtain the description of the experimental data which is
not too bad. We consider that it is much better than we could hope with such a
simple model.

|  |
---|---|---
$\sigma_{tot}(\pi-p)$ | $\sigma_{tot}(K-p)$ | $\sigma_{tot}(p-\bar{p})$
Figure 1:  |
---|---
$\sigma_{el}(p-\bar{p})$ | $d\sigma_{el}(p\bar{p})/dt$
Figure 2:  Model | $\sigma_{tot}$ | $\sigma_{el}$ | $\sigma_{diff}$ | $B_{el}$ | $<|S^{2}|>$
---|---|---|---|---|---
| mb | mb | mb | $GeV^{-2}$ |
Our A (DIS) | 83.0 | 23.54 | 10 | 16.67 | (0.24-0.89)%
Our B (Soft) | 101.3 | 28.84 | 10.5 | 18.4 | (0.24-0.57)%
GLMGLM1  | 110.5 | 25.3 | 11.6 | 20.5 | (0.7 -2)%
GLMMGLM2  | 91.7 | 20.9 | 11.8 | 17.3 | 0.21%
RMKRMK  | 88.0 (86.3) | 20.1 (18.1) | 13.3 (16.1) | 19 | (1.2 - 3.2)%
Table 1: Predictions for the LHC

One can see that our model reproduces quite well the energy and $t$ dependence
of total and elastic cross sections. The most essential test for the model,
however, is the dependence of the elastic slope versus energy. Indeed, in
traditional high energy phenomenology this dependence was related to the value
of $\alpha^{\prime}_{P}\approx 0.25\,GeV^{-2}$ which in our model is equal to
zero. Nevertheless, in spite of the fact that input
$\alpha^{\prime}_{P}\,=\,0$ the shadowing corrections ,taken into account in
our master formula of Eq.3, generate effective $\alpha^{\prime}_{p}\approx
0.2\,GeV^{-2}$. For the single diffraction cross section we take into account
the emission of an additional gluon. The low curve in Fig. 3 corresponds to
diffraction production of one dipole while the resulting curve includes the
extra gluon emission. In all figures we plotted our calculations for two sets
of parameterizations: set A (DIS) which gives very good description of DIS
data with $\chi^{2}$/d.o.f. = 1.02; and set B(soft) which leads to worse
description of DIS data ($\chi^{2}/d.o.f$\- 3.6) but predicts the soft
observables closer to the experimental data. Notations 2D and 3D specify the
model for nucleon with two or three dipoles.

|
---|---
Figure 3:

##  4 Conclusions

In the table we show our prediction for the LHC range of energy with
comparison with other models. The most pronounced difference between our
estimates and RMK modelRMK  in the value of the survival probability for the
central Higgs production. It is worth mentioning, that in our approach the
small value of $<|S^{2}|>$ stems from Good-Walker mechanism but not from the
enhanced diagrams as in GLMM model GLM2 .

We will be happy if you, thinking about Pomeron, would remember that (i) there
is a possibility to describe the data without introducing soft Pomeron; (ii)
the idea about the saturation momentum as the only dimensional parameter of
the strong interaction works; and (iii) all parameters for strong interaction
can be found from DIS.

We are grateful to Errol Gotsman, Uri Maor and Alex Palatnik for useful
discussions on the subject. This research was supported in part by the Israel
Science Foundation, founded by the Israeli Academy of Science and Humanities,
by BSF grant $\\#$ 20004019 and by a grant from Israel Ministry of Science,
Culture and Sport and the Foundation for Basic Research of the Russian
Federation.

## References

  * (1) A. Kormilitzin and E. Levin, “Soft processes at high energy without soft Pomeron: a QCD motivated model,” arXiv:0809.3886 [hep-ph].
  * (2) J. Bartels, E. Gotsman, E. Levin, M. Lublinsky and U. Maor, Phys. Rev. D 68 (2003) 054008 [arXiv:hep-ph/0304166]; Phys. Lett. B 556 (2003) 114 [arXiv:hep-ph/0212284].
  * (3) E. Levin, J. Miller and A. Prygarin, Nucl.Phys. A (in press); arXiv:0706.2944 [hep-ph].
  * (4) E. Gotsman, E. Levin and U. Maor, ”A Soft Interaction Model at Ultra High Energies: Amplitudes, Cross Sections and Survival Probabilities”, arXiv:0708.1506 [hep-ph]
  * (5) E. Gotsman, E. Levin, U. Maor and J. Miller, Eur. Phys. J. C,( in press); arXiv:0805.2799 [hep-ph].
  * (6) M. G. Ryskin, A. D. Martin and V. A. Khoze, Eur. Phys. J. C54 (2008) 199 [arXiv:0710.2494 [hep-ph]].


# Estimating Heterogeneous Causal Mediation Effects with Bayesian Decision
Tree Ensembles

Angela Ting and Antonio R. Linero
_Department of Statistics and Data Sciences_
_University of Texas at Austin_

###### Abstract

The causal inference literature has increasingly recognized that explicitly
targeting treatment effect heterogeneity can lead to improved scientific
understanding and policy recommendations. Towards the same ends, studying the
causal pathway connecting the treatment to the outcome can be also useful.
This paper addresses these problems in the context of _causal mediation
analysis_. We introduce a varying coefficient model based on Bayesian additive
regression trees to identify and regularize heterogeneous causal mediation
effects; analogously with linear structural equation models, these effects
correspond to covariate-dependent products of coefficients. We show that, even
on large datasets with few covariates, LSEMs can produce highly unstable
estimates of the conditional average direct and indirect effects, while our
_Bayesian causal mediation forests_ model produces estimates that are stable.
We find that our approach is conservative, with effect estimates “shrunk
towards homogeneity.” We examine the salient properties of our method using
both data from the Medical Expenditure Panel Survey and empirically-grounded
simulated data. Finally, we show how our model can be combined with posterior
summarization strategies to identify interesting subgroups and interpret the
model fit.

##  1 Introduction

Estimation of heterogeneous causal effects from observational data is a topic
of fundamental importance, with applications in personalized medicine
(Obermeyer and Emanuel,, 2016), policy recommendation (Athey,, 2017), and
social science (Yeager et al.,, 2019). A question of great recent interest in
the causal inference literature is how best to leverage state-of-the-art
prediction algorithms developed in the machine learning community to estimate
heterogeneous treatment effects (Künzel et al.,, 2019; Nie and Wager,, 2021;
Hahn et al.,, 2020). Much of this literature has focused on the question of
how best to modify the estimation strategies used in the machine learning
literature to be appropriate for inferring heterogeneous causal effects.

A complementary approach to making better policy recommendations is to learn
how a treatment of interest influences the outcome via its effects on
downstream variables that are themselves causally linked to the outcome; this
is referred to as _causal mediation analysis_ and the intermediate variables
are referred to as _mediators_ (Robins and Greenland,, 1992; Pearl,, 2001;
Rubin,, 2004). In addition to providing a sharper understanding of the causal
mechanisms at play, we will see that causal mediation analysis can in some
cases increase our power to detect causal effects. Similar questions about how
to effectively leverage predictive algorithms have emerged in this field, with
much of the focus on estimating _average_ , rather than heterogeneous,
mediation effects (Farbmacher et al.,, 2022; Linero and Zhang,, 2022; Zheng
and van der Laan,, 2012; Tchetgen and Shpitser,, 2012; Kim et al.,, 2017).

To the best of our knowledge, there has been limited work at the intersection
of these two settings, i.e., where one is interested in estimating treatment
effect heterogeneity at the level of direct and indirect causal mediation
effects using machine learning. The issue of estimating heterogeneity in
treatment effects in the context of mediation analysis is referred to as
_moderated mediation_ (Muller et al.,, 2005). This topic has garnered
significant attention in the social science literature, often utilizing linear
structural equation modeling (LSEM). For example, (Preacher et al.,, 2007) and
(Kershaw et al.,, 2010) applied moderated mediation using LSEMs to problems in
education and health psychology, respectively.

Estimating heterogeneous mediation effects in a nonparametric manner is a
challenging task that relies on both strong assumptions regarding confounding
and requires large amounts of data to reliably estimate the causal effects.
There are important challenges in this context that need to be addressed,
including: (i) determining how to properly regularize both the nuisance
parameters and parameters of interst to ensure sensible results; (ii)
developing methods to summarize the results of black-box fitting procedures in
a meaningful way; and (iii) establishing reliable techniques to identify
subgroups for which there is evidence of moderated mediation and to determine
which variables are acting as effect modifiers.

This paper proposes a two-layer extension of the Bayesian causal forests (BCF)
algorithm for estimating heterogeneous mediation effects, which combines a
standard BCF model for the mediator with a varying coefficient BART model for
the outcome (Hahn et al.,, 2020; Deshpande et al.,, 2020). Our approach is
motivated by the strong performance of BCFs in causal inference competitions
and in practical applications (Dorie et al.,, 2019). Our approach directly
parameterizes the models in terms of the direct and indirect effects of the
treatment on the outcome. This allows us to “shrink towards homogeneity,”
stabilizing the estimation of the mediation effects. Our approach performs
extremely well in regimes where treatment effects are nearly homogeneous, with
small root-mean squared errors for individual-level mediation effects and
credible intervals that attain close to the nominal rate of coverage for most
individuals. Hence, our proposed approach provides a powerful tool for
estimating heterogeneous mediation effects.

###  1.1 The Medical Expenditure Panel Survey

The Medical Expenditure Panel Survey (MEPS) is an ongoing large-scale survey
administered by the Agency for Healthcare Research and Quality that aims to
measure the healthcare system’s use by patients, hospitals, and insurance
companies. To demonstrate our proposed methodology, we employ the MEPS to
investigate the health consequences of smoking. Specifically, we aim to answer
the following questions: (i) does smoking have a causal impact on healthcare
expenses overall? (ii) to what extent is this impact mediated (or not) by
smoking’s effect on overall health? and (iii) are there any moderating
variables that affect the association between smoking and medical
expenditures?

In Section 4, we present an analysis of this dataset which yields a surprising
finding: the total causal effect of smoking on medical expenditures can be
masked by instability resulting from the estimation of the direct effect of
smoking on healthcare costs. Although one might intuitively assume that the
effect of smoking on medical expenditures is fully mediated by its impact on
health, our analysis under sequential ignorability shows that the estimated
direct effect of smoking on expenses is negative and largely counteracts the
positive indirect effect of smoking on expenses; this direct effect is likely
due to additional variables that we have not incorporated in the analysis.
Additionally, we identify several variables, with age being the most
important, that moderate the indirect effect of smoking on expenditures.

###  1.2 Outline

In Section 2 we review the potential outcomes framework for mediation, the
sequential ignorability assumption, the Bayesian additive regression trees
(BART) framework, and Bayesian causal forests (BCFs). In Section 3 we define
our Bayesian causal mediation forests model, and show how to use it to stably
estimate the direct and indirect effects. In Section 4 we use our methodology
to analyze data from the MEPS data to study mediation effect heterogeneity in
the effect of smoking on health care expenditures as mediated by the effect of
smoking on health, and conduct an empirically-designed simulation study to
show that our method performs well in terms of coverage and estimation error
for estimating both average and conditional average mediation effects. We
conclude in Section 5 with a discussion and possible extensions. Computational
detials and further simulation results are given in the supplementary
material.

##  2 Review of Mediation Analysis and BART

###  2.1 Overview of Mediation Analysis

Mediation refers to the process through which a treatment ($A$) influences an
outcome ($Y$) by acting through an intermediate _mediator_ variable ($M$),
which occurs between the treatment and the outcome; a graphical representation
is given in Figure 1. For example, let us consider the question of whether
smoking affects medical expenditures directly and indirectly through its
effect on health. Here, smoking status is a binary treatment ($A$), and the
outcome of interest is the logarithm of medical expenditure ($Y$). Our aim is
to break down the effect of smoking on medical expenditures into a _direct
effect_ of smoking and an _indirect effect_ that is mediated by smoking’s
effect on overall health (measured as an individual’s self-perceived quality
of health). A natural hypothesis is that smoking does not directly cause
higher medical expenditures but rather does so by reducing a person’s overall
health. Health is on the causal path between the treatment (smoking) and the
outcome (medical expenditures) and hence is a mediator.

Figure 1: A schematic representation of a treatment $A$, mediator $M$, outcome
$Y$, and confounders/effect modifiers $X$. Arrows depict the direction of
causality.

Mediation analysis has been applied in many scientific fields including
epidemiology, medicine, economics, and the social sciences (MacKinnon and
Dwyer,, 1993; Rubin,, 2004; MacKinnon,, 2008; Albert,, 2008; Imai et al.,,
2010; VanderWeele,, 2016). Much of this literature has focused on structural
equation models (SEMs) to quantify mediation effects as products of
coefficients in parametric models. In particular, linear structural equation
models (LSEMs) have been widely used (Baron and Kenny,, 1986; MacKinnon and
Dwyer,, 1993; MacKinnon,, 2008).

LSEMs have a major limitation in that the identification of the mediation
effects is tied to the choice and correct specification of a particular
parameteric model, limiting their applicability. To address this limitation,
Imai et al., (2010) proposed a nonparametric approach based on _potential
outcomes_ (Rubin,, 2004, 1974) that allows for the identification of average
causal mediation effects under the assumption of _sequential ignorability_.
This assumption states that the treatment is independent of all potential
values of the outcome and mediator given the covariates, and the observed
mediator is independent of all potential outcomes given the observed treatment
and covariates. By avoiding parametric assumptions, this framework provides a
general estimation procedure that is agnostic to the choice of model for the
outcome and mediator, making it applicable in a wide range of settings.

For individuals $i=1,...,n$ and treatment $a\in\\{0,1\\}$, define the
potential outcome $M_{i}(a)$ as the value of the mediator that would have been
observed had the individual received treatment $a$. Note that for each
individual, only one of $M_{i}(0)$ or $M_{i}(1)$ is actually observed. For
treated individuals ($A_{i}=1$), $M_{i}(0)$ is a _counterfactual_ , i.e., the
value of the mediator that would have been observed had the individual been
untreated instead. Similarly, the potential outcome $Y_{i}(a,m)$ is the value
of the outcome that would have been observed had the individual received
treatment $a$ and had a mediator at level $m$. For example,
$Y_{i}\\{0,M_{i}(1)\\}$ is the value of the outcome that would have been
observed if the individual was not treated and had a value of the mediator at
the same level they would have had if they were treated. We link the potential
outcomes to the observed data through the _consistency_ assumption, which
states that we observe the mediator $M_{i}=M_{i}(A_{i})$ and the outcome
$Y_{i}=Y_{i}\\{A_{i},M_{i}(A_{i})\\}$. Because the values of $Y_{i}$ and
$M_{i}$ are defined only in terms of the treatment $a$ potentially received by
individual $i$ (and not on the treatment received by other individuals), this
notation also implicitly states that there is no interference between units,
which is known as the _Stable Unit Treatment Value_ (SUTVA) assumption.

Using these potential outcomes, we can define the causal estimates of
interest. In causal mediation analysis, we are particularly interested in
estimating the natural direct and natural indirect effects (Pearl,, 2001;
Robins and Greenland,, 1992). The natural direct effect is defined as

$\displaystyle\zeta_{a}=E[Y_{i}\\{1,M_{i}(a)\\}-Y_{i}\\{0,M_{i}(a)\\}]$ (1)

and the natural indirect effect is defined as

$\displaystyle\delta_{a}=E[Y_{i}\\{a,M_{i}(1)\\}-Y_{i}\\{a,M_{i}(0)\\}].$ (2)

The natural direct effect isolates the effect of the treatment while keeping
the potential mediator fixed, and can be interpreted as the effect that the
treatment has directly on the outcome $Y_{i}$. Conversely, the natural
indirect effect isolates the effect of the potential mediator in response to
different treatment values while keeping the treatment fixed, and can be
interpreted as the effect that the treatment has indirectly on the outcome
$Y_{i}$ through the mediator $M_{i}$. The total effect of the treatment on the
outcome is a sum of the direct and indirect effects, and can be defined as

$\displaystyle\tau=\zeta_{0}+\delta_{1}=\zeta_{1}+\delta_{0}=E[Y_{i}\\{1,M_{i}(1)\\}-Y_{i}\\{0,M_{i}(0)\\}].$
(3)

We can similarly define _conditional average_ variants of both the direct and
indirect effects as

$\displaystyle\begin{split}\zeta_{a}(x)&=E[Y_{i}\\{1,M_{i}(a)\\}-Y_{i}\\{0,M_{i}(a)\\}\mid
X_{i}=x]\qquad\text{and}\qquad\\\
\delta_{a}(x)&=E[Y_{i}\\{a,M_{i}(1)\\}-Y_{i}\\{a,M_{i}(0)\\}\mid
X_{i}=x].\end{split}$ (4)

Most of our attention will be on the conditional average direct and indirect
effects, as defined in (4).

###  2.2 Assumptions

Let the statement $[A\perp\\!\\!\\!\perp B\mid C=c]$ mean that $A$ is
conditionally independent of $B$ given that $C=c$, let $\mathcal{X}$ denote
the sample space of $X_{i}$, and let $\mathcal{M}$ denote the sample space of
$M_{i}$. Following Imai et al., (2010), we make the following _sequential
ignorability_ assumption throughout, allowing for the identification of the
direct and indirect effects.

SI1

$\\{Y_{i}(a^{\prime},m),M_{i}(a)\\}\perp\\!\\!\\!\perp A_{i}\mid X_{i}=x$ for
$a,a^{\prime}=0,1$ and all $x\in\mathcal{X}$.

SI2

$Y_{i}(a^{\prime},m)\perp\\!\\!\\!\perp M_{i}(a)\mid A_{i}=a,X_{i}=x$ for
$a,a^{\prime}=0,1$ and all $x\in\mathcal{X}$.

SI3

$\text{Pr}(A_{i}=a\mid X_{i}=x)>0$ and $f\\{M_{i}(a)=m\mid
A_{i}=a,X_{i}=x\\}>0$ for $a=0,1$ and all $x\in\mathcal{X}$ and
$m\in\mathcal{M}$.

The first assumption states that, given the covariates, the treatment
assignment is ignorable, i.e., it is independent of potential outcomes and
potential mediators. This assumption is automatically satisfied when
individuals are randomly assigned to treatment and control groups, but is not
guaranteed to hold in observational studies, in which case researchers often
collect as many pre-treatment confounders as possible so that treatment
assignment ignorability is plausible after the differences in covariates
between treatment groups are accounted for. The second assumption states that,
given the observed treatment and covariates, the mediator is ignorable, i.e.,
it is independent of potential outcomes. This assumption, however, is not
guaranteed to hold even in randomized experiments. In general, it cannot be
directly tested from the data. The third assumption is a positivity assumption
for the treatment and mediator, stating that the probability of receiving the
treatment and control should be nonzero.

Under SI1–SI3, we can identify the distribution of any counterfactual outcome
$Y_{i}\\{a^{\prime},M_{i}(a)\\}$ nonparametrically as

$\displaystyle\begin{split}&f(Y_{i}\\{a,M_{i}(a^{\prime})\\}=y\mid X_{i}=x)\\\
&=\int_{\mathcal{M}}f(Y_{i}=y\mid M_{i}=m,A_{i}=a,X_{i}=x)\,f(M_{i}=m\mid
A_{i}=a,X_{i}=x)\ dm\end{split}$ (5)

for any $x\in\mathcal{X}$ and $a,a^{\prime}=0,1$ (Imai et al.,, 2010, Theorem
1). This allows us to make inferences about unobserved counterfactuals (left-
hand side) using observed outcomes and mediators (right-hand side). Moreover,
(5) is not dependent on a specific parametric model, and so can be applied to
flexible (nonparametric) models.

###  2.3 A Review of Bayesian Additive Regression Trees

We will use the Bayesian Additive Regression Trees (BART) model proposed by
Chipman et al., (2010). Consider an unknown function $r$ that predicts an
output $Y_{i}$ using a vector of inputs $X_{i}$

$\displaystyle Y_{i}=r(X_{i})+\epsilon_{i},\quad\epsilon_{i}\sim
N(0,\sigma^{2}).$ (6)

BART models $r(x)$ as a sum of $m$ regression trees
$r(x)=\sum_{j=1}^{m}g(x;T_{j},M_{j})$ where $T_{j}$ is a binary decision tree
consisting of interior node decision rules as well as a set of terminal nodes
and $M_{j}=\\{\mu_{j1},\ldots,\mu_{jb_{j}}\\}$ is a set of parameter values
associated with each of the $b_{j}$ terminal nodes of tree $T_{j}$. Each $x$
is associated with a single terminal node $k$ of $T_{j}$ and is then assigned
the value $g(x;T_{j},M_{j})=\mu_{jk}$. Under (6), $E(Y_{i}\mid X_{i}=x)$
equals the sum of all the terminal node $\mu_{jk}$’s assigned to $x$ by the
$g(x;T_{j},M_{j})$’s. For a comprehensive review of BART and its applications,
see Hill et al., (2020).

To apply BART it is necessary to specify a prior distribution over all the
parameters of the sum-of-trees model, i.e., $(T_{j},M_{j})$ for
$j=1,\ldots,m$. This prior should regularize the fit by keeping individual
tree effects from being disproportionately influential. The prior consists of
two components: a prior for each tree $T_{j}$ and a prior on the terminal
nodes $M_{j}\mid T_{j}$ where
$\pi(T_{j},M_{j})=\pi_{T}(T_{j})\pi_{M}(M_{j}\mid T_{j})$. The BART model then
sets $(T_{j},M_{j})\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\pi(T,M)$.

The prior $\pi_{T}$ is determined by three variables: (i) the probability that
a given node is an interior node, (ii) the distribution of the splitting
variable assignments at each interior node, and (iii) the distribution of the
splitting rule assignment in each interior node conditional on the splitting
variable. For (i), the probability that a node at depth $d$ is an interior
node is

$\displaystyle\alpha(1+d)^{-\beta},\quad\alpha\in(0,1),\beta\in[0,\infty)$ (7)

with $\alpha=0.95$ and $\beta=2$ being a default that favors small trees. For
(ii) and (iii), the distribution of the splitting variable assignments at each
interior node and the distribution of the splitting rule assignment in each
interior node conditional on the splitting variable are both given a uniform
prior.

For the prior on the terminal nodes, we assume $\pi(M_{j}\mid
T_{j})=\prod_{k=1}^{b_{j}}\pi_{\mu}(\mu_{jk})$. To specify $\sigma_{\mu}$, in
this paper we first shift and rescale the $Y_{i}$’s so that $Y_{i}$ has mean
$0$ and variance $1$. We then use the prior

$\displaystyle\pi_{\mu}(\mu_{jk})=N(\mu_{jk}\mid
0,\sigma^{2}_{\mu})\qquad\text{where}\qquad\sigma_{\mu}=\frac{3}{k\sqrt{m}}$

for a suitable value of $k$, with default $k=2$. Note that this prior shrinks
the terminal node values $\mu_{jk}$ towards zero and applies greater shrinkage
as the number of trees $m$ is increased, ensuring that each tree is a weak
learner in the ensemble of trees.

###  2.4 Bayesian Decision Tree Ensembles for Causal Inference

BART has been seen to perform particularly well in causal inference problems
for inferring heterogeneous and average treatment effects (Hill,, 2011;
Wendling et al.,, 2018; Dorie et al.,, 2019). For an outcome $Y_{i}$, binary
treatment $A_{i}$, and confounder/modifier variables $X_{i}$, Hill, (2011)
proposes the model

$\displaystyle Y_{i}(a)=\mu(X_{i},a)+\epsilon_{i},\quad\epsilon_{i}\sim
N(0,\sigma^{2}).$ (8)

The effect of receiving the treatment is therefore given by

$\displaystyle E\\{Y_{i}(1)-Y_{i}(0)\mid
X_{i}=x\\}=\tau(x)=\mu(x,1)-\mu(x,0).$

Given BART’s strong predictive performance, Hill, (2011) suggests using a BART
prior for $\mu(\cdot,\cdot)$ to flexibly model the outcome and hence obtain
flexible treatment effect estimates.

Hahn et al., (2020) note that successful predictive modeling depends largely
on careful regularization, and extend the work of Hill, (2011) by noting two
shortcomings of the model (8): first, the correlation between the propensity
score and $\mu(x,a)$ can induce _regularization induced confounding_ (RIC),
leading to highly biased causal estimates and, second, priors based on the
parameterization (8) encode prior information that treatment effects are
highly non-homogeneous. To mitigate RIC they develop a prior that depends on
an estimate of the propensity score $\widehat{\pi}_{i}$ as a 1-dimensional
summary of the covariates, while to address non-homogeneity they
reparameterize the regression as

$\displaystyle
Y_{i}(a)=\mu(X_{i},\widehat{\pi}_{i})+a\,\tau(X_{i})+\epsilon_{i}$

where $\mu(x,\widehat{\pi})$ captures the prognostic effect of the control
variables $X_{i}$ and $\tau(x)$ is exactly the treatment effect. Independent
BART priors are then placed on $\mu(\cdot)$ and $\tau(\cdot)$, with the prior
on $\tau(\cdot)$ encoding our prior beliefs about the degree of treatment
effect heterogeneity.

Linero and Zhang, (2022) consider estimation of direct and indirect effects in
causal mediation using BART models. Linked with the concept of RIC, they also
show that naively specified priors can be highly _dogmatic_ (Linero,, 2021) in
the sense of encoding a prior belief that, on average, the mediator and
outcome potential outcomes are unconfounded (and hence that inference for the
average mediation effects can proceed as though there were no confounding
present). Prior dogmatism induces regularization-induced confounding by giving
a strong prior preference to encourage the model to attribute causal effects
on the outcome as being due to the treatment rather than the confounders. To
address this issue, they include “clever covariates”
$\widehat{m}_{ai}=\widehat{E}(M_{i}\mid A_{i}=a,X_{i})$ into the outcome model
for $a\in\\{0,1\\}$. These clever covariates are analagous to the propensity
score estimate $\widehat{\pi}_{i}$ used to correct for RIC in the BCF. Linero
and Zhang, (2022) then introduce the _Bayesian causal mediation forests_
(BCMF) model

$\displaystyle\begin{split}Y_{i}(a,m)&=\mu_{y}(m,a,X_{i})+\epsilon_{i}\\\
M_{i}(a)&=\mu_{m}(a,X_{i})+\epsilon_{i}\end{split}$ (9)

where the functions $\mu_{y}(\cdot,0,\cdot)$, $\mu_{y}(\cdot,1,\cdot)$,
$\mu_{m}(0,\cdot)$, and $\mu_{m}(1,\cdot)$ are given independent BART priors
and the clever covariates $\widehat{m}_{0i}$ and $\widehat{m}_{1i}$ are
included as predictors into the BART model for the outcome.

While the model (9) accomplishes the goal of estimating average mediation
effects well (i.e., it solves the problem of RIC), it does not appropriately
control the degree of heterogeneity in the conditional average mediation
effects. A contribution of this work is to use the insights behind the
parameterization of BCFs to develop a model that applies seperate
regularization to the direct and indirect effects.

##  3 BART for Heterogeneous Mediation Effects

We now introduce our causal mediation analysis model; a “Bayesian backfitting”
algorithm for fitting this model is given in the supplementary material.
Analogous to BCFs, the models presented enable direct regularization of
$\delta_{a}(x)$ and $\zeta_{a}(x)$. This type of direct regularization has
been shown to be crucial in generating dependable estimates of heterogneous
causal effects in other contexts (Hahn et al.,, 2020; Nie and Wager,, 2021).

For numeric outcomes and mediators, we specify the models

$\displaystyle Y_{i}(a,m)$
$\displaystyle=\mu(X_{i})+a\,\zeta(X_{i})+m\,d(X_{i})+\epsilon_{i},$ (10)
$\displaystyle M_{i}(a)$
$\displaystyle=\mu_{m}(X_{i})+a\,\tau_{m}(X_{i})+\nu_{i}$ (11)

where independent BART priors are specified for
$(\mu,\zeta,d,\mu_{m},\tau_{m})$. The mediator model (11) simply corresponds
to a BCF model as proposed by Hahn et al., (2020), with $\tau_{m}(x)$
corresponding to a heterogeneous causal effect of the treatment on the
outcome. The outcome model (10), on the other hand, corresponds to a _varying
coefficient BART_ (VC-BART) model as proposed by Deshpande et al., (2020),
with the treatment $A_{i}$ and mediator $M_{i}$ entering linearly.

The model (10)–(11) is a varying coefficient version of commonly used LSEMs,
with the coefficients modeled nonparametrically as a function of $X_{i}$.
Because of this, the conditional average mediation effects are also
expressible as products of coefficients as

$\displaystyle\zeta_{a}(x)$ $\displaystyle=\zeta(x),\qquad\qquad\text{and}$
$\displaystyle\delta_{a}(x)$ $\displaystyle=\tau_{m}(x)\,d(x).$

Hence, this parameterization allows us to isolate the components $\zeta(x)$
and $\delta(x)$ and apply differing amounts of regularization to them. Note
that (10)–(11) assumes that no interaction exists between the mediator and
treatment in the outcome model, and hence $\delta_{a}(x)$ and $\zeta_{a}(x)$
do not depend on the treatment level $a$, i.e., $\delta_{0}(x)=\delta_{1}(x)$
and $\zeta_{0}(x)=\zeta_{1}(x)$.

For average effects, note that the marginal distribution of
$Y_{i}\\{a,M_{i}(a^{\prime})\\}$ is given by

$\displaystyle f\big{(}Y_{i}\\{a,M_{i}(a^{\prime})\\}=y\big{)}=\int
f\big{(}Y_{i}\\{a,M_{i}(a^{\prime})\\}=y\mid X_{i}=x\big{)}\,f(X_{i}=x)\ dx.$
(12)

It is therefore necessary to specify a model for the distribution of the
covariates. Often, when this distribution is not modeled explicitly, the
empirical distribution is used instead as an estimate, i.e.
$F_{X}(dx)=\sum_{i}\omega_{i}\,\delta_{X_{i}}(dx)$ where $\delta_{x}(\cdot)$
denotes a point-mass distribution at $x$ and $\omega_{i}=n^{-1}$. An
alternative to the empirical distribution is the _Bayesian Bootstrap_ (BB
Rubin,, 1981), which respects our inherent uncertainty in $F_{X}$ while
cleanly avoiding the need to model the distribution of the covariates. The BB
is similar to the empirical distribution, but instead of setting
$\omega_{i}=n^{-1}$ we use an improper prior
$\pi(\omega)=\prod_{i}\omega_{i}^{-1}$; this leads to the posterior
distribution $\omega\sim\text{Dirichlet}(1,\ldots,1)$ for the weights. Under
the BB the average effects are identified as
$\bar{\delta}=\sum_{i}\omega_{i}\,\delta(X_{i})\qquad\text{and}\qquad\bar{\zeta}=\sum_{i}\omega_{i}\,\zeta(X_{i}).$

###  3.1 Controlling Heterogeneity Through Prior Specification

An important advantage of the models (10)–(11) relative to the model of Linero
and Zhang, (2022) is that we can shrink the model fits towards _homogeneous_
mediator and treatment effects through judicious choice of hyperparameters;
after doing this, we can be confident that any heterogeneity we _do_ detect is
well-supported by the data, rather than being the result of instability due to
the use of nonparametric estimators.

The degree of heterogeneity of the direct effect can be controlled via the
prior specification for $\zeta(x)$. Specifically, we can shrink $\zeta(x)$ to
a constant function, with few effect moderators, by (i) setting the parameter
$\alpha$ in (7) to a small value (say, $\alpha=0.5$) so that most trees do not
include covariates and (ii) using a smaller number of trees (say, $m=20$).
Using the same strategies, we can control the degree of heterogeneity in
$\tau_{m}(x)$, which represents the causal effect of the treatment on the
mediator.

The considerations for the indirect effects are slightly more complicated, as
$\delta(x)=d(x)\,\tau_{m}(x)$ consists of two components. Note that
heterogeneity in $\delta(x)$ is inevitable if $\tau_{m}(x)$ is non-constant.
However, if $\tau_{m}(x)$ is constant then $\delta(x)$ can be made homogeneous
by shrinking $d(x)$ towards a constant function. Accordingly, we adopt the
same strategy for $d(x)$ as we adopt for $\zeta(x)$ and $\tau_{m}(x)$: using a
small number of trees and setting $\alpha$ small.

###  3.2 Modeling Non-Numeric Data

Binary mediators can also be easily incorporated by using the nonparametric
probit regression model

$\displaystyle[M_{i}(a)\mid X_{i}=x]$
$\displaystyle\sim\operatorname{Bernoulli}[\Phi\\{\mu_{m}(x)+a\,\tau_{m}(x)\\}],$

with BART priors again used for $(\mu_{m},\tau_{m})$. Under sequential
ignorability and (10), we can identify the mediation effects as

$\displaystyle\zeta_{a}(x)=\zeta(x)\qquad\text{and}\qquad\delta_{a}(x)=d(x)\left[\Phi\\{\mu_{m}(x)+\tau_{m}(x)\\}-\Phi\\{\mu_{m}(x)\\}\right]$

where $\Phi(\cdot)$ is the cumulative distribution function of a standard
normal random variable. Similar expressions for the direct and indirect
effects can also be computed when $[Y_{i}(a,m)\mid
X_{i}=x]\sim\operatorname{Bernoulli}[\Phi\\{\mu(x)+a\,\zeta(x)+m\,d(x)\\}]$
and $[M_{i}(a)\mid X_{i}=x]\sim
N\\{\mu_{m}(x)+a\,\tau_{m}(x),\sigma^{2}_{m}\\}$ by noting that

$\displaystyle E[Y_{i}\\{a,M_{i}(a^{\prime})\\}\mid
X_{i}=x]=\Phi\left(\frac{\mu(x)+a\,\zeta(x)+\\{\mu_{m}(x)+a^{\prime}\,\tau_{m}(x)\\}\,d(x)}{\sqrt{1+d^{2}(x)\,\sigma_{m}^{2}}}\right),$

which can be derived by noting that the probit model implies that
$E[Y_{i}\\{a,M_{i}(a^{\prime})\\}\mid
X_{i}=x]=\Pr(\epsilon_{i}\leq\mu(x)+a\,\zeta(x)+M_{i}(a^{\prime})\,d(x)\mid
X_{i}=x)$ where $\epsilon_{i}\sim N(0,1)$ and $\epsilon_{i}$ is independent of
$M_{i}(a^{\prime})$. This implies, for example, that when $Y_{i}$ is binary
and $M_{i}$ is continuous we have

$\displaystyle\delta_{a}(x)=\Phi\left(\frac{\mu(x)+a\,\zeta(x)+d(x)\,\\{\mu_{m}(x)+\tau_{m}(x)\\}}{\sqrt{1+d^{2}(x)\,\sigma^{2}_{m}}}\right)-\Phi\left(\frac{\mu(x)+a\,\zeta(x)+d(x)\,\mu_{m}(x)}{\sqrt{1+d^{2}(x)\,\sigma^{2}_{m}}}\right).$

###  3.3 Corrections for Regularization Induced Confounding

As shown by Linero and Zhang, (2022), Bayesian nonparametric models for
mediation are also subject to the same RIC phenomenon as models for
observational data. They make the following recommendations to combat this:

  1. 1.

Add an estimate of the propensity score
$\widehat{\pi}_{i}=\widehat{\Pr}(A_{i}=1\mid X_{i}=x)$ to both the outcome
model and mediator model.

  2. 2.

Add an estimate of the mediator regression function
$\widehat{m}_{ai}=\widehat{E}(M_{i}\mid A_{i}=a,X_{i}=x)$ to the outcome
regression for $a\in\\{0,1\\}$.

See Linero and Zhang, (2022) for an extensive discussion of why it is
necessary to include these variables and a thorough simulation experiment. In
principle it does not matter how
$(\widehat{\pi}_{i},\widehat{m}_{0i},\widehat{m}_{1i})$ are obtained, aside
from the fact that $\widehat{\pi}_{i}$ should depend only on $(A_{i},X_{i})$
and $\widehat{m}_{ai}$ should depend only on $(M_{i},A_{i},X_{i})$; we use
BART to estimate these quantities.

###  3.4 Summarizing the Posterior

In addition to the (conditional) average mediation effects, it is also of
interest to produce interpretable summaries of the fit of the BCMF model to
the data. These summaries can help identify subpopulations that respond
differently to the treatment, help us interpret the impact of the effect
moderators, and provide insight into BCMF’s predictive process.

Woody et al., (2021) propose a general framework for posterior summarization
based on projecting complex models onto interpretable surrogate models. For
example, we might project the samples of $\delta(x)$ onto an _additive
function_ $\gamma(x)=\alpha+\sum_{j=1}^{p}\gamma_{j}(x_{j})$, the idea being
that if $\gamma(x)$ is a good approximation to $\delta(x)$ then we can use the
interpretable structure of $\gamma(x)$ to understand how $\delta(x)$ makes
predictions.

For simplicity we focus on the indirect effect $\delta(x)$ and consider two
classes of summaries:

  * •

An additive function, constructed as
$\widehat{\gamma}=\arg\min_{\gamma}\sum_{i}\\{\delta(X_{i})-\gamma(X_{i})\\}^{2}+q_{\lambda}(\gamma)$
where $\gamma(x)=\alpha+\sum_{j}\gamma_{j}(x_{j})$ with each
$\gamma_{j}(x_{j})$ being a univariate spline. Here,
$q_{\lambda}(\gamma)=\sum_{j}q_{\lambda}(\gamma_{j})$ is a roughness penalty
for the individual additive components. This summary can be computed by
fitting a _generalized additive model_ (GAM, see Wood,, 2006, for a review)
with $\\{\delta(X_{i}):i=1,\ldots,n\\}$ as the outcome and
$\\{X_{i}:i=1,\ldots,n\\}$ as the predictors.

  * •

A decision tree summary, where $\widehat{\gamma}$ is constructed by running
the CART algorithm (Breiman et al.,, 1984), treating
$\\{\delta(X_{i}):i=1,\ldots,n\\}$ as the outcome and
$\\{X_{i}:i=1,\ldots,n\\}$ as the predictors.

CART summaries are useful for identifying subpopulations with substantially
different treatment effects, while GAM summaries are useful for understanding
the impact of the different predictors on the estimated effects in isolation.
See Section 4 for an illustation of this approach on the MEPS dataset.

As an overall measure of the quality of the summaries we use the squared
correlation between $\delta(x)$ and $\gamma(x)$ given by

$\displaystyle
R^{2}=1-\frac{\sum_{i}\\{\delta(X_{i})-\widehat{\gamma}(X_{i})\\}^{2}}{\sum_{i}\\{\delta(X_{i})-\widehat{\delta}\\}^{2}}$

where $\widehat{\delta}=n^{-1}\sum_{i}\delta(X_{i})$; Woody et al., (2021)
refer to $R^{2}$ as the “summary $R^{2}$.”

##  4 Medical Expenditure Panel Survey Data

We now apply our model to a subset of the the Medical Expenditure Panel Survey
(MEPS). We focus on the questions of whether (i) there is a causal effect of
smoking on an individual’s expected annual medical expenditures, (ii) there is
evidence that the effect is entirely mediated by the effect of smoking on an
health, and (iii) whether any of the proposed confounders also act as
modifiers of the indirect effect. We take the outcome $Y_{i}$ to be the
logarithm of an individual’s annual net medical expenditure reported in the
2012 survey, the treatment $A_{i}$ to be whether an individual is a smoker
(yes or no), and the mediator $M_{i}$ to be an ordinal measure of overall
self-perceived health (1: excellent, 2: very good, 3: good, 4: fair, 5: poor).

At the outset, we note that a naive two-sample $t$-test for a difference in
medical expendtures between smokers and non-smokers shows that there is strong
evidence ($P$-value $<0.0005$) that non-smokers pay _more_ in medical
expenditures than smokers. Accordingly, it is important to control for
confounders in assessing any causal relationships. Our model includes the
following patient attributes as possible confounders:

  * •

age: Age in years.

  * •

bmi: Body mass index, which may act as a post-treatment confounder of health
and medical expenditures.

  * •

education_level: Education in years.

  * •

income: Total family income per year.

  * •

poverty_level: Family income as percentage of the poverty line.

  * •

region: Northeast, West, South, or Midwest.

  * •

sex: Male or female.

  * •

marital_status: Married, divorced, separated, or widowed.

  * •

race: White, Pacific Islander, Indigenous, Black, Asian, or multiple races.

  * •

seatbelt: whether an individual wears a seatbelt in a car (always, almost
always, sometimes, never, seldom, never drives/rides in a car).

Figure 2: Posterior distribution for the average direct effect $\bar{\zeta}$
(top), average indirect effect $\bar{\delta}$ (bottom), and average total
effect $\bar{\tau}=\bar{\zeta}+\bar{\delta}$ (middle).

The posterior distribution of the average direct and indirect effect is shown
in Figure 2. We see that, under the sequential ignorability assumption, there
is evidence of both a direct and indirect effect of smoking on expenditures.
Interestingly, these effects are in _opposite directions_ and cancel each
other out to a large extent. As a result, the sign of the total effect is
uncertain. This illustrates an important potential benefit of a mediation
analysis: we can establish a causal relationship between smoking and medical
expenditures that we could not if we restricted attention strictly to the
total effect.

Figure 3: Posterior summarization of the indirect effect using a single
regression tree. Figure 4: Posterior density for the average indirect effect
within subgroups $G_{1},G_{2},...,G_{5}$ from the terminal nodes in Figure 3.

###  4.1 Posterior Summarization

We use the summarization strategies outlined in Section 3.4 to interpret the
model fit and better understand the covariates and interactions contributing
to the heterogeneity in the indirect effect; specifically, we project the
indirect effect function $\delta(x)$ onto a single regression tree and an
additive function.

We first consider a CART summary of the posterior mean of of $\delta(x)$,
which was obtained on a preliminary model fit to the MEPS dataset. According
to the regression tree summary in Figure 3, race, age, and sex are the most
significant effect modifiers for the indirect effect. Motivated by the
subgroups found in Figure 3, in Figure 4 we display the average indirect
effects within various subgroups formed by age and race. We find that the
largest indirect effects occur for white middle-aged individuals, while the
smallest effects are for non-white young adults.

Figure 5 and Figure 6 display the results obtained from the GAM summary for
continuous and discrete variables, respectively. These figures again highlight
the importance of age and race as effect modifiers, indicating that older and
white individuals have higher indirect effects on medical expenditures
mediated by perceived health status.

Figure 5: Posterior summarization of the indirect effect using a GAM for the
continuous variables. The projection of the posterior mean is given by the
dashed line while the shaded area gives a posterior 95% credible band of the
projection. Figure 6: Posterior distributions of the impact of categorical
variables on the indirect effect using the GAM projection. Solid lines give
the posterior mean of the associated coefficient in the GAM model, thick bars
are 66% credible bands, thin bars are 95% credible bands.

In summary, both the CART and GAM summaries reveal that age and race are
important effect modifiers. To measure the adequacy of the summary function
approximations, Figure 7 presents both the posterior distribution of $R^{2}$
obtained from fitting the summaries to each posterior sample of
$\delta(\cdot)$ and a single $R^{2}$ obtained from fitting the summaries to
the posterior mean of $\delta(\cdot)$. Our analysis shows that regression tree
is slightly better as a summary than a GAM, suggesting that the interactions
detected in Figure 3 provide important insight into the model’s predictive
process for $\delta(\cdot)$.

Figure 7: Posterior distribution of the summary $R^{2}$ for the regression
tree and GAM summaries of $\delta(x)$. The black line indicates the summary
$R^{2}$ for the posterior mean of $\delta(x)$.

###  4.2 Comparison of BCMF and a LSEM

To evaluate the practical usefulness of the BCMF model (10)–(11), we compare
its predictive performance to that of an LSEM with interactions between the
treatment, mediator, and covariates,

$\displaystyle\begin{split}Y_{i}(a,m)&=\beta_{0Y}+X_{i}^{T}\beta_{Y}+a(\gamma_{0Y}+X_{i}^{T}\gamma_{Y})+m(\xi_{0}+X_{i}^{T}\xi)+\epsilon_{i}\\\
M_{i}(a)&=\beta_{0M}+X_{i}^{T}\beta_{M}+a(\gamma_{0M}+X_{i}^{T}\gamma_{M})+\nu_{i}.\end{split}$
(13)

This model allows for heterogeneous mediation effects, but restricts them to
linear functions of the confounders. To quantify the uncertainty of the LSEM
estimates, we use the residual bootstrap. By comparing the predictive
performance of these two models, we can assess whether the added complexity of
the BCMF model is warranted.

Figure 8: Top: Boxplot displaying the distribution of the estimated individual
direct and indirect effects for the BCMF and LSEM models fit to the original
MEPS dataset. Bottom: Boxplot displaying the distribution of the estimated
individual effects for the BCMF and LSEM, stratified by race.

To understand the salient differences between the predictions made from the
BCMF and LSEM, we compare the estimates of $\delta(X_{i})$ and $\zeta(X_{i})$
for each individual, both as a whole and stratified by race, in Figure 8.
Figure 8 presents the estimates of $\delta(X_{i})$ and $\zeta(X_{i})$ for each
individual, both in aggregate and stratified by race. The effect estimates of
the BCMF are substantially less heterogeneous than the LSEM, with the LSEM
estimating a substantial number of both positive and negative effects for both
$\delta(X_{i})$ and $\zeta(X_{i})$. Additionally, we see substantially less
heterogeneity across race; for example, the LSEM makes counterintuitive
predictions about both the direct and indirect effect of smoking within the
group of Pacific Islanders. While the MEPS dataset is large, there are
relatively few Pacific Islanders in the data, and in the subset we analyzed
only 13 of them smoke. By applying regularization, the BCMF shrinks the direct
and indirect effects within this subpopulation closer to those of the other
races.

Variable |  $R_{\text{test}}$ (BART) |  $R_{\text{test}}$ (LSEM) |  $P$-value
---|---|---|---
phealth | 0.445 | 0.419 | $0.0001$
log(Y) | 0.454 | 0.431 | $0.0005$
Table 1: Held-out correlation for the mediator (phealth) and outcome (logY)
across all individuals for the BART and LSEM fits, and the $p$-value for a
paired Wilcoxon matched pairs signed-rank test comparing the predictive
performance on held-out data for the two models.

Next, we fit the BCMF and LSEM models to the same training set of $n=8056$
individuals and compute predictions
$(\widehat{Y}_{i,\textnormal{lsem}},\widehat{M}_{i,\textnormal{lsem}},\widehat{M}_{i,\textnormal{bcmf}},\widehat{Y}_{i,\textnormal{bcmf}})$
on the test set of $n=8057$ individuals using the fitted models. We use these
predictions on the test set to evaluate the performance of the model in three
ways. First, we consider the correlation between $(M_{i},Y_{i})$ and their
predictions on the test set. Second, we perform a paired Wilcoxon signed-rank
test comparing the squared difference
$(Y_{i}-\widehat{Y}_{i,\textnormal{lsem}})^{2}$ to
$(Y_{i}-\widehat{Y}_{i,\textnormal{bcmf}})^{2}$ (and similarly for $M_{i}$).
Results are given in Table 1, and we see both that the correlation is somewhat
higher for the BCMF model than the LSEM model, and that the difference in
performance was highly statistically significant according to the signed-rank
test.

Term | Estimate | Standard Error | Statistic |  $P$-value
---|---|---|---|---
$\widehat{M}_{\text{lsem}}$ | 0.1561 | 0.0803 | 1.9447 | 0.0518
$\widehat{M}_{\text{bcmf}}$ | 0.8477 | 0.0809 | 10.4835 | $<0.0001$
$\widehat{Y}_{\text{lsem}}$ | 0.1509 | 0.0725 | 2.0824 | 0.0373
$\widehat{Y}_{\text{bcmf}}$ | 0.9030 | 0.0747 | 12.0840 | $<0.0001$
Table 2: Coefficient estimates for the linear model that aggregates the linear
and BART fits on the MEPS test data.

Our third comparison considers _stacking_ (Wolpert,, 1992) the predictions of
the BCMF and LSEM by fitting the linear models
$Y_{i}=\beta_{0}+\beta_{1}\,\widehat{Y}_{i,\textnormal{lsem}}+\beta_{2}\,\widehat{Y}_{i,\textnormal{bcmf}}+\epsilon_{i}$
(and similarly for $M_{i}$). Results of the stacking procedure are given in
Table 2. From this fit, we see that the linear model relies much more heavily
on the predictions from the BCMF than the linear model, and that the BCMF
predictions are much more statistically significant than the predictions from
the LSEM (in the sense that there is strong evidence that the BCMF predictor
improves upon the LSEM predictor, while there is only weak evidence of the
converse). Interestingly, the LSEM predictions _are_ found to be statistically
significant, suggesting that a modification of the BCMF that also includes
_linear_ adjustments for the confounders (i.e., includes linear terms
$x^{\top}b$ in the functions $(\mu(x),\zeta(x),d(x),\mu_{m}(x),\tau_{m}(x))$)
may improve the fit of the model.

###  4.3 Simulation Study

We now conduct a simulation study to better understand the operating
characteristics of the BCMF model. Our study aims to answer the following
question: (i) Does the BCMF model perform better in terms of predictive
accuracy in estimating the mediation effects? (ii) Do the credible intervals
for $\delta(X_{i})$ and $\zeta(X_{i})$ attain coverage rates close to their
nominal levels? (iii) Can the BCMF model estimate the effects accurately
within the subgroups of the data identified by the CART summary?

#### Data Generating Mechanism

We use a data generating mechanism in which the confounders and treatment
assignment are sampled direct from the MEPS dataset, while the mediator and
outcome ground truths are obtained by fitting both our model and an LSEM to
the data. To assess the performance of both methods, we replicated each
simulation setting 200 times, with 8056 observations in the training set and
8057 in the testing set. We used the same training/testing split across all
simulated datasets to evaluate the coverage probability of the
confidence/credible intervals generated by each method.

A crucial difference between the LSEM and BCMF models is that the LSEM does
not regularize the mediation effects. Consequently, the LSEM produces a ground
truth for the mediation effects that is more heterogeneous than is expected in
practice, especially for subgroups of the population with a small sample size.
For example, since the MEPS dataset includes few Pacific Islanders, the LSEM’s
estimate of the effect of race as an effect modifier is unstable for this
group. To account for this, we also consider a third ground truth that is also
an LSEM but with the parameters of (13) instead estimated using the R-Learner
approach of Nie and Wager, (2021), which uses the lasso to reduce the amount
of heterogeneity.

#### Results: Individual Effects

We fit our model and the LSEM to each simulated dataset and measure point
estimates of the effects, the limits and width of 95% credible intervals, and
whether or not the interval captures the true parameter for each replication.
Using the 200 replications, we then compute the root mean square error,
absolute bias, average width of the intervals, and the coverage probability.

Figure 9: Individual simulation reslts for $\delta(x)$ under all combinations of fitting the BART/LSEM model under the BART/LSEM/RLEARN ground truths. Top left gives the coverage probability of nominal 95% credible intervals among all individuals, top right gives the root mean squared error, bottom left gives the absolute bias, and bottom right gives the average interval length. Figure 10: Individual simulation reslts for $\zeta(x)$ under all combinations of fitting the BART/LSEM model under the BART/LSEM/RLEARN ground truths. Top left gives the coverage probability of nominal 95% credible intervals among all individuals, top right gives the root mean squared error, bottom left gives the absolute bias, and bottom right gives the average interval length. Setting | Method | Coverage | RMSE | Bias | Length
---|---|---|---|---|---
BART | BART | 0.91 | 0.12 | 0.06 | 0.50
BART | LSEM | 0.93 | 0.23 | 0.07 | 0.83
LSEM | BART | 0.94 | 0.19 | 0.10 | 0.77
LSEM | LSEM | 0.94 | 0.22 | 0.01 | 0.84
RLEARN | BART | 0.88 | 0.13 | 0.07 | 0.50
RLEARN | LSEM | 0.94 | 0.21 | 0.03 | 0.80
Table 3: Coverage probability, root mean square error, and absolute bias for $\zeta(X_{i})$ across all individuals in the test set. Setting | Method | Coverage | RMSE | Bias | Interval Length
---|---|---|---|---|---
BART | BART | 0.95 | 0.03 | 0.02 | 0.14
BART | LSEM | 0.89 | 0.06 | 0.03 | 0.19
LSEM | BART | 0.96 | 0.05 | 0.02 | 0.21
LSEM | LSEM | 0.95 | 0.05 | 0.00 | 0.22
RLEARN | BART | 0.99 | 0.03 | 0.01 | 0.16
RLEARN | LSEM | 0.95 | 0.05 | 0.01 | 0.21
Table 4: Coverage probability, root mean square error, and absolute bias for
$\delta(X_{i})$ across all individuals in the test set.

We present the results of our simulation study in Figure 9 and Figure 10,
where we compare the BCMF and LSEM models under different combinations of
ground truth and fitted models. Table 3 and Table 4 summarize the results from
Figure 9 and Figure 10, respectively, across all individuals in the test set.
When the BCMF is fitted to the BCMF ground truth, it outperforms the LSEM in
terms of achieving close to the nominal coverage on interval estimates with
substantially lower interval lengths, root mean squared error, and absolute
bias. Interestingly, we also observed that the BCMF model is competitive in
terms of root mean squared error when the LSEM is used to generate the data.
We conjecture that this is due to the fact that the data generating mechanism
estimated by the LSEM fit to the original data, while still quite
heterogeneous, is homogeneous enough (and the effects are small enough) that
the benefits of the regularization of the BCMF outweigh the fact that the LSEM
is correctly specified. We observe similar behavior, which is even more
pronounced, when the R-Learner is used to generate the ground truth.

The results for the coverage of the LSEM and BCMF models reveal some
interesting results. Surprisingly, the LSEM appears to be robust in terms of
coverage, although it produces much larger intervals compared to the BCMF.
However, for the individual-level direct effects, the BCMF performed poorly
for some individuals in achieving nominal coverage. We further investigate
this behavior in the supplementary material and find that the BCMF did not
attain nominal coverage for individuals with highly heterogeneous effects,
i.e., those whose conditional average mediation effects
$(\delta(X_{i}),\zeta(X_{i}))$ deviate greatly from the average effects
$(\bar{\delta},\bar{\zeta})$. This behavior is expected, as the BCMF model is
explicitly designed to shrink towards a structure with a small degree of
heterogeneity. While this reduces the power to detect strongly heterogeneous
effects, it does not inflate the Type I error in detecting heterogeneity,
making the BCMF model conservative in detecting heterogeneity.

#### Results: Average and Subgroup Average Effects

The BCMF also produces reliable estimates of the average mediation effects
within subpopulations. We consider here both fixed and data-dependent
subgroups obtained under the BCMF ground truth. The fixed subgroups are the
groups identified by the terminal nodes in Figure 3: age $\geq$ 67, non-white
and 34 $\leq$ age $<$ 67, non-white and age $<$ 34, white and 34 $\leq$ age
$<$ 67, and white and age $<$ 34\. The data-determined subgroups are
determined through posterior projection summarization, by fitting a tree and
identifying the terminal node groups for each simulated dataset. A comparison
of the inferences for the average effects under each simulation scenario is
given in the supplementary material.

Group | Indirect Effect | Direct Effect
---|---|---
$\texttt{age}\geq 67$ | 0.99 | 0.92
$\texttt{non-white},34\leq\texttt{age}<67$ | 0.96 | 0.88
$\texttt{non-white},\texttt{age}<34$ | 0.86 | 0.94
$\texttt{white},34\leq\texttt{age}<67$ | 0.88 | 0.92
$\texttt{white},\texttt{age}<34$ | 0.96 | 0.96
Average | 0.93 | 0.93
Dynamic | 0.95 | 0.95
Table 5: Subgroup coverage probability of $\zeta_{a}(A_{i})$ and
$\delta_{a}(A_{i})$ using the subgroups in Figure 3.

Table 5 shows the results of the simulation for both the fixed subgroups and
for the data-dependent subgroups (labeled “Dynamic”). We see that the BCMF
produces intervals whose coverage is close to the nominal level, with slightly
poorer results in the non-white groups. Interestingly, the coverage for data-
dependent groups have _higher_ coverage for the credible intervals, and in
fact attain exact 95% coverage for both the direct and indirect effects. The
intervals for the average effects $\bar{\delta}$ and $\bar{\zeta}$ also attain
close to nominal coverage.

##  5 Discussion

In this paper we introduced a Bayesian causal mediation forest (BCMF) model
that can separately identify and regularize the conditional average natural
direct and indirect effects using varying coefficient models. Our approach is
reminiscent of LSEMs, making it easy to identify these effects as products of
varying coefficients. Additionally, we demonstrate that our model produces
lower prediction error than a comparable LSEM on both real and simulated MEPS
data. Furthermore, we argue that our model is conservative in estimating
heterogeneity since it assumes small and mostly homogeneous mediation effects.
We also provide posterior summarization methods for interpreting model fit and
subgroup detection.

To improve our methods and analysis, there are several directions one could
take. First, we can improve the models for the outcome and mediator. For
instance, log medical expenditure exhibits heteroskedasticity, with the
variance of $Y_{i}$ and $X_{i}$ having a complex relationship, as demonstrated
by Linero et al., (2020). Additionally, since the mediator in this problem is
ordinal, and empirically is well-approximated with a rounded normal
distribution; thus, we can improve our model by using a cumulative probit
model for $M_{i}$ rather than a normal model. The impact of using a continuous
model for $M_{i}$ rather than an ordinal model is unclear, and warrants
further investigation.

Exclusion of individuals with no medical expenditure from the analysis (which
we have done here) is problematic, as the likelihood of incurring medical
expenditure is likely to be linked with smoking status. As a further
improvement to our analysis, a better approach would be to use principal
stratification (Frangakis and Rubin,, 2002). This approach would estimate the
causal effect of smoking on medical expenditures within the strata of
individuals who incur medical expenditures, irrespective of their smoking
status. In such an analysis, it is assumed that all individuals who would
incur medical expenses if they did not smoke would also incur medical expenses
if they did smoke. This would enable a more honest evaluation of the causal
effect of smoking on medical expenditures.

Lastly, while our model performs well in terms of root mean squared error, for
some individuals it does not quite reach the nominal coverage level for
credible intervals. In the supplementary material, we show that our BCMF
under-covers for individuals whose conditional mediation effect differs
significantly from the average effects $\bar{\delta}$ and $\bar{\zeta}$.
Whether this is a problem that can be fixed or simply a consequence of using a
model that shrinks towards homogeneous effects warrants further investigation.
Code reproducing our analysis and simulation results is available at
www.github.com/vcbcmf/vcbcmf.

## References

  * Albert, (2008) Albert, J. M. (2008).  Mediation analysis via potential outcomes models.  Statistics in Medicine, 27(8):1282–1304.
  * Athey, (2017) Athey, S. (2017).  Beyond prediction: Using big data for policy problems.  Science, 355(6324):483–485.
  * Baron and Kenny, (1986) Baron, R. M. and Kenny, D. A. (1986).  The moderator–mediator variable distinction in social psychological research: Conceptual, strategic, and statistical considerations.  Journal of Personality and Social Psychology, 51(6):1173.
  * Breiman et al., (1984) Breiman, L., Friedman, J. H., Olshen, R. A., and Stone, C. J. (1984).  Classification and Regression Trees.  Chapman and Hall, New York, 1st edition.
  * Chipman et al., (2010) Chipman, H. A., George, E. I., and McCulloch, R. E. (2010).  Bart: Bayesian additive regression trees.  The Annals of Applied Statistics, 4(1):266–298.
  * Deshpande et al., (2020) Deshpande, S. K., Bai, R., Balocchi, C., Starling, J. E., and Weiss, J. (2020).  VCBART: Bayesian trees for varying coefficients.
  * Dorie et al., (2019) Dorie, V., Hill, J., Shalit, U., Scott, M., and Cervone, D. (2019).  Automated versus do-it-yourself methods for causal inference: Lessons learned from a data analysis competition.  Statistical Science, 34(1):43–68.
  * Farbmacher et al., (2022) Farbmacher, H., Huber, M., Lafférs, L., Langen, H., and Spindler, M. (2022).  Causal mediation analysis with double machine learning.  The Econometrics Journal, 25(2):277–300.
  * Frangakis and Rubin, (2002) Frangakis, C. E. and Rubin, D. B. (2002).  Principal stratification in causal inference.  Biometrics, 58(1):21–29.
  * Hahn et al., (2020) Hahn, P. R., Murray, J. S., and Carvalho, C. M. (2020).  Bayesian regression tree models for causal inference: Regularization, confounding, and heterogeneous effects (with discussion).  Bayesian Analysis, 15(3):965–1056.
  * Hill et al., (2020) Hill, J., Linero, A., and Murray, J. (2020).  Bayesian additive regression trees: A review and look forward.  Annual Review of Statistics and Its Application, 7(1):251–278.
  * Hill, (2011) Hill, J. L. (2011).  Bayesian nonparametric modeling for causal inference.  Journal of Computational and Graphical Statistics, 20(1):217–240.
  * Imai et al., (2010) Imai, K., Keele, L., and Tingley, D. (2010).  A general approach to causal mediation analysis.  Psychological Methods, 15(4):309.
  * Kershaw et al., (2010) Kershaw, K. N., Mezuk, B., Abdou, C. M., Rafferty, J. A., and Jackson, J. S. (2010).  Socioeconomic position, health behaviors, and c-reactive protein: a moderated-mediation analysis.  Health Psychology, 29(3):307.
  * Kim et al., (2017) Kim, C., Daniels, M. J., Marcus, B. H., and Roy, J. A. (2017).  A framework for Bayesian nonparametric inference for causal effects of mediation.  Biometrics, 73(2):401–409.
  * Künzel et al., (2019) Künzel, S. R., Sekhon, J. S., Bickel, P. J., and Yu, B. (2019).  Metalearners for estimating heterogeneous treatment effects using machine learning.  Proceedings of the national academy of sciences, 116(10):4156–4165.
  * Linero, (2021) Linero, A. R. (2021).  In nonparametric and high-dimensional models, Bayesian ignorability is an informative prior.
  * Linero et al., (2020) Linero, A. R., Sinha, D., and Lipsitz, S. R. (2020).  Semiparametric mixed-scale models using shared Bayesian forests.  Biometrics, 76(1):131–144.
  * Linero and Zhang, (2022) Linero, A. R. and Zhang, Q. (2022).  Mediation analysis using bayesian tree ensembles.  Psychological Methods.
  * MacKinnon, (2008) MacKinnon, D. P. (2008).  Introduction to statistical mediation analysis.
  * MacKinnon and Dwyer, (1993) MacKinnon, D. P. and Dwyer, J. H. (1993).  Estimating mediated effects in prevention studies.  Evaluation Review, 17(2):144–158.
  * Muller et al., (2005) Muller, D., Judd, C. M., and Yzerbyt, V. Y. (2005).  When moderation is mediated and mediation is moderated.  Journal of Personality and Social Psychology, 89(6):852.
  * Nie and Wager, (2021) Nie, X. and Wager, S. (2021).  Quasi-oracle estimation of heterogeneous treatment effects.  Biometrika, 108(2):299–319.
  * Obermeyer and Emanuel, (2016) Obermeyer, Z. and Emanuel, E. J. (2016).  Predicting the future—big data, machine learning, and clinical medicine.  The New England journal of medicine, 375(13):1216.
  * Pearl, (2001) Pearl, J. (2001).  Direct and indirect effects.  In Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence, page 411–420. Morgan Kaufmann Publishers Inc.
  * Preacher et al., (2007) Preacher, K. J., Rucker, D. D., and Hayes, A. F. (2007).  Addressing moderated mediation hypotheses: Theory, methods, and prescriptions.  Multivariate Behavioral Research, 42(1):185–227.
  * Robins and Greenland, (1992) Robins, J. M. and Greenland, S. (1992).  Identifiability and exchangeability for direct and indirect effects.  Epidemiology, pages 143–155.
  * Rubin, (1974) Rubin, D. B. (1974).  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of Educational Psychology, 66(5):688–701.
  * Rubin, (1981) Rubin, D. B. (1981).  The Bayesian bootstrap.  The Annals of Statistics, pages 130–134.
  * Rubin, (2004) Rubin, D. B. (2004).  Direct and indirect causal effects via potential outcomes.  Scandinavian Journal of Statistics, 31(2):161–170.
  * Tchetgen and Shpitser, (2012) Tchetgen, E. J. T. and Shpitser, I. (2012).  Semiparametric theory for causal mediation analysis: Efficiency bounds, multiple robustness, and sensitivity analysis.  The Annals of Statistics, 40(3):1816.
  * VanderWeele, (2016) VanderWeele, T. J. (2016).  Mediation analysis: A practitioner’s guide.  Annual Review of Public Health, 37:17–32.
  * Wendling et al., (2018) Wendling, T., Jung, K., Callahan, A., Schuler, A., Shah, N. H., and Gallego, B. (2018).  Comparing methods for estimation of heterogeneous treatment effects using observational data from health care databases.  Statistics in Medicine, 37(23):3309–3324.
  * Wolpert, (1992) Wolpert, D. H. (1992).  Stacked generalization.  Neural Networks, 5(2):241–259.
  * Wood, (2006) Wood, S. N. (2006).  Generalized Additive Models: An Introduction with R.  Chapman and Hall/CRC, New York, 1st edition.
  * Woody et al., (2021) Woody, S., Carvalho, C. M., and Murray, J. S. (2021).  Model interpretation through lower-dimensional posterior summarization.  Journal of Computational and Graphical Statistics, 30(1):144–161.
  * Yeager et al., (2019) Yeager, D. S., Hanselman, P., Walton, G. M., Murray, J. S., Crosnoe, R., Muller, C., Tipton, E., Schneider, B., Hulleman, C. S., Hinojosa, C. P., Paunesku, D., Romero, C., Flint, K., Roberts, A., Trott, J., Iachan, R., Buontempo, J., Yang, S. M., Carvalho, C. M., Hahn, P. R., Gopalan, M., Mhatre, P., Ferguson, R., Duckworth, A. L., and Dweck, C. S. (2019).  A national experiment reveals where a growth mindset improves achievement.  Nature, 573(7774):364–369.
  * Zheng and van der Laan, (2012) Zheng, W. and van der Laan, M. J. (2012).  Targeted maximum likelihood estimation of natural direct effects.  The International Journal of Biostatistics, 8(1):1–40.