##  Appendix N Microbenchmarks of Hercules and POSEIDON

TABLE VIII: Microbenchmarks of Hercules and POSEIDON Functionality | Execution time (s) | Comm. (MB)
---|---|---
Approx_Sigmoid |  0.017 — 0.017 | –
Appro_ReLu |  0.01 — 0.09 | –
Appro_Softmax |  0.07 — 0.07 | –
Average-pooling |  0.034—0.034 | –
Max-pooling |  5.73—1.91 |  23.5—6.9
DBootstrap |  0.09 — 0.09 |  4.5 — 1.5
FC layer |  0.097 — 0.0015 | –
FC layer-backprop |  0.14 — 0.002 | –
CV layer |  0.04 — 0.0006 | –
CV layer-backprop |  0.06 — 0.0009 | –
DKeySwitch |  0.06 — 0.06 |  22.14 — 22.14

We present the microbenchmark cost of Hercules and POSEIDON in performing
different functions. As shown in Table VIII, the bold font represents the
overhead of Hercules. For experimental configurations, We consider the number
of users $N=10$, the dimension of each sample is $h=32$, the number of neurons
in each layer is 64 or with kernel size $3\times 3$, and the dimension of
cyclotomic polynomial ring used in CKKS is set $\mathcal{N}=2^{13}$. All the
costs represent the processing of 1 sample per party.

We record the computation and communication costs (in an amortized way) of
calculating FC, CV, FC backpropagation, CV backpropagation, and different
activation functions in Hercules and POSEIDON. For the activation functions
sigmoid and sofmax, Hercules and POSEIDON exhibit the same overhead because
they both use the same function approximation method. For ReLU, since Hercules
uses a 20 bit-precision composite polynomial to fit the original function,
while POSEIDON simply uses a polynomial with a degree of 3 to approximate the
original ReLU (using the least squares method), it makes the cost of Hercules
larger than POSEIDON but derives a higher approximation accuracy. Note that
given the same error bound, the time required to calculate the polynomial
generated by the least square method under ciphertext is much longer than
ours. For max pooling, POSEIDON replaces it with Chebyshev function
interpolation, thereby obtaining a polynomial of degree 31 with a precision of
7 bits. Experimental results show that our 7-bit precision composite
polynomial is much less expensive than POSEIDON. This is mainly due to the
characteristics of composite polynomials. As discussed before, a composite
function with degree of $deg(G)$ can be calculated with a complexity of
$O(\log(deg(G)))$, while the computation complexity of calculating any
function $G$ is at least $\Theta(\sqrt{deg(G)})$.

The operations in FC, CV and their backpropagation are mainly composed of
matrix multiplication. As shown in Table VIII, Hercules is far superior to
POSEIDON in terms of calculating homomorphic matrix multiplication.
Specifically, POSEIDON adopts AP to achieve fast SIMD calculations. For the
multiplication of two $h\times h$-dimensional matrices, the complexity of the
homomorphic rotation operation required by AP is
$\max_{i\in[\mathbb{L}]}(\omega_{i}\times\log(h\times\omega_{i}))$. For
Hercules, as shown in Table I, the complexity required for the matrix
multiplication is only $3h+5\sqrt{h}$. Moreover, AP requires multiple copies
and zero padding operations for each row or column of the input matrix,
(depending on the number of neurons in each hidden layer, and the absolute
value of the difference between the row or column dimension of the matrix and
the corresponding hidden layer neuron. On the contrary, our method does not
require additional element copying except for a small amount of zero padding
in the initial stage to facilitate calculations. Therefore, Hercules obviously
exhibits lower communication overhead.

##  Appendix O Storage Overhead

TABLE IX: Comparison of storage cost Datasets | Storage cost (GB)
---|---
Hercules | POSEIDON
BCW | $0.0012$ | 0.095
ESR | $0.184$ | $11.778$
CREDIT | $0.061$ | $3.948$
MNIST | $4.906$ | $314.025$
SVHN | $54.9$ | $3515.6$
CIFAR-10 | $16.47$ | $2109.37$
CIFAR-100 | $16.47$ | $2109.37$
|  |

As shown in TABLE IX, we also count the storage overhead required for training
each scheme on different datasets, where the storage overhead is mainly
dominated by the size of the ciphertext that each user needs to save locally.
We observe that the storage overhead required by Hercules on each dataset is
much lower than that of POSEIDON. Since POSEIDON designs a alternating packing
(AP) method to pack ciphertext in the encryption process, this leads to
multiple copies and zero padding operations for each ciphertext, depending on
the number of neurons in each hidden layer. On the contrary, our method does
not require additional element copying except for a small amount of zero
padding in the initial stage to facilitate calculations.

##  Appendix P Comparison with Other Prior Works

Beyond the HE-based schemes, a wealth of MPC solutions based on interaction
between multiple servers have been proposed to design privacy-preserving
neural network training and prediction frameworks. These solutions may be not
practical in the real-world setting without multiple servers. Nevertheless, we
still make a rough comparison between these works and Hercules to demonstrate
the merits and demerits of each scheme. As shown in Table VII in the Appendix,
MPC protocol with multiple servers relies on splitting the training task into
two or more servers, where the servers are usually assumed to be non-
colluding. Then, the state-of-the-art secret sharing methods, including
arithmetic sharing [19], boolean sharing [8], and Yao’s garbled circuit [21]
are carefully integrated to efficiently implement various mathematical
operations under the ciphertext. It is computationally cost-effective,
avoiding the high communication overhead among large-scale users. For example,
given the MNIST dataset, we roughly compare the computational cost of
Hercules, ABY3 [8], and XONN [58] in training the same network architecture.
Here we use a three-layer architecture with 128 neurons in each layer. The
number of users in Hercules is set to 3 and the global training epoch is set
to 15. Experimental results show that Hercules requires a total of 24.3 hours
to train such a model, while ABY3 and XONN need 1.02 and 0.58 hours,
respectively. We remind that Hercules is operated in a different scenario (FL)
and threat model from those MPC-based schemes. It supports more participants
with collusion, and the cost only grows linearly with the number of users. In
contrast, those MPC-based solutions require to outsource training tasks among
limited computing servers, which may be impractical in some scenarios.


# Borrowing from Supplemental Sources to Estimate Causal Effects from a
Primary Data Source

Jeffrey A. Boatman
Division of Biostatistics
University of Minnesota
Minneapolis, MN 55455
<EMAIL_ADDRESS>
&David M. Vock
Division of Biostatistics
University of Minnesota
Minneapolis, MN 55455
<EMAIL_ADDRESS>
&Joseph S. Koopmeiners
Division of Biostatistics
University of Minnesota
Minneapolis, MN 55455
<EMAIL_ADDRESS>

###### Abstract

The increasing multiplicity of data sources offers exciting possibilities in
estimating the effects of a treatment, intervention, or exposure, particularly
if observational and experimental sources could be used simultaneously.
Borrowing between sources can potentially result in more efficient estimators,
but it must be done in a principled manner to mitigate increased bias and Type
I error. Furthermore, when the effect of treatment is confounded, as in
observational sources or in clinical trials with noncompliance, causal effect
estimators are needed to simultaneously adjust for confounding and to estimate
effects across data sources. We consider the problem of estimating causal
effects from a primary source and borrowing from any number of supplemental
sources. We propose using regression-based estimators that borrow based on
assuming exchangeability of the regression coefficients and parameters between
data sources. Borrowing is accomplished with multisource exchageability models
and Bayesian model averaging. We show via simulation that a Bayesian linear
model and Bayesian additive regression trees both have desirable properties
and borrow under appropriate circumstances. We apply the estimators to
recently completed trials of very low nicotine content cigarettes
investigating their impact on smoking behavior.

_K_ eywords Bayesian Additive Regression Trees $\cdot$ Bayesian Linear Model
$\cdot$ Bayesian Model Averaging $\cdot$ Borrowing $\cdot$ Causal Inference
$\cdot$ Multisource Exchangeability Models.

##  1 Introduction

The Center for Evaluation of Nicotine in Cigarettes (CENIC) conducts
randomized controlled trials (RCTs) investigating how marked decreases in the
nicotine content of cigarettes affect the use and effects of tobacco in
current smokers. Their work is part of a body of research meant to inform the
United States Food and Drug Administration (FDA) as it contemplates strategies
to reduce the public health burden of smoking. Not surprisingly, because
smokers enrolled in trials of very low nicotine content (VLNC) cigarettes
still have access to normal commercial cigarettes, noncompliance to randomized
treatment in these trials is very common (Nardone _and others_ , 2016). The
intention to treat (ITT) analysis reported for these trials thus estimates the
effect of nicotine reduction if smokers still have access to commercial
cigarettes. But ITT does not estimate the effect of a potential FDA-mandated
reduction in the nicotine content of cigarettes: in this hypothetical future,
normal nicotine cigarettes would not be legally available, and all smokers
would be forced to comply with a nicotine reduction. Causal effect estimators
addressing noncompliance must therefore be used in conjunction with ITT
estimators to fully understand the possible effects of FDA regulations on the
nicotine content of cigarettes. Unfortunately, estimating causal effects can
be inefficient if noncompliance is high or in the presence of strong
confounding.

One way to improve efficiency is to combine data sources for causal inference.
This is a particularly attractive option for CENIC given similarities in
entrance criteria and interventions across studies. CENIC project 1, study 2
(CENIC-P1S2) was a 6-week, randomized, $2\times 2$ factorial design to
evaluate the effect of VLNC versus normal nicotine content (NNC) cigarettes in
the presence or absence of transdermal nicotine (TDN)(Smith _and others_ ,
2019). Each treatment group is small ($n\approx 60$/group), posing further
challenges for estimating causal effects, but extensive supplementary data are
available from CENIC project 1, study 1 (CENIC-P1S1) (Donny _and others_ ,
2015) and from CENIC project 2 (CENIC-P2)(Hatsukami _and others_ , 2018).
CENIC-P1 was a 6-week RCT that randomized smokers to one of seven groups
consisting of usual brand cigarettes or experimental cigarettes with nicotine
content ranging from NNC to VLNC cigarettes (0.4 mg nicotine per gram
tobacco). From this trial, data to supplement CENIC-P1S2 include the usual
brand cigarette group as control and the VLNC groups as treatment (high and
low tar groups combined). CENIC-P2 was a 20-week RCT that randomized smokers
to one of three groups, 1) maintenence on NNC cigarettes, 2) immediate
reduction to VLNC cigarettes, of 3) gradual reduction to VLNC cigarettes. From
this trial, data to supplement CENIC-P1S2 include the NNC group as control and
the immediate-reduction group as treatment using outcomes collected at week 8.

One approach to combining data is simply to pool data from multiple sources to
increase the overall sample size. This naive approach would increase the
precision of the estimated effect of an intervention, but it leads to both
bias and inflated type 1 error rates in the presence of inter-source
heterogeneity. A more statistically rigorous approach is to use estimators
that “shrink" the estimated treatment effect from a primary data source
towards estimates of the treatment effect from supplementary sources, thus
borrowing strength and increasing precision. For example, we will have more
confidence in our estimate of the treatment effect in a small trial if it is
consistent with the treatment effect observed in the supplemental data source.

Statistical methods that facilitate the borrowing of information across data
sources arise most naturally in the Bayesian paradigm. These methods have been
used in the meta-analysis literature and have more recently been adapted to
the analysis of RCTs. Bayesian approaches to incorporating supplemental
information into the analysis of RCTs include the derivation of predictive
prior distributions using the standard hierarchical model (Gelman, 2006);
power priors, which use a pre-specified parameter that down-weights
supplemental data (Ibrahim and Chen, 2000); and commensurate priors, which
utilize a hierarchical modeling framework with a commensurability parameter to
control borrowing from supplemental sources (Hobbs _and others_ , 2011).

Kaizer _and others_ (2018) recently introduced multi-source exchangeability
models (MEMs), a novel approach to integrating supplementary data into the
analysis of a primary data source. MEMs offer a flexible approach to cross-
source data integration by including parameters that separately control
borrowing from each supplementary data source and can be estimated from the
data. This can be achieved by specifying models that represent all possible
ways that the supplementary data can be combined with the primary data source
and using Bayesian model averaging to obtain a weighted average of the
estimated treatment effect across all models. Simulation results suggest that
MEMs have favorable statistical properties when compared to existing
hierarchical modeling strategies (Kaizer _and others_ , 2018).

In this manuscript, we show how to use MEMs to borrow from supplemental
sources in the estimation of causal effects using a primary data source. Our
method uses regression models to estimate causal effects in the presence of
confounding. In Section 2, we discuss two models for estimating causal
effects, a Bayesian Linear Model (BLM), and Bayesian Additive Regression Trees
(BART)(Chipman _and others_ , 2010), and we should how to use these models to
borrow from supplemental data to estimate causal effects from a primary data
source in the MEMs framework. In Section 3, we present results of a simulation
study to assess the performance of the estimators. In Section 4 we estimate
the causal effect of VLNC cigarettes from CENIC-P1S2 allowing borrowing from
CENIC-P1S1 and CENIC-P2, and we conclude in Section 5.

##  2 Methods

###  2.1 Observed Data

We consider data from a single primary source $P$ and $H$ supplemental sources
with sample sizes $n_{P},n_{1},\dots,n_{H}$. The data are assumed to be random
samples from each source population. We assume the following data are
collected for the $i$th participant: a source variable $S_{i}$ that can take
values $P,1,\dots,H$, indicating the data source; the quantitative outcome
$Y_{i}$; a binary indicator of treatment assignment $A_{i}$, with untreated
(or control) participants having $A_{i}=0$ and $A_{i}=1$ indicating assignment
to treatment; and a vector of patient characteristics
$\boldsymbol{X}_{i}=\left(X_{i1},\dots,X_{ip}\right)^{T}$ associated with
$Y_{i}$ and possibly with $A_{i}$ if treatment is not randomly assigned. For
data source $s$, let
$\mathcal{D}_{s}=\left\\{\left(Y_{i},A_{i},\boldsymbol{X}_{i},S_{i}\right):S_{i}=s\right\\}$
represent the observed data.
$\mathcal{D}=\left(\mathcal{D}_{P},\mathcal{D}_{1},\dots,\mathcal{D}_{H}\right)$
denotes the collection of data from all sources.

###  2.2 Potential Outcomes and Target of Inference

Let $Y^{*}(a,s)$ denote the outcome of a randomly selected participant if,
possibly contrary to fact, $A=a$ and $S=s$. Because for each participant we do
not observe $Y^{*}(a,s)$ for all $a$ and $s$, $Y^{*}(a,s)$ is known as a
potential outcome. Our target of inference is the population average treatment
effect (PATE)
$\Delta_{P}=E_{\boldsymbol{X}|S=P}\left[E\left\\{Y^{*}(1,P)-Y^{*}(0,P)|\boldsymbol{X},S=P\right\\}\right]$,
the expected difference difference for being treated versus remaining
untreated in the primary source population, where the outer expectation is
with respect to the conditional distribution of $\boldsymbol{X}|S=P$.

###  2.3 Identifying Assumptions

We make the following assumptions to relate the distribution of the observed
data to the distribution of the potential outcome (Robins and Hernán, 2009).
First, we assume that the potential outcome under treatment $a$ in source $s$
is equal to the observed outcome given that $A=a$ and $S=s$, which can be
restated as
$Y=\sum_{s=P,1,\dots,H}\left\\{AY^{*}\left(1,s\right)+\left(1-A\right)Y^{*}\left(0,s\right)\right\\}I\left(S=s\right)$.
This is referred to as the consistency assumption. Second, we assume that we
have collected sufficient data on participant characteristics so that
treatment assignment is conditionally ignorable, or $\left.\left\\{A\perp
Y^{*}(0,S),Y^{*}(1,S)\right\\}\;\middle|\;\boldsymbol{X},S\right.$. Third, we
assume that $0<Pr\left(A=1|\boldsymbol{X},S\right)<1$ for all
$\boldsymbol{X}$. This is referred to as the positivity assumption.

###  2.4 Estimators

Under the assumptions stated above, we have the following equalities:
$E\left(Y|A=1,S,\boldsymbol{X}\right)=E\left\\{Y^{*}(1,S)|S,\boldsymbol{X}\right\\}$,
and
$E\left(Y|A=0,S,\boldsymbol{X}\right)=E\left\\{Y^{*}(0,S)|S,\boldsymbol{X}\right\\}$.
This allows us to work with the observed data only and to estimate
$\Delta_{P}$, provided that we have an estimator for
$E\left(Y|A,S,\boldsymbol{X}\right)$. Although many estimators are possible,
we describe two, a Bayesian linear model (BLM) and Bayesian additive
regression trees (BART). For now, we consider the simple case where the only
data available are from the primary source, and there is no borrowing. We
discuss borrowing from supplemental sources in Section 2.6.

####  2.4.1 Bayesian Linear Model

The first estimator we consider is the Bayesian linear model (BLM) with a
normal-inverse-gamma prior. Define the model matrix $\boldsymbol{D}$ with the
$i$th row equal to $\left(1,A_{i},\boldsymbol{X}_{i}^{T}\right)^{T}$, and let
$\boldsymbol{Y}=\left(Y_{i},\dots Y_{n_{P}}\right)^{T}$ denote the outcome
vector. The BLM is

$\boldsymbol{Y}=\boldsymbol{D}\boldsymbol{\beta}_{P}+\boldsymbol{\epsilon},\quad\boldsymbol{\epsilon}\sim
N\left(0,\sigma^{2}_{P}\boldsymbol{I}\right),$ (1)

with $\sigma^{2}_{P}\sim
IG\left(a_{P},b_{P}\right),\boldsymbol{\beta}_{P}|\sigma^{2}_{P}\sim
N\left(\boldsymbol{\mu}_{P},\sigma^{2}_{P}\boldsymbol{V}_{P}\right)$, where
$a_{P},b_{P},\boldsymbol{\mu}_{P},\text{ and }\boldsymbol{V}_{P}$ are known
hyperparameters. The subscript $P$ on the coefficients and parameters is meant
to clearly indicate that these are for modelling data in the primary source
$P$.

Our hyperparameter specification is similar to those of Raftery _and others_
(1997). We let $\boldsymbol{\mu}_{P}=\left(\bar{y},0,\dots,0\right)^{T}$ where
$\bar{y}$ is the sample mean of $Y_{1},\dots,Y_{n_{P}}$, and
$\boldsymbol{V}_{P}=\left(\frac{1}{n_{P}}\boldsymbol{D}^{T}\boldsymbol{D}\right)^{-1}$.
We choose $a_{P}$ and $b_{P}$ so that $\sigma^{2}_{P}$ has mean
$\hat{\sigma}^{2}$ and variance $2\hat{\sigma}^{4}$, where $\hat{\sigma}^{2}$
is from the least squares regression of $\boldsymbol{Y}$ on $\boldsymbol{D}$.

The primary advantage of the BLM is the resulting efficiency of the estimator
for $\Delta_{P}$. If Model (1) were the correct model, and if
$\boldsymbol{\beta}$ were estimated using the ordinary least squares (OLS)
estimator $\hat{\boldsymbol{\beta}}$, then the estimator for $\Delta_{P}$
implied by $\hat{\boldsymbol{\beta}}$ would be consistent and the most
efficient (Lunceford and Davidian, 2004). Although the efficiency of the
Bayesian estimator for $\Delta_{P}$ using Model (1) depends on the prior
specification, the estimator will have performance asymptotically equivalent
to the OLS estimator provided that a suitable posterior summary statistic
(e.g., the mean) is chosen as the estimator for $\Delta_{P}$. The primary
disadvantage of using of Model (1) is that specifying a correct, or
approximately correct, model is challenging if the dimension of
$\boldsymbol{X}$ is large, particularly if interactions or nonlinear
relationships are present.

####  2.4.2 Bayesian Additive Regression Trees

BART (Chipman _and others_ , 2010) is a nonparametric tree ensemble model
primarily used for prediction. Prior to fitting the model, Chipman _and
others_ (2010) shift and scale $Y$ so that the sample mean is 0 and the range
is 1; we assume throughout this section that $Y$ is on the transformed scale.
BART consists of two components, a sum-of-trees model, and a regularization
prior. The sum-of-trees model with $m$ trees is

$Y_{i}=\sum_{j=1}^{m}g\left(A_{i},\boldsymbol{X}_{i};T_{Pj},M_{Pj}\right)+\epsilon_{i},\quad\epsilon_{i}\sim
N\left(0,\sigma^{2}_{P}\right).$ (2)

where $T_{Pj}$ denotes the $j$th tree and
$M_{Pj}=\left(\mu_{Pj1},\dots,\mu_{Pjk_{j}}\right)$ denotes the terminal node
values of $T_{P_{j}}$. The function $g$ sends $A$ and $\boldsymbol{X}$ down
branches of each of the $m$ trees and sorts it based on interior node binary
decision rules of the form $x\leq c$ vs. $x>c$ until a terminal node is
reached. Each terminal node then assigns each observation one value of
$M_{Pj}$. Thus the conditional mean value of $Y$ is the sum of the terminal
node values to which the observation has been assigned based on values of $A$
and $\boldsymbol{X}$. The regularization prior ensures that each tree is a
“weak learner” and contributes a small amount to the overall fit. The terminal
node values are a priori independent and identically distributed as
$N\left(0,\tau^{2}\right)$, where $\tau^{2}$ is a constant. The prior for
$\sigma^{2}_{P}$ is an inverse gamma, $\sigma^{2}_{P}\sim
IG\left(\frac{\nu}{2},\frac{\nu\lambda}{2}\right)$, with $\nu$ and $\lambda$
both fixed. Chipman _and others_ (2010) use a data-informed approach to select
$\lambda$, and they recommend a default value of $\nu=3$. See Chipman _and
others_ (2010) for further details on the regularization prior and
specification of hyperparameters.

Hill (2011) proposed using BART for causal inference due its ease of use, its
flexibility in modelling complex response surfaces, and its ability to
identify heterogeneous treatment effects. Similar to the BLM, the analyst must
specify the outcome, treatment, and confounders. But unlike the BLM, the
analyst need not specify a functional form for the relationship between the
outcome and the predictors. Any interactions and non-linear effects that could
escape the attention of the analyst using a linear model can be captured quite
naturally through the fitted sum-of-trees model. The primary disadvantage of
BART is that it cannot capitalize on linear relationships to accurately
predict the potential outcome response surface where very little data are
available. BART is therefore not expected to perform well in scenarios with
strong confounding. The posterior credible intervals are naturally wider in
these areas, however, which partially ameliorates this problem (Hill, 2011).

###  2.5 Posterior Inference

Here we describe posterior inference on $\Delta_{P}$ assuming we have already
fit a model to the primary data source using either BART or the BLM for
$E\left(Y|A,S=P,\boldsymbol{X},\boldsymbol{\theta}_{P}\right)$ with model
parameters $\boldsymbol{\theta}_{P}$. For the BLM,
$\boldsymbol{\theta}_{P}=\left(\boldsymbol{\beta}_{P},\sigma^{2}_{P}\right)$;
for BART,
$\boldsymbol{\theta}_{P}=\left\\{\left(T_{Pj},M_{Pj}\right)_{j=1,\dots,m},\sigma^{2}_{P}\right\\}$.
Suppose we take $B$ draws from the posterior distribution of
$\boldsymbol{\theta}_{P}$, with $\boldsymbol{\theta}_{P}^{(b)}$ denoting the
$b$th draw. For each draw from the posterior, we compute for each
$\boldsymbol{X}$ in the primary data source the conditional average treatment
effect,

$E\left\\{Y|A=1,S=P,\boldsymbol{X},\boldsymbol{\theta}_{P}^{(b)}\right\\}-E\left\\{Y|A=0,S=P,\boldsymbol{X},\boldsymbol{\theta}_{P}^{(b)}\right\\},$
(3)

which is the difference in expected potential outcomes if the observation were
assigned treatment versus control. The PATE for the $b$th posterior draw is
then

$\Delta_{P}^{(b)}=E_{\boldsymbol{X}|S=P}\left[E\left\\{Y|A=1,S=P,\boldsymbol{X},\boldsymbol{\theta}_{P}^{(b)}\right\\}-E\left\\{Y|A=1,S=P,\boldsymbol{X},\boldsymbol{\theta}_{P}^{(b)}\right\\}\right].$

The outer expectation requires an integration over $\boldsymbol{X}|S=P$, the
distribution of covariates in the primary source population. To avoid
specifying a multivariate distribution for $\boldsymbol{X}|S=P$ with priors,
we approximate the integral using the Bayesian bootstrap (Rubin, 1981), which
has previously been used to integrate over the distribution of covariates
(Wang _and others_ , 2015; Nethery _and others_ , 2019). For each value of the
conditional average treatment effects in (3), we perform a single bootstrap
iteration as follows. Let
$\boldsymbol{p}^{(b)}=\left\\{p_{1}^{(b)},\dots,p_{n_{P}}^{(b)}\right\\}$
denote the vector of sampling probabilities for the $b$th posterior draw as
described by Rubin (1981). We let
$\Delta_{P}^{(b)}=\sum_{i=1}^{n_{P}}p_{i}^{(b)}\left[E\left\\{Y|A=1,S=P,\boldsymbol{X},\boldsymbol{\theta}_{P}^{(b)}\right\\}-E\left\\{Y|A=0,S=P,\boldsymbol{X},\boldsymbol{\theta}_{P}^{(b)}\right\\}\right]$.
The resulting collection $\\{\Delta_{P}^{(1)},\dots,\Delta_{P}^{(B)}\\}$ is
taken to be the posterior distribution of $\Delta_{P}$.

###  2.6 Multisource Exchangeability Models

Having described estimators for the $\Delta_{P}$ using the primary source
only, our goal is now to incorporate information from supplemental data
sources into making inference about $\Delta_{P}$. To do this, we use
Multisource Exchangeability Models (MEMs), a Bayesian model averaging (BMA)
technique first proposed by Kaizer _and others_ (2018). We begin by selecting
a model, either BLM or BART. Let $\boldsymbol{\theta}_{h}$ denote the source-
specific model parameters for the supplemental data source $h$. We can borrow
information from supplemental data source $h$ by assuming that the data source
is exchangeable with the primary data source. In this context, data source $h$
is exchangeable with $P$ if the model parameters are identical to the model
parameters for the primary data source $P$. As each of the $H$ data sources
can be exchangeable, or not, with the primary data source, there are $2^{H}$
possible patterns of exchangeability. Let $Z_{h}$ denote an indicator of
whether source $h$ is exchangeable with source $P$. The $q$th pattern of
exchangeability, corresponding to the $q$th MEM, is denoted as
$\boldsymbol{\Omega}_{q}=\left(Z_{1}=z_{1q},\dots,Z_{H}=z_{Hq}\right)$. The
model parameters for source $h$ under $\boldsymbol{\Omega}_{q}$ are
$z_{1q}\boldsymbol{\theta}_{P}+\left(1-z_{1q}\right)\boldsymbol{\theta}_{h}$.
For simplicity, the priors for
$\boldsymbol{\theta}_{P},\boldsymbol{\theta}_{1},\dots,\boldsymbol{\theta}_{H}$
are independent so that the prior for $\boldsymbol{\Theta}_{q}$, the
parameters of $\boldsymbol{\Omega}_{q}$, is
$p\left(\boldsymbol{\Theta}_{p}|\Omega_{p}\right)=p\left(\boldsymbol{\theta}_{P}\right)\prod_{h=1}^{H}p\left(\boldsymbol{\theta}_{h}\right)^{\left(1-z_{hp}\right)}$.

The marginal likelihood is obtained by averaging the likelihood over the
prior:

$p\left(\mathcal{D}|\boldsymbol{\Omega}_{q}\right)=\int
p\left(\mathcal{D}|\boldsymbol{\Theta}_{q},\boldsymbol{\Omega}_{q}\right)p\left(\boldsymbol{\Theta}_{q}|\boldsymbol{\Omega}_{q}\right)d\boldsymbol{\Theta}_{q}.$
(4)

For simplicity, we do not specify a model for the joint distribution for
$Y,A,S,\boldsymbol{X}$, but only for the conditional distribution of $Y$. Thus
the likelihood
$p\left(\mathcal{D}|\boldsymbol{\Theta},\boldsymbol{\Omega}_{p}\right)$ in
Equation (4) is simply a product of conditional normal densities given the
model for $Y$.

The posterior probability that $\boldsymbol{\Omega}_{q}$ is the correct model,
given the data, is

$\omega_{q}=p\left(\boldsymbol{\Omega}_{q}|\mathcal{D}\right)=\frac{p\left(\mathcal{D}|\boldsymbol{\Omega}_{q}\right)p\left(\boldsymbol{\Omega}_{q}\right)}{\sum_{w=1}^{2^{H}}p\left(\mathcal{D}|\boldsymbol{\Omega}_{w}\right)p\left(\boldsymbol{\Omega}_{w}\right)},$

where $p\left(\boldsymbol{\Omega}_{p}\right)$ is the prior probability that
$\boldsymbol{\Omega}_{p}$ is the true model. We explore prior probabilities on
the MEMs in simulation studies of Section 3.

Finally, for each MEM, data sources that share $\boldsymbol{\theta}_{P}$ are
simply combined for estimation as they would be in a traditional analysis.
Borrowing is thus accomplished, under each MEM, by combining the exchangeable
data sources to derive the MEM-specific model parameter posterior
distributions for
$\boldsymbol{\theta}_{P}^{\boldsymbol{\Omega}_{1}},\dots,\boldsymbol{\theta}_{P}^{\boldsymbol{\Omega}_{2^{H}}}$
and the concomitant MEM-specific PATE posteriors for
$\Delta_{P}^{\boldsymbol{\Omega}_{1}},\dots,\Delta_{P}^{\boldsymbol{\Omega}_{2^{H}}}$.The
posterior for $\Delta_{P}$ is then a weighted average of the posterior
distributions given each MEM,

$p\left(\Delta_{P}|\mathcal{D}\right)=\sum_{q=1}^{2^{H}}\omega_{q}p\left(\Delta_{P}^{\boldsymbol{\Omega}_{q}}|\mathcal{D},\boldsymbol{\Omega}_{q}\right).$
(5)

Draws from the posterior of
$p\left(\Delta_{P}^{\boldsymbol{\Omega}_{q}}|\mathcal{D},\boldsymbol{\Omega}_{q}\right)$
proceed as described above, with computation of conditional expected outcomes
and integration over $\boldsymbol{X}|S=P$ using only data from the primary
source and the posterior for
$\boldsymbol{\theta}_{P}^{\boldsymbol{\Omega}_{q}}$, the parameters for the
primary source under the exchangeability assumption for
$\boldsymbol{\Omega}_{q}$.

###  2.7 Marginal Likelihoods

Because the likelihood
$p\left(\mathcal{D}|\boldsymbol{\Theta}_{q},\boldsymbol{\Omega}_{q}\right)$ is
a product of conditional densities of $Y$ rather than a joint distribution for
the observed data, derivation of the marginal likelihood in Equation (4)
treats $A,S$, and $\boldsymbol{X}$ as fixed. For simplicity we derive the
marginal likelihood for the primary data source only for the MEM with
$Z_{h}=0$ for all $h$, indicating no borrowing between data sources. For MEMs
with borrowing, data sources with shared parameters $\boldsymbol{\theta}_{P}$
are combined for computing (4).

####  2.7.1 Bayesian Linear Model

For Model (1) with the normal-inverse gamma prior, the marginal likelihood has
the multivariate $t$-density
$p\left(\boldsymbol{Y}|\boldsymbol{D}\right)\sim
t_{2a_{P}}\left\\{\boldsymbol{D}\boldsymbol{\mu}_{P},\frac{b_{P}}{a_{P}}\left(\boldsymbol{I}+\boldsymbol{D}\boldsymbol{V}_{P}\boldsymbol{D}^{T}\right)\right\\}$,
where $\boldsymbol{I}$ is the $n_{P}\times n_{P}$ identity matrix.

####  2.7.2 BART

Marginal Likelihood Under Default Prior. To simplify the derivation, let
$\mathcal{T}_{P}=\\{T_{P1},\dots,T_{Pm}\\}$ and
$\mathcal{M}_{P}=\\{M_{P1},\dots,M_{Pm}\\}$ denote the set of trees and
terminal nodes, respectively. The marginal likelihood given model (2) under
the default prior specification is

$\displaystyle p\left(\boldsymbol{Y}|\boldsymbol{D}\right)=$
$\displaystyle\sum\int\int
p\left(\boldsymbol{Y}|\boldsymbol{D},\mathcal{M}_{P},\mathcal{T}_{P},\sigma^{2}_{P}\right)p\left(\mathcal{M}_{P}|\mathcal{T}_{P}\right)p\left(\mathcal{T}_{P}\right)p\left(\sigma^{2}_{P}\right)d\mathcal{M}_{P}d\sigma^{2}_{P},$

where the summation is over the discrete distribution of trees. To integrate
out $\mathcal{M}_{P}$, note that the terminal nodes are a priori independent
and identically distributed as $N\left(0,\tau^{2}\right)$. This implies the
prior $E\left(\boldsymbol{Y}|\boldsymbol{D},\mathcal{T}_{P}\right)\sim
N\left(\boldsymbol{0},m\tau^{2}\boldsymbol{R}_{P}\right)$, where
$\boldsymbol{R}_{P}$ is a correlation matrix with off-diagonal elements
$r_{kl}$ equal to the proportion of terminal nodes shared between
$E\left(Y_{k}|A_{k},\boldsymbol{X}_{k}\right)$ and
$E\left(Y_{l}|A_{l},\boldsymbol{X}_{l}\right)$. Hence,
$\boldsymbol{Y}|\boldsymbol{D},\mathcal{T}_{P},\sigma^{2}_{P}\sim
N\left(\boldsymbol{0},\boldsymbol{\Sigma}_{P}\right)$, where
$\boldsymbol{\Sigma}_{P}=m\tau^{2}\boldsymbol{R}_{P}+\sigma^{2}_{P}\boldsymbol{I}$.
Using this result, we can re-write the marginal likelihood as

$p\left(\boldsymbol{Y}|\boldsymbol{D}\right)\propto\sum\int\lvert\boldsymbol{\Sigma}_{P}\rvert^{-\frac{1}{2}}e^{-\frac{1}{2}\boldsymbol{Y}^{T}\boldsymbol{\Sigma}_{P}^{-1}\boldsymbol{Y}}\left(\frac{1}{\sigma^{2}_{P}}\right)^{\frac{\nu}{2}+1}e^{-\frac{\nu\lambda}{2\sigma^{2}_{P}}}p\left(\mathcal{T}_{P}\right)d\sigma^{2}_{P},$

where we have substituted the kernels of the normal probability density
function (pdf) for
$p\left(\boldsymbol{Y}|\boldsymbol{D},\mathcal{T}_{P},\sigma^{2}_{P}\right)$
and the inverse gamma pdf for $p\left(\sigma^{2}_{P}\right)$. At this point we
would integrate over $\sigma^{2}_{P}$, but this is not (analytically)
straightforward:
$\boldsymbol{\Sigma}_{P}=m\tau^{2}\boldsymbol{R}_{P}+\sigma^{2}_{P}\boldsymbol{I}$
is a function of $\sigma^{2}_{P}$, but $\sigma^{2}_{P}$ cannot be factored out
of $\boldsymbol{\Sigma}_{P}$. For this reason, we revert to the prior on the
terminal node values
$\mu_{jk}\overset{iid}{\sim}N\left(0,\frac{\sigma^{2}_{P}}{\gamma}\right)$
originally used by Chipman _and others_ (1998) in their Bayesian
implementation of classification and regression trees (see also Hernández _and
others_ (2018)).

Marginal Likelihood Under Modified Prior. Assume that the terminal nodes
values are a priori independent and identically distributed as
$N\left(0,\frac{\sigma^{2}_{P}}{\gamma}\right)$, where $\gamma$ is a constant.
Under this modified prior, the marginal likelihood is

$p\left(\boldsymbol{Y}|\boldsymbol{D}\right)\propto\sum\int\lvert\boldsymbol{\Sigma}_{P}^{*}\rvert^{-\frac{1}{2}}e^{-\frac{1}{2}\boldsymbol{Y}^{T}\boldsymbol{\Sigma}_{P}^{*-1}\boldsymbol{Y}}\left(\frac{1}{\sigma^{2}_{P}}\right)^{\frac{\nu}{2}+1}e^{-\frac{\nu\lambda}{2\sigma^{2}_{P}}}p\left(\mathcal{T}_{P}\right)d\sigma^{2}_{P},$

where
$\boldsymbol{\Sigma}_{P}^{*}=\frac{m\sigma^{2}_{P}}{\gamma}\boldsymbol{R}_{P}+\sigma^{2}_{P}\boldsymbol{I}$.
Note that $\sigma^{2}_{P}$ can be factored out of
$\boldsymbol{\Sigma}_{P}^{*}$:
$\boldsymbol{\Sigma}_{P}^{*}=\sigma^{2}\boldsymbol{U}_{P}$, where
$\boldsymbol{U}_{P}=\frac{m}{\gamma}\boldsymbol{R}_{P}+\boldsymbol{I}$. The
marginal likelihood can now be simplified to

$\displaystyle p\left(\boldsymbol{Y}|\boldsymbol{D}\right)\propto$
$\displaystyle\sum\int\lvert\sigma^{2}_{P}\boldsymbol{U}_{P}\rvert^{-\frac{1}{2}}e^{-\frac{1}{2\sigma^{2}_{P}}\boldsymbol{Y}^{T}\boldsymbol{U}_{P}^{-1}\boldsymbol{Y}}\left(\frac{1}{\sigma^{2}_{P}}\right)^{\frac{\nu}{2}+1}e^{-\frac{\nu\lambda}{2\sigma^{2}_{P}}}p\left(\mathcal{T}_{P}\right)d\sigma^{2}_{P}$
$\displaystyle=$
$\displaystyle\sum\lvert\boldsymbol{U}_{P}\rvert^{-\frac{1}{2}}\left\\{\int\left(\frac{1}{\sigma^{2}_{P}}\right)^{\frac{n_{P}+\nu}{2}+1}e^{-\frac{1}{2\sigma^{2}_{P}}\left(\boldsymbol{Y}^{T}\boldsymbol{U}_{P}^{-1}\boldsymbol{Y}+\nu\lambda\right)}d\sigma^{2}_{P}\right\\}p\left(\mathcal{T}_{P}\right)$
$\displaystyle\propto$
$\displaystyle\sum\left[\lvert\boldsymbol{U}_{P}\rvert^{-\frac{1}{2}}\left\\{1+\frac{1}{\nu}\boldsymbol{Y}^{T}\left(\lambda\boldsymbol{U}_{P}\right)^{-1}\boldsymbol{Y}\right\\}^{-\frac{\nu+n_{P}}{2}}\right]p\left(\mathcal{T}_{P}\right)$

The term in the square brackets is the density of the central multivariate $t$
distribution with $\nu$ degrees of freedom and shape matrix
$\lambda\boldsymbol{U}_{P}$. Note that in the first equality,
$\det{\left(\boldsymbol{U}_{P}\right)}^{-\frac{1}{2}}$ cannot be treated as a
proportionality constant, because $\boldsymbol{R}_{P}$ depends on the
structure of the $m$ trees. Simplifying the above expression,

$p\left(\boldsymbol{Y}|\boldsymbol{D}\right)=\sum
t_{\nu}\left(\boldsymbol{Y}|\boldsymbol{0},\lambda\boldsymbol{U}_{P}\right)p\left(\mathcal{T}_{P}\right).$

Under the modified prior, the marginal density of $\boldsymbol{Y}$ is a
weighted average of multivariate $t$ densities.

Although this modified prior is not required for estimation or for computing
the marginal likelihood, it does make computing the marginal likelihood more
straightforward. Furthermore, BART is traditionally viewed as a “black box”
method in that in that it takes input and produces output with little
intuition for the intermediate steps (but see Tan and Roy (2019) for a lucid
description of BART). At minimum, this result provides some intuition for how
the posterior weights come about, particularly because the $t$-distribution is
well-characterized and lends itself to analysis.

##  3 Simulation

We conducted a simulation study to assess the performance of the BLM and BART
while borrowing from a supplemental source. Part 1 demonstrates the borrowing
and performance of the estimators under several simple scenarios. Part 2 is
designed to test the performance of the estimators under a variety of
realistic scenarios. All models were fit using the `borrowr R` package
(Boatman _and others_ , 2019), which includes an implementation of the BLM and
BART using a modified version of the `R BART` package to draw from the
posterior distributions implied by the modified prior on the terminal nodes
discussed in Section 2.4.2 (Hernández _and others_ , 2018).

###  3.1 Simulation Part 1

We consider data from a primary source and a single supplemental source, each
with sample size 100, under three scenarios. For all scenarios,
$\Pr\left(A=1\right)=0.5$ (marginally). For the primary data source the
outcome is generated according to the model
$Y=A+f\left(X\right)+\epsilon,\epsilon\sim N(0,1)$, with $f$ and $X$ defined
below for each scenario. In all scenarios the PATE is fixed at 1. In the
supplemental data source the outcome is generated as
$Y=\left(1+\delta\right)A+f\left(X\right)+\epsilon,\epsilon\sim N(0,1)$, with
$\delta$ varied along a grid of values from -2.5 to 2.5 in increments of 1/2.
In all scenarios, the BLM regresses $Y$ on $A$ and $X$ with no interactions or
non-linear terms. For BART, we used the default priors $m=200,k=2,\nu=3$, with
100 burn-in posterior draws, and 100 draws from the posterior for inference.
The marginal likelihood for BART was estimated by averaging the likelihood
over 100 draws from the prior distribution. For each iteration $\gamma$ was
set to $\frac{1}{16m\hat{\sigma}^{2}}$, where $\hat{\sigma}^{2}$ is estimated
from a least squares regression of $Y$ on $A$ and $X$, so that
$\gamma\sigma^{2}$ is comparable to the unmodified BART prior variance of
$\frac{1}{16m}$. We include a “no borrowing” (NB) estimator that adjusts for
confounding but does not borrow from supplemental sources. This is a BLM that
regresses $Y$ on $A$ and $X$ but uses only the primary data source. We
investigated two priors on models. With only one supplemental source, the
prior is defined by the prior probability of exchangeability, $\Pr(Z_{1}=1)$.
The first prior is $\Pr(Z_{1}=1)=\frac{1}{2}$, and the second is
$\Pr(Z_{1}=1)=\left(\frac{1}{2}\right)^{r}$, where $r$ is the number of
predictors in the model ($+1$ for the intercept for the BLM). The former is
shown in figures with no subscript (“BLM" and “BART"). The latter is meant to
discourage borrowing as the dimension of $X$ grows and is denoted in figures
with the subscript $r$ (“BLMr" “and BARTr"). For each scenario we report the
estimated bias, root mean squared error (MSE), and the posterior weight
corresponding to the MEM in favor of borrowing (that is, parameters between
the two data sources are shared) for 1,000 Monte Carlo simulations for each
value of $\delta$ within each scenario.

####  3.1.1 Scenario 1: Parallel Linear Response Surfaces

For this scenario we generate $X\sim N(0,1)$, $A|X\sim
Ber\left(\left(1+e^{-x}\right)^{-1}\right)$, and use $f(x)=x$. The top row of
Figure 1 shows the results. All estimators except NB borrowing across $\delta$
values of roughly -1.5 to 1.5, with BART borrowing the most and BLMP the
least. All estimators are minimally biased when not borrowing, and show small
amounts of bias when borrowing with $\delta\neq 0$. All estimators show
roughly the same reduction in MSE when borrowing close to $\delta=0$. As
$\delta$ moves from 0, all estimators begin to show increased MSE relative to
MSE when not borrowing, but this penalty is noticeably smaller for BLMr than
for BLM, and the same patter holds for BARTr compared to BART. In general the
root MSE for BART is greater than for the BLM.

Part 1, Scenario 1  Part 1, Scenario 2  Part 1, Scenario 3

Figure 1: Simulation Part 1 Results. Scenario 1: parallel, linear response
surfaces. Scenario 2: different mean of confounder in supplemental source.
Scenario 3: over-specified outcome model. $\delta$: Difference in treatment
effect between sources (supplemental - primary) on effect size scale. BLM:
Bayesian linear model; BLMP: Bayesian linear model with power prior; NB: no
borrowing.

####  3.1.2 Scenario 2: Parallel Linear Response Surfaces with Different
Confounder Distribution in the Supplemental Source

We generate data for the primary source as for Scenario 1, except that in the
supplemental source we generate $X\sim N\left(3,1\right)$. The purpose of the
scenario is to investigate the extent to which an investigator can benefit
from including supplemental sources that are drawn from sub-populations with
differing distributions of confounders than the primary data source. The
middle row of Figure 1 shows the results. The performace of BLM is similar to
Scenario 1, but BART shows almost no MSE reduction for $\delta$ close to 0,
and almost no MSE increase for $\delta$ values farther from 0.

####  3.1.3 Scenario 3: Non-Linear Response Surfaces with Different
Conditional Variances in Confounder

This scenario was informed by previous work (Rubin, 1973) indicating that
linear estimators can perform very poorly when $f(x)$ is non-linear and the
conditional variances of $X|A=0$ and $X|A=1$ are different. For both data
sources, we generate $A\sim Ber(0.5)$ and $X$ according to the conditional
distributions $X|A=0\sim N\left(0,\frac{4}{3}\right)$ and $X|A=1\sim
N\left(0,\frac{2}{3}\right)$, where the normal distributions are parameterized
with the variance. The outcome is generated using $f(x)=e^{x}$. As shown in in
the bottom row of Figure 1, this is a challenging scenario for both
estimators. The posterior weights indicate that BART borrows across a wide
range of $\delta$, while the weights for the BLM are lower at $\delta=0$ than
in both previous scenarios. The BLM is biased and the MSE is higher than BART
for all values of $\delta$. Although negative bias is evident for BART at
$\delta=0$, the smaller magnitude bias and MSE of BART clearly indicates that
BART performs better than the BLM in this scenario.

###  3.2 Simulation Part 2: Atlantic Causal Inference Competition Simulation

In the second part of the simulation, we applied the estimators to data from
the 2016 Atlantic Causal Inference Conference (ACIC) as described by Dorie
_and others_ (2019). The goal of the part of the simulation is to assess the
performance of the estimators under a variety of data-generating scenarios
more realistic than those presented in part 1. Briefly, ACIC generates
simulated data according to 77 scenarios, with simulation parameters varied
across 6 “knobs": degree of non-linearity in the response surface; degree of
non-linearity in the treatment surface; proportion treated; overlap, defined
as the degree to which areas of the covariate space have only untreated
individuals; alignment, defined as the degree to which covariates appear in
both the treatment and the response model; and treatment effect heterogeneity,
or departures from parallel response surfaces between treated and untreated
observations. Covariates are from a real data set and include 3 categorical, 5
binary, 27 count, and 23 continuous variables. See Dorie _and others_ (2019)
for additional details. The design matrix with dummy coding for categorical
predictors had 81 predictors excluding the intercept term. For each iteration
we randomly selected 400 observations and assigned 200 to be primary and 200
to be supplemental sources. The first 100 in each were treated, and the
remaining were control.

In addition to the model priors described above, we included for BART only
additional model priors that discouraged borrowing to varying degrees as a
function of $r$, the number of predictors: BARTlog2r, with model prior
$\Pr(Z_{1}=1)=\left(\frac{1}{2}\right)^{\log_{2}r}=\frac{1}{r}$; and BARTr2,
with model prior $\Pr(Z_{1}=1)=\left(\frac{1}{2}\right)^{\frac{r}{2}}$. Pilot
work suggested that BARTr, with model prior
$\Pr(Z_{1}=1)=\left(\frac{1}{2}\right)^{r}$ discouraged borrowing too strongly
in this high-dimensional setting, so we excluded BART with this model prior.
Futhermore, the performance of this implmentaion of the BLM in this simulation
was uniformly inferior to BART, so we show results considered BLMr with no
other model priors. We again included a “no borrowing” estimator that uses
only the primary day. For this part the no borrowing estimator is denoted
BARTNB. All potential predictors were included for BART. The BLM regressed the
outcome on all variables without considering possible interactions, non-
linear, or polynomial terms. We increased the number of posterior draws to
1,000 due to the more complicated response surfaces compared to part 1.

Figure 2 shows the results, with each point in the boxplots showing the mean
over 50 Monte Carlo iterations. BART is minimally biased when not borrowing
and when borrowing under each of the model priors, but the BLM is signifcantly
positively biased. All BART estimator that borrow show reduced root MSE
compared to BARTNB, with smaller reductions seen for model priors that more
strongly discourage borrowing. In contrast, the root MSE of BLMr is quite high
and performs much worse than all BART estimators.

Part 1, Scenario 2

Figure 2: ACIC Simulation Results. Each point is the mean for 50 Monte Carlo
iterations for one data-generating scenario. NB: no borrowing. Other
subscripts indicate different model prior probabilities.

We also varied the treatment effect in the supplemental source to to evaluate
the performance of the model priors. Similar to part 1, we varied the
treatment effect in the supplemental source by adding $\delta\times sd(Y)$,
where $sd(Y)$ is the standard deviation of the outcome before modifying the
treatment effect. $\delta$ was varied along a grid from -1.5 to 1.5 increments
of 0.75. Because the BLM performed poorly in Figure 2, we consider only BART
with various model priors. Figure 3 shows the results. Note that the boxplots
at 0 simply duplicate the results in Figure 2. When the standardized
difference in treatment effect is $\pm 1.5$, there is no borrowing and the
performance of all estimators is roughly equivalent. At $\pm 0.75$, increased
bias and higher root MSE is evident. The magnitude of bias and root MSE are
higher using model priors that do not discourage borrowing based on the the
number of predictors (BART), and lowest for models that do (BARTr2), with
BARTlog2r showing intermediate magnitude of bias and root MSE.

Figure 3: ACIC Simulation Results. Each point is the mean for 50 Monte Carlo
iterations for one data-generating scenario. NB: no borrowing. Other
subscripts indicate different model prior probabilities.

###  3.3 Simulation Summary

The simulation results suggest a simple conclusion. Part 1 suggests that if we
are reasonably confident we can specify a correct, or approximately correct,
outcome model, BLMr is the best estimator. Part 2 suggests that if the
dimension of the predictors is large, or if we’re not confident we can
correctly specify the outcome model, BART is the best estimator. Furthermore,
part 2 suggests better performance using model priors that discourage
borrowing as the dimension of the predictors grows. Although the potential
root MSE reduction is limited with these priors, the potential root MSE
penalty is much lower as well.

##  4 Application

We analyzed data from 3 RCTs evaluating the effect of VLNC cigarettes on
cigarettes smoked per day (CPD), using CENIC project 1, study 2 (P1S2) as the
primary source and project 1, study 1 (P1S1) and project 2 (P2) as
supplemental sources. P1S1 was a 6-week trial that randomly assigned
participants to one of seven groups, consisting of a usual brand group and
groups assigned experimental cigarettes with nicotine content ranging from
15.8 mg nicotine/gram tobacco to 0.4 mg nicotine/gram tobacco, or a group
assigned high-tar cigarettes with 0.4 mg nicotine/gram tobacco (hereafter
abbreviated mg/g) (Donny _and others_ , 2015). For this application the 15.8
mg/g group is considered control and the 0.4 mg/g groups, high and low tar
combined, are considered treatment. P1S2 was a 6-week, 2x2 factorial RCT that
randomly assigned current smokers to smoke experimental cigarettes with
nicotine content of 15.8 mg/g or 0.4 mg/g, with or without nicotine
replacement therapy (NRT) (Smith _and others_ , 2019). For this application,
the 15.8 mg/g group without NRT is the control group, and the 0.4 mg/g group
without NRT is the treatment group. In contrast to P1S1 and P1S2, P2 was a
20-week RCT designed to test gradual vs. immediate nicotine reduction
(Hatsukami _and others_ , 2018). Control participants were assigned to smoke
experimental cigarettes with nicotine content of 15.5 mg/g; gradual reduction
participants were assigned to smoke experimental cigarettes with nicotine
content gradually reduced from 15.5 mg/g to 0.4 mg/g over the trial, and the
immediate reduction group was assigned to immediately begin smoking cigarettes
with nicotine content of 0.4 mg/g. Participants were assessed mid-trial at 8
weeks rather than 6 as in P1S1 and P1S2, so we use the outcome data from week
8. For this application the 15.5 mg/g group was treated as control, and the
immediate reduction group was treated as the treatment group.

In these 3 trials, treatment was randomly assigned. However, causal effect
estimators that adjust for confounding are still necessary due to
noncompliance: participants self-reported whether they smoked non-study
cigarettes, and if so, were considered noncompliant. Although treatment
assignment is randomized, compliance to treatment is not, and treatment-
control comparisons among compliant participants are therefore confounded. We
considered participants compliant if they reported smoking 0 non-study
cigarettes per day and if their total nicotine equivalents was less than 6.41,
as values exceeding this threshold are inconsistent with compliance (Denlinger
_and others_ , 2016). Participants in control groups are smoking
investigational cigarettes, but these contain normal amounts of nicotine.
Participants in control groups who smoke non-study cigarettes are noncompliant
from a strict perspective, but as the study cigarettes are similar to usual
brand cigarettes, we consider all participants in the control groups compliant
and take their outcome to be study CPD plus non-study CPD.

The estimators described can accommodate this situation with minimal changes.
Let $Q$ be the number of non-study cigarettes smoked per day. A participant if
considered compliant if $Q=0$ and non-compliant otherwise. Consider expanded
notation for the potential outcomes, with $Y^{*}(a,s,q)$ indicating the
potential outcome if, possibly contrary to fact, $A=a,S=s$, and $Q=q$. Also
define the compliance indicator $C=I\left(Q=0\right)$. Our target of inference
is the expected difference in the number of study cigarettes smoked per day
for treatment versus control in the primary data source assuming all
participants were to be compliant, or
$\Delta_{P}=E_{\boldsymbol{X}|S=P}\left[E\left\\{Y^{*}(1,P,0)-Y^{*}(0,P,0)|\boldsymbol{X},S=P\right\\}\right]$.
Valid inference requires the following modifications to the assumptions in
Section 2.3. We assume that if $A=a,S=s$ and $Q=0$, then $Y=Y^{*}(a,s,0)$,
that is, if a participant is compliant with the study protocol, then her
observed outcome is identical to the potential outcome for $q=0$ (consistency
assumption). We further assume that we have collected sufficient detail on
participant characteristics so that compliance is conditionally ignorable, or
$\left.\left\\{C\perp
Y^{*}(0,S,0),Y^{*}(1,S,0)\right\\}\;\middle|\;\boldsymbol{X},S\right.$
(conditional ignorability and no unmeasured confounders). Under these
assumptions we have
$E\left\\{Y^{*}(1,S,0)|S,\boldsymbol{X}\right\\}=E(Y|A=1,C=1,S,\boldsymbol{X})$
and
$E\left\\{Y^{*}(0,S,0)|S,\boldsymbol{X}\right\\}=E(Y|A=0,C=1,S,\boldsymbol{X})$.
The expectations on the right side of the equalities can be estimated with the
BLM and BART as described above.

To identify variables associated with compliance, we fit a lasso regression
model (Tibshirani, 1996) with treatment indicators for treatment ($A$) and
compliance ($C$); dummy variables for data source, age, race, educational
attainment; and a variety of other potential confounders $X$ including
variables measuring level of addiction, level of dependence, and satisfaction
with and craving for usual brand cigarettes. The model indicated non-zero
coefficients for $A$, $C$, source, age, baseline CPD, and dependence
(Fagerstrom Test for Nicotine Dependence; FTND). To identify possible
interactions, we next fit a least squares model with additive linear terms for
these variables as well as all pairwise interactions excluding terms involving
a treatment-compliance interaction. (An interaction between $A$ and $C$ cannot
be evaluated in this case because all participants in the control group were
treated as compliant.) The model indicated significant interactions between
$A$ and age, between $C$ and baseline CPD, and between $A$ and baseline CPD.
The Bayesian model used additive linear terms for these variables as well as
the interactions identified in the last squares model. For BART, we used
$A,C$, and all baseline variables with no screening for variable selection.

Prior to implementing the proposed estimators, we fit the BLM and BART
separately for each data source and estimated the causal effect of treatment
within each source. This step is not required, but we report these results to
provide intuition for whether borrowing from P1S1 and P2 is appropriate. Table
1 summarizes the samples sizes (after removing observations with missing data
for outcomes and potential confounders), BLM coefficients, and causal effects
estimated separately within each data source. The model coefficients are
similar between data sources, indicating that borrowing from both
supplementary sources may be appropriate. In addition, the causal effects
estimated within each source are also similar. Although the BLM and BART
borrow based on assuming conditional exchangeability, and not on the marginal
causal effects, the similarity in the causal effects nevertheless makes us
more confident that borrowing from P1S1 and P2 is appropriate.

Table 1: Sample Sizes, Posterior mean (sd) of the Bayesian linear model coefficients, and estimated treatment effects ($\Delta$) estimated for each data source without borrowing. | P1S2 | P1S1 | P2
---|---|---|---
Sample Sizes |  |  |
Control | 54 | 105 | 216
Treatment | 48 | 218 | 375
Total | 102 | 323 | 591
Coefficients |  |  |
Intercept | 18.22(0.65) | 16.27(0.44) | 16.52(0.29)
$A$ | -16.95(7.18) | -13.39(4.07) | -11.14(2.56)
$C$ | -8.41(4.94) | -7.20(3.04) | -1.89(1.69)
Age | 0.00(0.08) | -0.01(0.06) | 0.01(0.04)
Baseline CPD | 0.52(0.27) | 0.33(0.22) | 0.80(0.11)
FTND | 0.63(0.41) | 0.49(0.24) | 0.51(0.17)
$A\times$ Age | 0.20(0.12) | 0.07(0.07) | 0.16(0.05)
$C\times$ Baseline CPD | 0.58(0.24) | 0.64(0.20) | 0.25(0.09)
$A\times$ Baseline CPD | 0.06(0.22) | 0.33(0.21) | -0.09(0.08)
$\bf{\boldsymbol{\Delta_{P}}}$ (Treatment - Control) |  |  |
BLM | -6.35(2.19) | -5.37(1.43) | -5.03(0.75)
BART | -5.91(2.02) | -6.04(1.29) | -4.96(0.75)

For models with borrowing we assumed the model prior
$\Pr(Z_{h}=1)=\left(\frac{1}{2}\right)^{r}$ with $r=39$ predictors. Table 2
shows the results of the primary analysis. For the BLM, the the most weight is
given to models assuming both P1S1 and P2 are exchangeable with P1S2, and
assuming only P2 is. In contrast, BART essentially indicates that only P2 is
exchangeable with P12. From the BLM, the posterior mean(SD) reduction in study
cigarettes smoked per day for being assigned to treatment is 5.08(0.68) as
estimated from the BLM; the BART estimate is very similar at 4.91(0.70),
suggesting that the BLM is sufficiently complex to model the response surface.
Both estimators show dramatic reductions in the posterior variance of the
causal effect compared to no borrowing. It is counterintuitive that the BART
estimated causal effect is not within the range of effects separated
separately within each source, but as borrowing is on conditional means and
not on marginal means, it will not necessarily be the case that the marginal
effect will be shrunk toward either marginal treatment effect.

Table 2: Posterior model weights ($\omega$) and mean(sd) of the treatment
effect ($\Delta_{P}$) for the Bayesian linear model (BLM) and Bayesian
additive regression trees (BART). _Exchangeability and Posterior MEM Weights_
---
| Exchangeable with P1S2 |  | $\omega$
MEM | P1S1 | P2 |  | BLM | BART
1 | Yes | Yes |  | 0.6343 | 0.0003
2 | No | Yes |  | 0.3652 | 0.9996
3 | Yes | No |  | 0.0005 | 0.0000
4 | No | No |  | 0.0000 | 0.0001
_Weighted Posterior Mean(SD) Treatment Effects_
|  |  |  | BLM | BART
$\Delta_{P}$ (Treatment - Control) |  | -5.08(0.68) | -4.91(0.70)

##  5 Discussion

In this manuscript we demonstrated how to borrow from supplementary data
sources to estimate causal effects from a primary data source. The simulation
demonstrates that both the BLM and BART are viable estimators with good
properties depending on the level of confounding and on the complexity of the
response surface. The primary advantages of the BLM are that it results in
interpretable coefficients, and that the variance of the causal estimator may
be lower than for BART due to its frequentist operating characeristics. The
primary advantage of BART is that the analyst is not required to specify a
function for for the covariate-response association, meaning that BART can
naturally capture complicated relationships among the covariates and outcome.

The primary drawback of our proposed estimators is that they rely on
regression models to estimate causal effects. A major disadvantage of this is
that, in cases with incomplete overlap in the covariate distribution between
treatment and control groups, the estimated response surfaces in many areas
are based on extrapolation (Lunceford and Davidian, 2004). For this reason,
regression models are sometimes avoided for causal inference, and other
methods, particularly those based on propensity scores, are favored. One of
the primary advantages of the propensity score is that it allows the analyst
to condition on the propensity score in the outcome model rather than
specifying a complete covariate-outcome model. This is problematic when
attempting to borrow from supplemental sources, however. Borrowing between
data sources happens most naturally in a Bayesian setting, but propensity
score use in the Bayesian paradigm is not straightforward. Conditioning on a
propensity score in the outcome model in a Bayesian design does not correspond
to a valid use of Bayes theorem with correct posterior inference (Zigler,
2016).

We chose to borrow based on the similarity between regression coefficients and
parameters, and not on the causal effect itself. If we had access to the
individual-level causal effects, we could avoid the use of regression models
entirely and base our inference and borrowing on the distribution of
difference scores. Unfortunately, such data are not available as we can only
observe one potential outcome, and never both. Furthermore, individual-level
unobserved potential outcomes cannot be predicted without making assumptions
about unidentifiable parameters. Nevertheless, it may be possible to use a
model to estimate mean causal effects and base the borrowing on the mean
potential outcome difference scores. We leave this possibility open for future
consideration. We also developed model priors that discourage borrowing as the
number of predictors grows, but optimizing this prior to minimize potential
MSE increases requires additional investigation.

In summary, we proposed showed how to use two existing modelling strategies
and showed how to borrow from supplemental data sources using MEMs. The
estimators can even be used to borrow from observation data to estimate causal
effects in an RCT provided that the share predictors of the outcome. Both the
BLM and BART estimators performed well under simulation scenarios, and the
estimators gave sensible results in the application, dramatically reducing the
variance of the causal effect estimator. Throughout, we assume that we have a
primary data source and supplemental sources. We plan to extend our estimators
to cases with symmetrical borrowing, where no source is considered primary,
but we wish to borrow between sources as appropriate to estimate causal
effects.

## Acknowledgements

This research was partially funded by NIH under grants P30-CA077598,
R01-CA214824, and R01-CA225190 from the National Cancer Institute and
R03-DA041870, R01-DA046320 and U54-DA031659 from the National Institute on
Drug Abuse and FDA Center for Tobacco Products (CTP). The content is solely
the responsibility of the authors and does not necessarily represent the
official views of the NIH or FDA CTP.

## Simulation Code and Software

Simulation code is available at https://github.com/jeffrey-boatman/asymmetric-
borrowing-simulation. The `R` package `borrowr` implements the estimators
described in this manuscript and is available on the Comprehensive R Archive
Network. Research data are not shared.

## References

  * Boatman _and others_ (2019) Boatman, Jeffrey A., Vock, David M. and Koopmeiners, Joseph S. (2019).  borrowr: Estimate Causal Effects with Borrowing Between Data Sources.  R package version 0.1.1.
  * Chipman _and others_ (1998) Chipman, Hugh A., George, Edward I. and McCulloch, Robert E. (1998).  Bayesian CART model search.  Journal of the American Statistical Association 94(443), 935–948.
  * Chipman _and others_ (2010) Chipman, Hugh A., George, Edward I. and McCulloch, Robert E. (2010).  BART: Bayesian additive regression trees.  Annals of Applied Statistics 4(1), 266–298.
  * Denlinger _and others_ (2016) Denlinger, R. L., Smith, T. T., Murphy, S. E., Koopmeiners, J. S., Benowitz, N. L., Hatsukami, D. K., Pacek, L. R., Colino, C., Cwalina, S. N. and Donny, E. C. (2016).  Nicotine and anatabine exposure from very low nicotine content cigarettes.  Tobacco Regulatory Science 2(2), 186–203.
  * Donny _and others_ (2015) Donny, Eric C., Denlinger, Rachel L., Tidey, Jennifer W., Koopmeiners, Joseph S., Benowitz, Neal L., Vandrey _and others_. (2015).  Randomized Trial of Reduced-Nicotine Standards for Cigarettes.  New England Journal of Medicine 373(14), 1340–1349.
  * Dorie _and others_ (2019) Dorie, Vincent, Hill, Jennifer, Shalit, Uri, Scott, Marc and Cervone, Dan. (2019).  Automated versus Do-It-Yourself Methods for Causal Inference: Lessons Learned from a Data Analysis Competition.  Statist. Sci. 34(1), 43–68.
  * Gelman (2006) Gelman, Andrew. (2006).  Prior distributions for variance parameters in hierarchical models(Comment on Article by Browne and Draper).  Bayesian Analysis 1(3), 515–534.
  * Hatsukami _and others_ (2018) Hatsukami, Dorothy K., Luo, Xianghua, Jensen, Joni A., Al’Absi, Mustafa, Allen, Sharon S., Carmella, Steven G. _and others_. (2018, sep).  Effect of immediate vs gradual reduction in nicotine content of cigarettes on biomarkers of smoke exposure a randomized clinical trial.  JAMA - Journal of the American Medical Association 320(9), 880–891.
  * Hernández _and others_ (2018) Hernández, Belinda, Raftery, Adrian E., Pennington, Stephen R. and Parnell, Andrew C. (2018).  Bayesian Additive Regression Trees using Bayesian model averaging.  Statistics and Computing 28, 869–890.
  * Hill (2011) Hill, Jennifer L. (2011).  Bayesian Nonparametric Modeling for Causal Inference.  Journal of Computational and Graphical Statistics 20(1), 217–240.
  * Hobbs _and others_ (2011) Hobbs, Brian P., Carlin, Bradley P., Mandrekar, Sumithra J. and Sargent, Daniel J. (2011).  Hierarchical Commensurate and Power Prior Models for Adaptive Incorporation of Historical Information in Clinical Trials.  Biometrics 67, 1047–1056.
  * Ibrahim and Chen (2000) Ibrahim, Joseph G. and Chen, Ming-Hui. (2000).  Power prior distributions for regression models.  Statistical Science 15(1), 46–60.
  * Kaizer _and others_ (2018) Kaizer, Alexander M., Koopmeiners, Joseph S. and Hobbs, Brian P. (2018).  Bayesian hierarchical modeling based on multisource exchangeability.  Biostatistics 19(2), 169–184.
  * Lunceford and Davidian (2004) Lunceford, J. K. and Davidian, M. (2004).  Stratification and weighting via the propensity score in estimation of causal treatment effects: a comparative study.  Statistics in Medicine 23, 2937–2960.
  * Nardone _and others_ (2016) Nardone, Natalie, Donny, Eric C., Hatsukami, Dorothy K., Koopmeiners, Joseph S., Murphy, Sharon E., Strasser, Andrew A., Tidey, Jennifer W., Vandrey, Ryan and Benowitz, Neal L. (2016).  Estimations and predictors of non-compliance in switchers to reduced nicotine content cigarettes.  Addiction 111(12), 2208–2216.
  * Nethery _and others_ (2019) Nethery, Rachel C., Mealli, Fabrizia and Dominici, Francesca. (2019).  Estimating population average causal effects in the presence of non-overlap: The effect of natural gas compressor station exposure on cancer mortality.  Annals of Applied Statistics 13(2), 1242–1267.
  * Raftery _and others_ (1997) Raftery, Adrian E, Madigan, David and Hoeting, Jennifer A. (1997).  Bayesian Model Averaging for Linear Regression Models.  Journal of the American Statistical Association 92(437), 179–191.
  * Robins and Hernán (2009) Robins, James M. and Hernán, Miguel A.l. (2009).  Longidtudinal Data Analysis, Chapter 23: Estimation of the Causal Effects of Time-Varying Exposures. Boca Raton: Chapman & Hall/CRC.
  * Rubin (1973) Rubin, Donald B. (1973).  The Use of Matched Sampling and Regression Adjustment to Remove Bias in Observational Studies.  Biometrics 29, 185–203.
  * Rubin (1981) Rubin, Donald B. (1981).  The Bayesian Bootstrap.  The Annals of Statistics 9(1), 130–134.
  * Smith _and others_ (2019) Smith, Tracy T., Koopmeiners, Joseph S., Tessier, Katelyn M., Davis, Esa M., Conklin, Cynthia A., Denlinger-Apte, Rachel L. _and others_. (2019, oct).  Randomized Trial of Low-Nicotine Cigarettes and Transdermal Nicotine.  American Journal of Preventive Medicine 57(4), 515–524.
  * Tan and Roy (2019) Tan, Yaoyuan Vincent and Roy, Jason. (2019, aug).  Bayesian additive regression trees and the General BART model.  Statistics in Medicine 38, 5048–5069.
  * Tibshirani (1996) Tibshirani, Robert. (1996, jan).  Regression Shrinkage and Selection Via the Lasso.  Journal of the Royal Statistical Society: Series B (Methodological) 58(1), 267–288.
  * Wang _and others_ (2015) Wang, Chi, Dominici, Francesca, Parmigiani, Giovanni and Zigler, Corwin Matthew. (2015).  Accounting for uncertainty in confounder and effect modifier selection when estimating average causal effects in generalized linear models.  Biometrics 71, 654–665.
  * Zigler (2016) Zigler, Corwin Matthew. (2016).  The Central Role of Bayes’ Theorem for Joint Estimation of Causal Effects and Propensity Scores.  The American Statistician 70(1), 47–54.