# Multiple Domain Causal Networks

Tianhui Zhou1, William E. Carson IV2, Michael Hunter Klein3, David
Carlson1,3,4
1Department of Biostatistics and Bioinformatics
2Department of Biomedical Engineering
3Department of Electrical and Computer Engineering
4Department of Civil and Environmental Engineering
Duke University
Durham, NC 27708
{tianhui.zhou, wec14, michael.klein413<EMAIL_ADDRESS>
###### Abstract

Observational studies are regarded as economic alternatives to randomized
trials, often used in their stead to investigate and determine treatment
efficacy. Due to lack of sample size, observational studies commonly combine
data from multiple sources or different sites/centers. Despite the benefits of
an increased sample size, a naïve combination of multicenter data may result
in incongruities stemming from center-specific protocols for generating
cohorts or reactions towards treatments distinct to a given center, among
other things. These issues arise in a variety of other contexts, including
capturing a treatment effect related to an individual’s unique biological
characteristics. Existing methods for estimating heterogeneous treatment
effects have not adequately addressed the multicenter context, but rather
treat it simply as a means to obtain sufficient sample size. Additionally,
previous approaches to estimating treatment effects do not straightforwardly
generalize to the multicenter design, especially when required to provide
treatment insights for patients from a new, unobserved center. To address
these shortcomings, we propose Multiple Domain Causal Networks (MDCN), an
approach that simultaneously strengthens the information sharing between
similar centers while addressing the selection bias in treatment assignment
through learning of a new feature embedding. In empirical evaluations, MDCN is
consistently more accurate when estimating the heterogeneous treatment effect
in new centers compared to benchmarks that adjust solely based on treatment
imbalance or general center differences. Finally, we justify our approach by
providing theoretical analyses that demonstrate that MDCN improves on the
generalization bound of the new, unobserved target center.

##  1 Introduction

Recent advancements in deep learning have facilitated characterization and
modeling of complex relationships [36, 40]. Accordingly, deep learning has
since been applied to problems of personalized medicine, with the goal to
accurately quantify each individual’s heterogeneous response to a given
treatment [24]. Due to the high cost of conducting experiments, researchers
often turn to observational data to find clues through building predictive
models or making quantitative assessments [10]. To enlarge the overall sample
size, data collected from multiple different centers (clusters) are combined.
This practice is commonly referred to as a “multicenter observational study”
[14]. Naïvely combining data from multiple sources or ignoring the underlying
cluster structure can be problematic. For example, individuals from different
centers could have distinctive feature distributions with minimum overlap.
Moreover, the mechanisms underlying treatment selection bias may differ from
center to center. The combination of these issues present a challenge to
researchers attempting to quantitatively model the effects of treatment
assignment of patients from a new center. Limited work has been conducted for
the causal estimation in a multicenter setting [46], as the vast majority of
existing approaches address only part of the problem. To bridge this gap, our
goal is to design a framework that jointly addresses the structural
considerations of multicenter data and robust estimation of heterogeneous
treatment effects for individuals from new centers.

Different centers or different treatment groups within a same center could
mathematically be described as feature spaces $X\in\mathcal{X}$ of different
distributions. This aligns with the definition of “domain” in the domain
adaptation literature where the similarity among different domains are
leveraged to improve a model’s generalization [6]. Inspired by the idea, we
propose Multiple Domain Causal Networks (MDCN): a novel approach that takes
into account both selection bias in treatment assignment and discrepancies
inherent to data collected at different centers. Succinctly, we employ domain
adaptation to simultaneously increase the overlap between treatment groups and
match centers with more similarities to gain robustness in inference. Moving
forward, we will use “domain” to refer to a center or cluster from which data
was collected. Not only can the proposed approach be applied to data collected
from multiple medical centers, but it can also can be used in broader contexts
(e.g., considering individuals as separate domains). For instance, we motivate
our method using neural data collected from mice, where behavioral assays can
be framed as treatment regimes and we expect heterogeneous treatment effects
between individual mice. This approach is highly relevant to observational
studies that evaluate the impact of pharmaceuticals [12] or other
interventions based on brain measurements (e.g. neurostimulation [9]). Our
main contributions include the following:

  1. 1.

We formulate the heterogeneous treatment effect estimation problem under the
context of multicenter observational studies.

  2. 2.

We propose a new representation learning approach that accounts for both
domain-level discrepancies and selection biases in treatment assignment.

  3. 3.

We provide supporting theoretical proof demonstrating that the error on the
new, unobserved target center is bounded with the proposed method.

##  2 Related Work

To date, few works address out-of-domain treatment effect estimation,
especially on multicenter data. One example is generalized mixed effect models
that handle cluster-structured data derived from multiple domains [13, 34];
however, this approach is an additive linear model that requires specific
likelihood functions for its mean model and variance model, which restricts
customization and is less capable of handling heterogeneous relationships.
Instead, we will build on recent advances in using machine learning for causal
estimation and in domain adaption, briefly described below.

There are many recent advances in machine learning for causal inference. A
variety of methods have been used, including methods that match individuals
[3, 8, 38] or parts of the covariate space (e.g. tree-based methods [2, 22]. A
critical issue in causal inference is mismatch of the covariate space,
motivating weighting methods [1, 28]). Rather than weighting all samples, one
could disentangle the confounding effect from the feature space [26, 53], find
a more balanced embedding between treatment groups [41], or embed the
treatment assignment mechanism in a new latent representation [20, 43].
Despite demonstrable improvement in model generalization for causal inference,
these models do not account for multiple domain setups and their extensions to
out-of-domain prediction are not straightforward.

Domain adaptation methods provide elegant solutions to account for the
inherent distribution shift across different domains [4, 48, 49, 54]. Many
domain adaptation methods aim to find a latent variable representation where
the distributions of multiple domains are matched [55]. Matched latent domains
can be learned by penalizing with statistical measures such as the Maximum
Mean Discrepancy [7] or Wasserstein Distance [42], or by matching second order
statistics [47]. Another family of techniques make use of adversarial methods
to make the latent sample’s source domain identity indistinguishable [16, 23,
30]. It has been shown that too strong of a latent domain matching penalty can
harm performance, especially when source domains are unrelated to each other
[33] or suffer from label shift [31]. To account for this issue, several
techniques have been developed for relaxing the latent domain matching penalty
using other objective loss functions [23] or by explicitly re-weighting the
importance of matching domains through Multiple Domain Matching Networks
(MDMN) [29]. MDMN uses a Wasserstein distance approximation to dynamically re-
weight the importance of domain matching depending on the proximity of
domains, choosing to use only the most relevant source domains. While these
methods leverage advanced domain adaptation techniques, it is crucial in our
situation to also examine the within-domain differences in treatment
assignment.

In summary, existing work only partially solves the problem of predicting out-
of-domain treatment effects with data from multiple sources.

##  3 Problem Setup

Suppose we have data from $S$ domains (centers). Let $\\{1,\cdots,S-1\\}$
represent source domains, and let $S$ represent the unobserved target domain.
We assume a binary treatment condition with label $T\in\\{0,1\\}$. The feature
space and the potential outcome space are represented as
$X\in\mathcal{X}\subset\mathbb{R}^{p_{x}}$ and
$\\{Y(0),Y(1)\\}\in\mathcal{Y}\subset\mathbb{R}^{p_{y}}$, respectively. For
any domain $s$, the selection bias in observational studies causes $X$ to have
distinct distributions in the control group ($T=0$), treatment group ($T=1$),
and overall domain.

$\text{Control: }X|_{T=0}\sim D_{s,0};~{}\text{Treatment: }X|_{T=1}\sim
D_{s,1};~{}\text{Overall: }X\sim D_{s}.$

We use $D_{s,0},D_{s,1}$, and $D_{s}$ to represent their distribution
functions, and their relationship can be described by
$(1-p_{s}^{T=1})D_{s,0}+p_{s}^{T=1}D_{s,1}=D_{s}$, with $p_{s}^{T=1}$
indicating the marginal probability of treatment assignment. Likewise,
$\\{g_{s,0},g_{s,1}\\}$ represents the unknown outcome functions for
$\\{Y(0),Y(1)\\}$.

The observable information of an individual data sample of any source domain
$s$ is represented via a tuple, $\\{T,X,Y=Y(T)\\}$. Given a target domain $S$
with $X\sim D_{S}$, our goal is to accurately infer both potential outcomes,
so that we can more accurately estimate treatment outcomes and therefore
assign individuals to the treatments that are likely to benefit them most.

There are two obstacles that are likely to hurt generalization to the target
domain. The first is selection bias, since treatment assignment is often not
randomized in observational studies. As a result, the acquired samples on
which the model is trained cannot objectively reflect the error on the full
domain, as $D_{s,0}\neq D_{s,1}\neq D_{s},\forall s$. The second consideration
is domain-level discrepancies, which may inhibit the effects of reducing the
error over the combined source domains. More formally, for different domains
$a$ and $b$, we could describe the discrepancies as the shifts in feature
space distributions and underlying outcome functions, represented as,
$D_{a}\neq D_{b}$ and $\\{g_{a,0},g_{a,1}\\}\neq\\{g_{b,0},g_{b,1}\\}$,
respectively. Here, we make similar assumptions to those made in the causal
literature [21, 37]; however, we make modifications that account for the
extension to our multi-domain scenario.

######  Assumption 1 (Within-domain positivity).

For domain $s$ where $X\sim D_{s}$, the probability of assignment to any
treatment group is bounded away from zero: $0<Pr(T=1|X,s)<1$.

######  Assumption 2 (Within-domain consistency).

For domain $s$, the observed outcome with assigned treatment is equal to its
potential outcome: $Y|T,X\sim D_{s}=Y(T)|T,X\sim D_{s}$.

######  Assumption 3 (Within-domain ignorability).

For domain $s$, the potential outcomes are jointly independent of the
treatment assignment conditional on $X$: $[Y(0),Y(1)]\perp T|X\sim D_{s}$.

Considering treatment assignment as probabilistic ensures that studying the
treatment difference forms a meaningful target for each domain. The original
assumption does not involve $s$. This could create a problem where the
treatment assignment is probabilistic overall ($0<Pr(T=1|X)<1$) but
deterministic at domain level ($Pr(T=1|X,s)=1$ or $Pr(T=1|X,s)=0$).
Assumptions 2 and 3 are now also stated with respects to each domain to better
account for differences at the domain level. As a result, individuals with
identical features but from different domains may have treatment assignments
and responses that are domain-dependent and thus different overall.

##  4 Multiple Domain Causal Networks

Multiple Domain Causal Networks (MDCN) is a novel framework that both
addresses selection bias in treatment assignment and leverages similarity from
different centers. MCDN has three main components, all implemented as neural
networks: a feature embedding network
$\phi:~{}\mathbb{R}^{p_{x}}\rightarrow\mathbb{R}^{q}$ that learns new
representations, and two outcome networks
$h_{0}:~{}\mathbb{R}^{q}\rightarrow\mathbb{R}$ and
$h_{1}:~{}\mathbb{R}^{q}\rightarrow\mathbb{R}$ that infer potential outcomes.
The roles of these three components are depicted in Figure 1. Using the
taxonomy in Künzel et al. [27], our outcome networks can be considered as a
“T-learner”.

First, the original feature space $X$ is embedded in $\phi(X)$. This embedding
is incorporated into three loss terms, the purpose of each being to endow it
with specific properties. Through $L(\phi,h_{0},h_{1})$, $\phi$ is injected
with information predictive of the outcomes. $L_{BT}(\phi)$ and $L_{CD}(\phi)$
are two regularization terms. The former makes $\phi$ more robust against the
selection bias between treatment groups, and the latter reduces the cross-
domain differences so that generalization is improved by learning from
examples from similar domains.

$L_{all}(\phi,h_{0},h_{1})=L(\phi,h_{0},h_{1})+\alpha L_{BT}(\phi)+\beta
L_{CD}(\phi)$ (1)

Figure 1: MDCN, $\phi,h_{0},h_{1}$ are optimized together with respective
losses.

The full loss is summarized in (1), with $\alpha$ and $\beta$ as tuning
parameters. In all our experiments, we fix $\alpha=$5e-4 and $\beta$=1e-3, as
chosen from the range of recommended values in Shalit et al. [41] and Li et
al. [29]. Our model has empirically robust performance with these values,
although they could be tuned in practice. We give the form of
$L(\phi,h_{0},h_{1})$ below. The design of $L_{BT}$ and $L_{CD}$ are explained
in Sections 4.1 and 4.2. $L(\phi,h_{0},h_{1})$ trains the potential outcome
models $h_{0}$ and $h_{1}$ to predict the observed outcomes given their
assigned treatments. We use $l(\cdot,\cdot)$ to represent an arbitrary loss
function, the form of which can be chosen based on the requirements of the
modeling problem (e.g., cross-entropy loss for binary or categorical outcomes,
squared loss for continuous outcomes). (2) represents the expectation of this
loss as a summation over all source domains,

$L(\phi,h_{0},h_{1})=\sum\nolimits_{s=1}^{S-1}\biggl{[}E_{x\sim
D_{s,0},y,t}l(h_{0}(\phi(x),y(t=0))+E_{x\sim
D_{s,1},y,t}l(h_{1}(\phi(x),y(t=1))\biggl{]}$ (2)

We note that $h_{0}$ and $h_{1}$ are used to predict the potential outcomes
universally across all domains. This strategy is more scalable compared to the
alternative of learning a unique function for each domain. Additionally, we
believe that when given enough capacity, $h_{0}$ and $h_{1}$ can still discern
the differences at the domain level when combined with the extra properties
endowed in $\phi$ learned via the losses $L_{BT}$ and $L_{CD}$. However, when
we have a small number of domains, we could construct domain-specific
$h_{0,s}$ and $h_{1,s}$. The prediction on the target domain could be based on
a weighted summation of the outcome functions from its neighboring sources
domains. We also view it as a future goal of varying the design of the
potential outcome functions. In Appendix B.1, Algorithm 1 describes the
pseudo-code for the full implementation.

###  4.1 Between-Treatment Adjustment

Our approach to between-treatment adjustment is motivated by counterfactual
regression (CFR) [41]. An embedding $\phi(x)$ with increased overlap between
different treatment groups is learned to reduce the effects of selection bias,
as supported by the theory of Ben-David et al. [6] that generalization between
more similar spaces have a lower error bound. To measure and encourage this
overlap, we penalize with the Kantorovich-Rubinstein form of the Wasserstein-1
distance [51].

$\sup_{\|{f_{bt}}\|_{L}\leq
1}\mathbb{E}_{x\sim{D}_{0}}[f_{bt}(\phi(x))]-\mathbb{E}_{x\sim{D}_{1}}[f_{bt}(\phi(x))]$
(3)

$D_{0}$ and $D_{1}$ in (3) characterize how $X$ is distributed on the
treatment group and control group. $f_{bt}$ is drawn from the family of
Lipschitz functions with constant 1, which can be approximated by neural
networks with a gradient penalty [19]. Through this objective, $\phi(x)$ will
be encouraged to reduce the distance between-treatment groups. However, this
only resolves the imbalance within a single domain. The direct application of
CFR to balance all samples does not address the difficulties inherent to data
collected from multiple domains. While overall balance is encouraged and
enforced, the balancing of specific domains may be neglected. Moreover, the
relative closeness of domains in the original space could also become
distorted with this crude adjustment. As a modification, we adopt a separate
${f_{bt}^{s}}$ for each domain. This aims at reducing the distance between-
treatment groups for all domains. As ${f_{bt}^{s}}$ is domain-specific, the
resulting adjustment is therefore also domain-specific. Additionally, there is
less impact on cross-domain relationships. In Section 6.1, we provide
visualizations to showcase the advantage of this approach over CFR. With this
new strategy, the between group loss is,

$L_{BT}(\phi,f_{bt})=\sum\nolimits_{s=1}^{S-1}\biggl{[}\sup_{\|{f_{bt}^{s}}\|_{L}\leq
1}\mathbb{E}_{x\sim{D}_{i,0}}[f_{bt}^{s}(\phi(x))]-\mathbb{E}_{x\sim{D}_{i,1}}[f_{bt}^{s}(\phi(x))]\biggl{]}.$
(4)

We represent $L_{BT}(\phi,f_{cd})$ with notation $L_{BT}(\phi)$, as $f_{bt}$
is an auxiliary function to learn distances. We combine
$f_{bt}=\\{f_{bt}^{1},\cdots,f_{bt}^{S-1}\\}$ using a single neural network
with $S-1$ outputs, which saves computational cost and encourages information
sharing between domains, and use the methods of Gulrajani et al. [19] for
learning.

###  4.2 Cross-Domain Adjustment

Figure 2: Depiction of cross-domain adjustment with three domains. Wider
arrows indicate larger weights/stronger relationships for similar domains.

We implement a cross-domain adjustment to encourage similarities at the domain
level. This is achieved by extending the Multiple Domain Matching Network
(MDMN) [29]. Akin to the between-treatment adjustment, we measure and reduce
the distance based on the Wasserstein-1 distance. Instead of enforcing the
distance penalty statically between any two domains, MDMN adopts a weighting
scheme which assigns larger weight to more similar domains in (5). For domains
that are far apart or dissimilar, forcefully reducing their differences may
result in excessive loss of predictive information. As part of the predictive
information is shared by all domains, the rest may be specific to only a few.
The domain-specific information is reflected by these differences, which
cannot be covered by the domain-invariant representation alone but is helpful
in capturing the shift in outcome functions. Therefore, smaller weights are
allocated in these instances. Ideally, each domain is matched to a few highly-
correlated domains in space and forms clusters. For any domain $s$, this is
formulated as,

$d(D_{s},D_{/s})=\sup_{\|{f_{cd}^{s}}\|_{L}\leq
1}\mathbb{E}_{x\sim{D}_{s}}[f_{cd}^{s}(\phi(x))]-\sum\limits_{i\neq
s}\mathbb{E}_{x\sim{D}_{i}}[w_{s,i}f_{cd}^{s}(\phi(x))].$ (5)

Specifically, the weight $w_{s,i}$ is calculated according to the pairwise
differences between any two domains in (6), and $\sum_{i\neq s}w_{s,i}=1$.
While this distance is minimized by picking the single closest domain, we
encourage robustness by smoothing over nearby domains. The weighted summation
of all other domains can be seen as a pseudo domain, $D_{/s}$, with the well-
defined distribution, $D_{/s}=\sum\nolimits_{i\neq s}w_{s,i}D_{i}$:

$\begin{split}&\boldsymbol{l}_{s}=\\{l_{s,i}\\}_{i\neq
s},~{}l_{s,i}=\mathbb{E}_{x\sim{D}_{s}}[f_{cd}^{s}(\phi(x))]-\mathbb{E}_{x\sim{D}_{i}}[f_{cd}^{s}(\phi(x))],\\\
&\boldsymbol{w}_{s}=\\{w_{s,i}\\}_{i\neq
s}=\text{softmax}(-\boldsymbol{l}_{s}).\end{split}$ (6)

In practice, $\boldsymbol{l}_{s}$ is often accompanied with a temperature term
to increase its stability. $\boldsymbol{l}_{s}$ at $t^{\prime}$th iteration
could be represented as:
$\boldsymbol{l}_{s}^{t}=.9\boldsymbol{l}_{s}^{t-1}+.1\boldsymbol{l}_{s}^{c}$,
with $\boldsymbol{l}_{s}^{c}$ being the estimates from the current batch [29].
Aggregating over all domains gives us the cross domain adjustment term:

$L_{CD}(\phi,f_{cd})=\sum\limits_{s=1}^{S}\bigg{[}\sup_{\|{f_{cd}^{s}}\|_{L}\leq
1}\mathbb{E}_{x\sim{D}_{s}}[f_{cd}^{s}(\phi(x))]-\sum\limits_{i\neq
s}\mathbb{E}_{x\sim{D}_{i}}[w_{s,i}f_{cd}^{s}(\phi(x))]\bigg{]}.$ (7)

Likewise, we can also use notation $L_{CD}(\phi)$ for $L_{CD}(\phi,f_{cd})$.
To control for the computational cost, we implement
$f_{cd}=\\{f_{cd}^{1},\cdots,f_{cd}^{S}\\}$ in the form of a multi-output
function. However, due to the pairwise comparisons in weight calculation, the
complexity $\mathcal{O}(S^{2})$ is needed as $S$ increases. We train this loss
and between-treatment adjustment simultaneously on $\phi$ again using the
gradient penalties of Gulrajani et al. [19], which is depicted in Figure 2.

##  5 Theoretical Analysis

Here, we provide theoretical analyses to explain the design of MDCN. The
overarching goal is to bound the error on the unlabeled target domain via
labeled source domain data. Full proofs for the theoretical results presented
in this section can be found in Appendix A. We give a few preliminaries below
before presenting the final bound. To quantify the difference between
hypothesis (or outcome) functions, we state the probabilistic discrepancy
defined in Li et al. [29]:

######  Definition 1 (Probabilistic discrepancy).

For two hypotheses $h$ and $h^{\prime}:\mathbb{R}^{P}\rightarrow\mathbb{R}$,
their difference given a probabilistic distribution $D$ over $\mathcal{X}$ is
defined as $\gamma(h,h^{\prime}|D)=\mathbb{E}_{x\sim D}|h(x)-h^{\prime}(x)|$.

Next, we define a family of our hypothesis functions in Definition 2. We limit
our proposed hypothesis class to the Lipschitz family with parameter
$\lambda$. $\\{h_{0},h_{1}\\}\in F_{\lambda}$.

######  Definition 2 (Lipschitz continuity).

A function $f:\mathbb{R}^{P}\rightarrow\mathbb{R}$ is Lipschitz continuous
with parameter $\lambda$ if
$|f(x_{1})-f(x_{2})|\leq\lambda||x_{1}-x_{2}||_{2}$ holds for any vectors
$x_{1},x_{2}\in\mathcal{X}$. We denote the family as $\mathcal{F}_{\lambda}$.

We assume that the true outcome functions are also included in this family but
with a different smoothness parameter $\lambda^{*}$,
$\\{g_{s,0},g_{s,1}\\}_{s=1}^{S}\in F_{\lambda^{*}}$. This encourages smooth
transitions in the outcome labels with regards to changes in the feature
space. Though smoothness may not hold in practice, Lipschitz continuous
functions provide a reasonable approximation to a wide range of functions
[44]. With all the preliminaries, we can present the overall bound on the
target domain in Theorem 1,

######  Theorem 1.

For any positive weights $w=\\{w_{s}\\}_{s=1}^{S-1}$ with
$\sum_{s=1}^{S-1}w_{s}=1$, the discrepancy between the true hypothesis
functions $\\{g_{0,S},g_{1,S}\\}$ and the proposed hypothesis functions
$\\{h_{0},h_{1}\\}$ on target domain $S$ is bounded by,

$\begin{aligned}
\gamma(h_{0},g_{S,0}|D_{S})+\gamma(h_{1},g_{S,1}|D_{S})&\leq(\lambda+\lambda^{*})[2W_{1}(D_{S},\sum\nolimits_{s=1}^{S-1}w_{s}D_{s})+\sum\nolimits_{s=1}^{S-1}w_{s}W_{1}(D_{s,0},D_{s,1})]\\\
&+\sum\nolimits_{s=1}^{S-1}w_{s}[\gamma(h_{0},g_{s,0}|D_{s,0})+\gamma(h_{1},g_{i,1}|D_{s,1})]+\gamma_{0}^{*}+\gamma_{1}^{*}.\end{aligned}$

(8)

$[\gamma(h_{0},g_{s,0}|D_{s,0})+\gamma(h_{1},g_{i,1}|D_{s,1})]$ represents the
probabilistic discrepancy between $\\{h_{0},h_{1}\\}$ and the true outcome
functions $\\{g_{s,0},g_{s,1}\\}$ on the observable part of the data:
$D_{s,0},D_{s,1}$. Minimizing this corresponds to $L(\phi,h_{0},h_{1})$ in
(2). In practice, we could modify the discrepancy mildly such as by optimizing
the squared loss in regressions instead for more stability.
$W_{1}(D_{s,0},D_{s,1})$ measures the distance between the treatment space and
the control space in domain $s$, and it matches with the structure of
$L_{BC}(\phi)$. Likewise, $L_{CD}(\phi)$ is inspired by
$W_{1}(D_{S},\sum\nolimits_{s=1}^{S-1}w_{s}D_{s})$. Though the bound is
dependent on the values of $\lambda$ or $\lambda^{*}$, their explicit values
are not needed in practice, as they can be incorporated into tuning parameters
$\alpha$ and $\beta$. Lastly, $\gamma_{0}^{*}+\gamma_{1}^{*}$ depicts the
fundamental difference in the true outcome functions
$\\{g_{s,0},g_{s,1}\\}_{s=1}^{S}$ across domains, and cannot be optimized. A
large $\gamma_{0}^{*}+\gamma_{1}^{*}$ signals the existence of a huge shift in
outcome functions. In that case, we are less guaranteed to achieve accurate
predictions on the target domain. The weight $w=\\{w_{s}\\}_{s=1}^{S-1}$ also
plays a key role in this bound. In (8), increased emphasis is placed on source
domains that are close to the target domain. Samples from these similar source
domains may serve as good reference cases and provide valuable insights in
improving the performance on the target domain.

##  6 Experiments

We use two examples to demonstrate how MDCN improves estimation of
heterogeneous treatment effects on the target domain. Efficacy is evaluated
according to estimation of the conditional average treatment effect (CATE)
[41]. We modify CATE by conditioning the quality on the domain label $S$ to
take into consideration the shift in outcome functions:

$\tau(x,S)=\mathbb{E}[Y(1)|x,S]-\mathbb{E}[Y(0)|x,S].$ (9)

We use the precision in estimation of heterogeneous treatment effects (PEHE)
from Hill [22] to measure the distance between CATE and its estimates:
$\text{PEHE(S)}=\sqrt{\mathbb{E}_{x\sim
D_{s}}[\tau(x,S)-\hat{\tau}(x,S)]^{2}}$. As evaluating CATE requires access to
ground truth values, we create synthetic and semi-synthetic data for this
purpose. Causal forest (CF) is included as a benchmark due to its previous
application to multicenter observational study data [2, 46]. Other benchmarks
include different variants of MDCN, which also serves as an ablation study
that helps us understand the benefits of different components of MDCN. Apart
from CF, we note that all variants have identical architectures for the
outcome model and feature embedding model to ensure a fair comparison. CFR is
included as a naïve benchmark that only addresses treatment group imbalance
globally [41]. The proposed between-treatment adjustment $L_{BT}(\phi)$ in (4)
upgrades CFR by performing the novel adjustment within each domain, which we
call the domain CFR (DCFR). Performance of MDMN is also compared to illustrate
the inadequacy of applying the domain-level adjustment alone. Lastly, we
combine MDMN and CFR to form MDMNCFR. Though it adjusts at both cross-domain
and between-treatment levels, its design is inferior to MDCN which is
rigorously backed by supporting theories. Additional details can be found in
Appendix B.1.

###  6.1 Synthetic Circular Data

In this section, we create a synthetic case with shifts in both feature spaces
and outcome functions across domains. We define the dimension of the feature
space $\mathcal{X}\in\mathbb{R}^{3}$ and the number of domains $S=10$. These
domains are labeled by $s\in\\{0,\cdots,9\\}$. Any domain $s$ versus the rest
of the combined domains can be regarded as the setup for source versus target
domains. To simulate the data, we first use $s$ to create an angle parameter:
$\angle_{s}=s\times\pi/10$, which evenly splits a circle. Through
$\angle_{s}$, we vary domains and outcome functions with the following
procedures.

Domain shift: $\displaystyle
X|s=\\{x_{1},x_{2},x_{3}\\}\sim\mathcal{N}([4\sin(\angle_{s}),4\cos(\angle_{s}),0]^{\prime},\mathbb{I}_{3});$
Potential outcomes: $\displaystyle
Y(0)|X,s\sim\mathcal{N}(1.5[\sin(x_{1}+\angle_{s})+\cos(x_{2}+x_{3}+\angle_{s})],1),$
$\displaystyle
Y(1)|X,s\sim\mathcal{N}(1.5[\cos(x_{1}+\angle_{s})+\sin(x_{2}+x_{3}+\angle_{s})],1);$
Treatment Assignment :
$\displaystyle\mathbb{P}(T=1|X)=1/[1+\exp(-0.5x_{1}-0.5x_{2}-2x_{3})].$

For source domains, the observed outcomes correspond to the outcomes with
assigned treatments $Y=Y(0)(1-T)+Y(1)T$. We simulate 2,000 samples for each
domain and repeat this procedure 10 times for variability assessment. The
distributions for 10 domains in their first two dimensions are visualized in
Figure 3(a). The location shift of 10 domains forms a circle, and each domain
$s$ has 2 close neighbors on its two sides with labels: $s-1\pmod{10}$ and
$s+1\pmod{10}$.

(a) Visualization of 10 Domains

(b) CFR Domain Similarity

(c) MDMN Domain Similarity

Figure 3:  3(a) visualizes the locations of the 10 domains. 3(b) and 3(c) are
the heat maps measuring the similarity between domains. For MDMN-based
approaches (MDCN included), we use the learned weight
$(w_{s^{\prime},s}+w_{s,s^{\prime}})/2$ to fill in the
grid${}_{s,s^{\prime}}$. For other methods, we use $\phi(x)$’s 1d t-SNE [32]
to estimate the distance $\boldsymbol{l_{s}}$ and weight $\boldsymbol{w_{s}}$
in (6) to plot the similarity (details in Appendix C). Darker colors indicate
more similarity. Ideally, the whole diagonal and two off-diagonal corners
should be dark with more decaying for grids further away from them. CFR
strongly connects some distant domains, such as 1 and 7 or 2 and 6.
Essentially, it reduces the distances in the diagonal to reduce the treatment
group imbalance marginally, which contradicts with how similarity is shared
between near domains. MDCN fully recovers the true similarity. Other methods
are visualized in Figure S1. Table 1: PEHE on Circular Data. MDMN is the most
competitive of all benchmarks. Adjusting treatment group imbalance within each
domain (DCFR) consistently outperforms CFR in all domains. Wilcoxon signed-
rank test is used to report the significance between our newly proposed
approaches and CFR. $P$ value $<5\cdot 10^{-2}$ ($\dagger$), $P$ value <
$1\cdot 10^{-4}$ ($\dagger\dagger$).

Target Domain | CF | MLP | CFR | MDMN | MDMNCFR†† | DCFR†† | MDCN††
---|---|---|---|---|---|---|---
0 | 1.24 $\pm$ .02 | 1.69 $\pm$ .08 | 1.59 $\pm$ .07 | 1.31 $\pm$ .05 | 1.33 $\pm$ .05 | 1.37 $\pm$ .07 | 1.28 $\pm$ .07
1 | 2.55 $\pm$ .03 | 2.01 $\pm$ .05 | 1.99 $\pm$ .08 | 1.42 $\pm$ .04 | 1.34 $\pm$ .04 | 1.48 $\pm$ .09 | 1.37 $\pm$ .04
2 | 1.43 $\pm$ .03 | 1.33 $\pm$ .06 | 1.32 $\pm$ .03 | 1.16 $\pm$ .04 | 1.17 $\pm$ .05 | 1.19 $\pm$ .03 | 1.18 $\pm$ .04
3 | 1.44 $\pm$ .03 | 1.27 $\pm$ .04 | 1.25 $\pm$ .03 | 1.01 $\pm$ .04 | 1.00 $\pm$ .03 | 1.02 $\pm$ .02 | 0.99 $\pm$ .03
4 | 1.45 $\pm$ .01 | 1.31 $\pm$ .06 | 1.29 $\pm$ .05 | 1.20 $\pm$ .02 | 1.23 $\pm$ .03 | 1.19 $\pm$ .04 | 1.16 $\pm$ .03
5 | 1.52 $\pm$ .02 | 1.34 $\pm$ .07 | 1.25 $\pm$ .05 | 1.28 $\pm$ .05 | 1.31 $\pm$ .08 | 1.25 $\pm$ .06 | 1.24 $\pm$ .07
6 | 1.41 $\pm$ .03 | 1.47 $\pm$ .07 | 1.32 $\pm$ .06 | 1.34 $\pm$ .07 | 1.26 $\pm$ .02 | 1.25 $\pm$ .03 | 1.32 $\pm$ .05
7 | 2.20 $\pm$ .02 | 1.82 $\pm$ .08 | 1.87 $\pm$ .12 | 1.64 $\pm$ .10 | 1.80 $\pm$ .07 | 1.49 $\pm$ .07 | 1.60 $\pm$ .06
8 | 1.69 $\pm$ .03 | 2.15 $\pm$ .05 | 2.10 $\pm$ .07 | 1.69 $\pm$ .04 | 1.71 $\pm$ .04 | 1.65 $\pm$ .07 | 1.62 $\pm$ .04
9 | 1.49 $\pm$ .03 | 2.17 $\pm$ .07 | 2.31 $\pm$ .08 | 1.65 $\pm$ .06 | 1.59 $\pm$ .08 | 1.61 $\pm$ .09 | 1.49 $\pm$ .06
Overall | 1.64 $\pm$ .04 | 1.66 $\pm$ .04 | 1.63 $\pm$ .04 | 1.37 $\pm$ .03 | 1.37 $\pm$ .03 | 1.35 $\pm$ .03 | 1.33 $\pm$ .03

The results on PEHE are summarized in Table 1. For lower dimensional
manifolds, tree-based approaches are advantageous as samples can easily
saturate the space. However, CF lacks the structure to address the shifts in
domains and outcome functions, so it only slightly outperforms MLP. Compared
with MLP, other methods with adjustments demonstrate clear improvements in
accuracy. DCFR minimizes the treatment imbalance within each domain instead of
marginally like CFR. This lessens the impact of domain-level relationships,
thus resulting in better performance. This result is visualized in Figure
3(b), as CFR distorts some of the cross-domain similarities to alleviate the
marginal treatment imbalance (e.g., regarding domain 2 and 6 to be closely
connected). DCFR visualized in S1 better maintains the original domain-level
connectivity. MDCM in Figure 4(c) can clearly recover the true similarity. The
combination of CFR and MDMN (MDNMCFR) does not improve MDMN. MDCN has a clear
advantage over all other benchmarks, achieving top marks in 5 out of 10
domains as well as best overall performance.

###  6.2 Semi-Synthetic Mouse Data

In this experiment, we make use of a mouse local field potential (LFP) dataset
collected in Gallagher et al. [15]. These data were collected from 21 mice in
accordance with guidelines provided by the Institutional Animal Care and Use
Committee (IACUC). LFPs were recorded while each mouse was exposed to two
conditions: “home cage” (resting state, low stress) and “open field” (mild
increased stress) [5]. We use these two states to mimic the control (home
cage) and treatment (open field) groups, thus providing an analog to
multicenter observational study data. Individual mice, similarly to individual
sites, are expected to demonstrate heterogeneity. We define an observation as
each 1 second-long window of brain signal recordings which is associated with
a continuous outcome, such as a behavioral outcome. Each mouse has less than
500 observations for each state. This treatment structure mimics approaches in
closed-loop brain stimulation, which chooses treatments based on current
neural representations, and is currently in use in depression [39] and seizure
prevention [45].

Table 2: PEHE on Mice Data. MDCN gives the overall best performance, and is
the best approach in 9 out of 21 domains. $P$ value $<5\cdot 10^{-2}$
($\dagger$), $P$ value < $1\cdot 10^{-4}$ ($\dagger\dagger$).

Target Domain | CF | MLP | CFR | MDMN | MDMNCFR† | DCFR†† | MDCN††
---|---|---|---|---|---|---|---
0 | 7.64 $\pm$ .07 | 1.61 $\pm$ .02 | 1.50 $\pm$ .02 | 1.52 $\pm$ .05 | 1.56 $\pm$ .04 | 1.45 $\pm$ .06 | 1.52 $\pm$ .06
1 | 3.73 $\pm$ .06 | 4.21 $\pm$ .12 | 4.12 $\pm$ .12 | 3.31 $\pm$ .19 | 3.35 $\pm$ .09 | 3.24 $\pm$ .12 | 3.45 $\pm$ .08
2 | 5.17 $\pm$ .04 | 5.42 $\pm$ .11 | 4.94 $\pm$ .13 | 3.29 $\pm$ .16 | 3.95 $\pm$ .26 | 3.70 $\pm$ .20 | 3.15 $\pm$ .23
3 | 6.92 $\pm$ .04 | 2.07 $\pm$ .05 | 2.11 $\pm$ .05 | 2.01 $\pm$ .07 | 2.16 $\pm$ .05 | 2.10 $\pm$ .06 | 2.17 $\pm$ .07
4 | 4.09 $\pm$ .12 | 2.61 $\pm$ .07 | 2.47 $\pm$ .06 | 2.57 $\pm$ .08 | 2.68 $\pm$ .10 | 2.65 $\pm$ .08 | 2.63 $\pm$ .07
5 | 6.68 $\pm$ .04 | 1.49 $\pm$ .12 | 1.37 $\pm$ .08 | 1.39 $\pm$ .03 | 1.34 $\pm$ .05 | 1.31 $\pm$ .04 | 1.40 $\pm$ .07
6 | 3.77 $\pm$ .05 | 2.19 $\pm$ .06 | 2.40 $\pm$ .14 | 2.26 $\pm$ .06 | 2.26 $\pm$ .06 | 2.33 $\pm$ .08 | 2.26 $\pm$ .08
7 | 6.06 $\pm$ .21 | 4.85 $\pm$ .08 | 5.11 $\pm$ .07 | 5.20 $\pm$ .08 | 5.22 $\pm$ .10 | 5.27 $\pm$ .09 | 5.11 $\pm$ .05
8 | 2.84 $\pm$ .05 | 2.27 $\pm$ .16 | 2.23 $\pm$ .10 | 1.71 $\pm$ .11 | 1.87 $\pm$ .17 | 1.73 $\pm$ .10 | 1.70 $\pm$ .09
9 | 2.53 $\pm$ .11 | 1.88 $\pm$ .16 | 1.76 $\pm$ .15 | 2.04 $\pm$ .22 | 1.57 $\pm$ .14 | 1.48 $\pm$ .14 | 1.31 $\pm$ .07
10 | 4.54 $\pm$ .25 | 1.72 $\pm$ .10 | 1.50 $\pm$ .08 | 1.46 $\pm$ .05 | 1.50 $\pm$ .03 | 1.40 $\pm$ .06 | 1.41 $\pm$ .03
11 | 7.54 $\pm$ .16 | 1.66 $\pm$ .03 | 1.72 $\pm$ .08 | 1.44 $\pm$ .07 | 1.41 $\pm$ .05 | 1.38 $\pm$ .06 | 1.24 $\pm$ .04
12 | 2.70 $\pm$ .14 | 2.81 $\pm$ .13 | 2.79 $\pm$ .08 | 2.90 $\pm$ .11 | 2.99 $\pm$ .08 | 2.87 $\pm$ .09 | 2.96 $\pm$ .10
13 | 2.60 $\pm$ .06 | 1.67 $\pm$ .19 | 1.56 $\pm$ .12 | 1.98 $\pm$ .12 | 2.26 $\pm$ .18 | 1.91 $\pm$ .15 | 2.23 $\pm$ .24
14 | 4.40 $\pm$ .42 | 1.48 $\pm$ .02 | 1.38 $\pm$ .02 | 1.39 $\pm$ .02 | 1.47 $\pm$ .02 | 1.43 $\pm$ .03 | 1.43 $\pm$ .03
15 | 2.49 $\pm$ .15 | 3.20 $\pm$ .15 | 2.95 $\pm$ .13 | 2.38 $\pm$ .11 | 2.37 $\pm$ .10 | 2.47 $\pm$ .14 | 2.54 $\pm$ .11
16 | 3.62 $\pm$ .17 | 2.66 $\pm$ .03 | 2.51 $\pm$ .13 | 2.60 $\pm$ .17 | 2.34 $\pm$ .09 | 2.46 $\pm$ .12 | 2.21 $\pm$ .09
17 | 2.66 $\pm$ .44 | 1.68 $\pm$ .03 | 1.57 $\pm$ .02 | 1.59 $\pm$ .05 | 1.51 $\pm$ .02 | 1.55 $\pm$ .05 | 1.51 $\pm$ .03
18 | 5.84 $\pm$ .17 | 6.14 $\pm$ .20 | 5.95 $\pm$ .16 | 4.90 $\pm$ .34 | 5.12 $\pm$ .23 | 5.17 $\pm$ .38 | 4.75 $\pm$ .27
19 | 4.42 $\pm$ .24 | 2.03 $\pm$ .07 | 1.81 $\pm$ .07 | 1.41 $\pm$ .05 | 1.41 $\pm$ .05 | 1.38 $\pm$ .03 | 1.37 $\pm$ .05
20 | 2.29 $\pm$ .15 | 2.49 $\pm$ .09 | 2.36 $\pm$ .12 | 2.38 $\pm$ .10 | 2.28 $\pm$ .09 | 2.29 $\pm$ .08 | 2.28 $\pm$ .09
overall | 4.41 $\pm$ .13 | 2.77 $\pm$ .09 | 2.58 $\pm$ .09 | 2.37 $\pm$ .08 | 2.41 $\pm$ .08 | 2.36 $\pm$ .08 | 2.32 $\pm$ .07

(a) DCFR

(b) MDMN

(c) MDCN

Figure 4:  T-SNE plot of $\phi(x)$. The domain 9 is chosen as the target
domain and indicated with an orange "T". For convenience, the three nearest
neighbors of the target domain are labeled as “0”, “1”, and “2”. The weight
matrix in MDCN and MDMN both identify the same three neighbors. Solid
transparent shades represent the control and treatment groups, respectively.
In 4(a), the target domain lies on the outer perimeter of these three
neighbors indicating lesser similarity. In 4(b), there exists a visible
separation between the treatment and the control groups in the neighbor “0”.
In 4(c), MDCN balances the treatment and control groups with respects to its
neighbors, and the target domain is properly surrounded by the three
neighbors. Other visualizations are provided in Appendix C.

Features meant to capture neural dynamics are generated from LFP data,
including power spectral density [52], coherence, and Granger causality [18]
for a total of with 9,856 features. To artificially increase sample size, we
employ an auto-encoder method that simultaneously interpolates the space
between samples and reduces dimensionality [35]. Further experimental details
are provided in Appendix B.2. Following preprocessing, each domain contains
20,000 samples with the feature space dimension reduced to a dimensionality
15. The new feature space still retains domain-specific information, as mouse
identity can be accurately predicted ($>90\%$). Next, we simulate the
potential outcomes. We randomly draw 15 correlated neural networks with the
input and output size: $\mathbb{R}^{15}\rightarrow\mathbb{R}$, for both
states, $\\{f_{1}^{0},\cdots,f_{15}^{0}\\}$,
$\\{f_{1}^{1},\cdots,f_{15}^{1}\\}$, that are designed to vary with the center
of each domain $s$,
$\bar{x}^{s}=\\{\bar{x}^{s}_{1},\cdots,\bar{x}^{s}_{15}\\}$ (details provided
in Appendix B.2). Essentially, domains that are close in center tend to have
similar outcome functions. We simulate the potential outcomes for domain $s$
as follows,
$Y(0)|X,s\sim\mathcal{N}(\sum\nolimits_{j=1}^{15}\sigma(\bar{x}^{s}_{1})f_{j}^{0}(X),1)$;
$Y(1)|X,s\sim\mathcal{N}(\sum\nolimits_{j=1}^{15}\sigma(\bar{x}^{s}_{1})f_{j}^{1}(X),1)$.

We evenly split the data into 10 folds for variability assessment with results
summarized in Table 2. The drawbacks of the tree-based CF approach are
demonstrated as the increased feature dimension exacerbates decrease of domain
overlap. MDCN demonstrates the best performance in 9 domains, outperforming
benchmarks which conduct either domain-level adjustment (MDMN) or treatment-
level adjustment (DCFR) but not both. DCFR outperforms CFR with statistical
significance, which further supports adjusting for treatment imbalance
differentially in different domains. This is also reflected in the combination
of CFR and MDMN (MDMNCFR), which results in a drop in performance when
compared to MDMN. We visualize the learned feature embeddings $\phi(x)$ in
Figure S2, shedding light on how MDCN predicts on unlabeled target domains. We
first find the three closest “neighborhoods” to the target domain using the
MDCN weight matrix and then plot their relative positions. The target domain
in MDCN is located at the center of its three neighbors in Figure 4(c), which
enables it to borrow information from their respective outcome predictions.
Additionally, the treatment and control groups from the three neighbors are
all well-balanced, making MDCN more robust against treatment group imbalances.
In Figure 4(a), DCFR regards two of its neighbors as irrelevant domains,
preventing it from learning from more reference cases. In Figure 4(b), the
treatment group imbalance is not well controlled in MDMN, thus resulting in
less competitive generalization to the target domain.

##  7 Discussion

Here, we propose MDCN, an approach that provides robust estimation of the
conditional treatment effect for data collected from multi-center
observational studies with particular emphasis in accurately inferring
treatment effects from individuals from a new center. Our approach addresses a
lack of methods research with respects to this specific but not uncommon
problem. More importantly, our model demonstrates the potential to provide
patients from new, unobserved centers with better references for treatment.
Additionally, we also elaborate on the design of MDCN and underpin it with
supporting theory. In empirical evaluations, MDCN consistently outperforms
state-of-the-art domain adaptation methods and causal methods. One limitation
of MDCN lies in the scalability with respects to domain size $S$ which comes
at a cost of $\mathcal{O}(S^{2})$. To address this scalability issue, one
could combine similar domains via some similarity heuristic. Likewise,
scalability becomes an issue when handling multiple or continuous treatment
regimes. As the number of treatment groups increases, learning a balanced
treatment embedding tends to be more challenging. Adequately addressing these
issues provides interesting directions for future research in this area.
Moving forward, our techniques may be applied to multicenter observational
studies and neuromodulation; the ethics of neuromodulation are complex, and we
do not view our study as directly impacting these discussions. Finally, as
with any developed methodology with relevance to healthcare, the methods
proposed here should not be solely relied upon to guide clinical practice,
suggest treatments, or recommend any health-related actions without
consultation and close collaboration with medical professionals.

## Acknowledgments

Research reported in this manuscript was supported by the National Institute
of Biomedical Imaging and Bioengineering and the National Institute of Mental
Health through the National Institutes of Health BRAIN Initiative under Award
Number R01EB026937 and the National Institute Of Mental Health of the National
Institutes of Health under Award Number R01MH125430. The content is solely the
responsibility of the authors and does not necessarily represent the official
views of the National Institutes of Health.

## References

  * Assaad et al. [2021] Serge Assaad, Shuxi Zeng, Chenyang Tao, Shounak Datta, Nikhil Mehta, Ricardo Henao, Fan Li, and Lawrence Carin.  Counterfactual representation learning with balancing weights.  _Artificial Intelligence and Statistics_ , 2021.
  * Athey and Imbens [2016] Susan Athey and Guido Imbens.  Recursive partitioning for heterogeneous causal effects.  _Proceedings of the National Academy of Sciences_ , 113(27):7353–7360, 2016.
  * Awan et al. [2020] Usaid Awan, Marco Morucci, Vittorio Orlandi, Sudeepa Roy, Cynthia Rudin, and Alexander Volfovsky.  Almost-matching-exactly for treatment effect estimation under network interference.  In _International Conference on Artificial Intelligence and Statistics_ , pages 3252–3262, 2020.
  * Azizzadenesheli et al. [2019] Kamyar Azizzadenesheli, Anqi Liu, Fanny Yang, and Animashree Anandkumar.  Regularized learning for domain adaptation under label shifts.  In _International Conference on Learning Representations_ , 2019.
  * Bailey and Crawley [2009] Kathleen R Bailey and Jacqueline N Crawley.  Anxiety-related behaviors in mice.  _Methods of Behavior Analysis in Neuroscience_ , 2009.
  * Ben-David et al. [2010] Shai Ben-David, John Blitzer, Koby Crammer, Alex Kulesza, Fernando Pereira, and Jennifer Wortman Vaughan.  A theory of learning from different domains.  _Machine Learning_ , 79(1-2):151–175, 2010.
  * Borgwardt et al. [2006] Karsten M. Borgwardt, Arthur Gretton, Malte J. Rasch, Hans-Peter Kriegel, Bernhard Schölkopf, and Alex J. Smola.  Integrating structured biological data by kernel maximum mean discrepancy.  _Bioinformatics_ , 22(14):e49–e57, 2006.
  * Chang and Dy [2017] Yale Chang and Jennifer Dy.  Informative subspace learning for counterfactual inference.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 31, 2017.
  * De Raedt et al. [2015] Rudi De Raedt, Marie-Anne Vanderhasselt, and Chris Baeken.  Neurostimulation as an intervention for treatment resistant depression: From research on mechanisms towards targeted neurocognitive strategies.  _Clinical Psychology Review_ , 41:61–69, 2015.
  * Deaton [2010] Angus Deaton.  Instruments, randomization, and learning about development.  _Journal of Economic Literature_ , 48(2):424–55, 2010.
  * Djork-Arné et al. [2016] Clevert Djork-Arné, Thomas Unterthiner, and Sepp Hochreiter.  Fast and accurate deep network learning by exponential linear units (elus).  In _Proceedings of the International Conference on Learning Representations_ , volume 6, 2016.
  * Drysdale et al. [2017] Andrew T Drysdale, Logan Grosenick, Jonathan Downar, Katharine Dunlop, Farrokh Mansouri, Yue Meng, Robert N Fetcho, Benjamin Zebley, Desmond J Oathes, Amit Etkin, et al.  Resting-state connectivity biomarkers define neurophysiological subtypes of depression.  _Nature Medicine_ , 23(1):28–38, 2017.
  * Duckworth et al. [2010] Angela Lee Duckworth, Eli Tsukayama, and Henry May.  Establishing causality using longitudinal hierarchical linear modeling: An illustration predicting achievement from self-control.  _Social Psychological and Personality Science_ , 1(4):311–317, 2010.
  * Ferrer et al. [2009] Ricard Ferrer, Antonio Artigas, David Suarez, Eduardo Palencia, Mitchell M Levy, Angel Arenzana, Xose Luis Pérez, and Josep-Maria Sirvent.  Effectiveness of treatments for severe sepsis: A prospective, multicenter, observational study.  _American Journal of Respiratory and Critical Care Medicine_ , 180(9):861–866, 2009.
  * Gallagher et al. [2017] Neil Gallagher, Kyle R Ulrich, Austin Talbot, Kafui Dzirasa, Lawrence Carin, and David E Carlson.  Cross-spectral factor analysis.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 30, 2017.
  * Ganin et al. [2016] Yaroslav Ganin, Evgeniya Ustinova, Hana Ajakan, Pascal Germain, Hugo Larochelle, François Laviolette, Mario March, and Victor Lempitsky.  Domain-adversarial training of neural networks.  _Journal of Machine Learning Research_ , 17(59):1–35, 2016.
  * Glorot and Bengio [2010] Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In _Proceedings of the International Conference on Artificial Intelligence and Statistics_ , pages 249–256, 2010.
  * Granger [1969] C. W. J. Granger.  Investigating causal relations by econometric models and cross-spectral methods.  _Econometrica_ , 37(3):424–438, 1969.
  * Gulrajani et al. [2017] Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville.  Improved training of wasserstein gans.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Hassanpour and Greiner [2020] Negar Hassanpour and Russell Greiner.  Variational auto-encoder architectures that excel at causal inference.  _Advances in Neural Information Processing Systems_ , 2020.
  * Hernan and Robins [2020] Miguel A Hernan and James M Robins.  _Causal Inference_.  Boca Raton: Chapman & Hall/CRC, 2020.
  * Hill [2011] Jennifer L Hill.  Bayesian nonparametric modeling for causal inference.  _Journal of Computational and Graphical Statistics_ , 20(1):217–240, 2011.
  * Hoffman et al. [2018] Judy Hoffman, Eric Tzeng, Taesung Park, Jun-Yan Zhu, Phillip Isola, Kate Saenko, Alexei Efros, and Trevor Darrell.  CyCADA: Cycle-consistent adversarial domain adaptation.  In Jennifer Dy and Andreas Krause, editors, _Proceedings of the International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pages 1989–1998. PMLR, 2018.
  * Jain [2002] Kewal K Jain.  Personalized medicine.  _Current Opinion in Molecular Therapeutics_ , 4(6):548–558, 2002.
  * Kingma and Ba [2015] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _Proceedings of the International Conference on Learning Representations_ , 2015.
  * Kuang et al. [2017] Kun Kuang, Peng Cui, Bo Li, Meng Jiang, Shiqiang Yang, and Fei Wang.  Treatment effect estimation with data-driven variable decomposition.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 31, 2017.
  * Künzel et al. [2019] Sören R Künzel, Jasjeet S Sekhon, Peter J Bickel, and Bin Yu.  Metalearners for estimating heterogeneous treatment effects using machine learning.  _Proceedings of the National Academy of Sciences_ , 116(10):4156–4165, 2019.
  * Li et al. [2018a] Fan Li, Kari Lock Morgan, and Alan M Zaslavsky.  Balancing covariates via propensity score weighting.  _Journal of the American Statistical Association_ , 113(521):390–400, 2018a.
  * Li et al. [2018b] Yitong Li, michael Murias, geraldine Dawson, and David E Carlson.  Extracting relationships by multi-domain matching.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018b.
  * Li et al. [2019a] Yitong Li, Michael Murias, Samantha Major, Geraldine Dawson, and David Carlson.  On target shift in adversarial domain adaptation.  In _Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics_ , volume 89 of _Proceedings of Machine Learning Research_ , pages 616–625, 2019a.
  * Li et al. [2019b] Yitong Li, Michael Murias, Samantha Major, Geraldine Dawson, and David Carlson.  On target shift in adversarial domain adaptation.  In _International Conference on Artificial Intelligence and Statistics_ , pages 616–625. PMLR, 2019b.
  * Maaten [2008] Laurens van der Maaten.  Visualizing datausing t-sne.  _Journal of Machine Learning Research_ , 9:2579–2605, 2008\.
  * Mansour et al. [2008] Yishay Mansour, Mehryar Mohri, and Afshin Rostamizadeh.  Domain adaptation with multiple sources.  In D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou, editors, _Advances in Neural Information Processing Systems_ , volume 21. Curran Associates, Inc., 2008.
  * McGilchrist [1994] CA McGilchrist.  Estimation in generalized mixed models.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 56(1):61–69, 1994.
  * Oring et al. [2021] Alon Oring, Zohar Yakhini, and Yacov Hel-Or.  Autoencoder image interpolation by shaping the latent space.  In _Proceedings of the International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pages 8281–8290. PMLR, 2021.
  * Reed and MarksII [1999] Russell Reed and Robert J MarksII.  _Neural smithing: Supervised learning in feedforward artificial neural networks_.  Mit Press, 1999.
  * Rosenbaum and Rubin [1983] Paul R Rosenbaum and Donald B Rubin.  The central role of the propensity score in observational studies for causal effects.  _Biometrika_ , 70(1):41–55, 1983.
  * Rubin and Thomas [2000] Donald B Rubin and Neal Thomas.  Combining propensity score matching with additional adjustments for prognostic covariates.  _Journal of the American Statistical Association_ , 95(450):573–585, 2000.
  * Scangos et al. [2021] Katherine W Scangos, Ankit N Khambhati, Patrick M Daly, Ghassan S Makhoul, Leo P Sugrue, Hashem Zamanian, Tony X Liu, Vikram R Rao, Kristin K Sellers, Heather E Dawes, et al.  Closed-loop neuromodulation in an individual with treatment-resistant depression.  _Nature medicine_ , 27(10):1696–1700, 2021.
  * Schmidhuber [2015] Jürgen Schmidhuber.  Deep learning in neural networks: An overview.  _Neural Networks_ , 61:85–117, 2015.
  * Shalit et al. [2017] Uri Shalit, Fredrik D Johansson, and David Sontag.  Estimating individual treatment effect: generalization bounds and algorithms.  In _International Conference on Machine Learning_ , pages 3076–3085. PMLR, 2017.
  * Shen et al. [2018] Jian Shen, Yanru Qu, Weinan Zhang, and Yong Yu.  Wasserstein distance guided representation learning for domain adaptation.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 32(1), 2018.
  * Shi et al. [2019] Claudia Shi, David Blei, and Victor Veitch.  Adapting neural networks for the estimation of treatment effects.  In _Advances in Neural Information Processing Systems_ , pages 2507–2517, 2019.
  * Sohrab [2003] Houshang H Sohrab.  _Basic Real Analysis_ , volume 231.  Springer, 2003.
  * Stanslaski et al. [2012] Scott Stanslaski, Pedram Afshar, Peng Cong, Jon Giftakis, Paul Stypulkowski, Dave Carlson, Dave Linde, Dave Ullestad, Al-Thaddeus Avestruz, and Timothy Denison.  Design and validation of a fully implantable, chronic, closed-loop neuromodulation device with concurrent sensing and stimulation.  _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , 20(4):410–421, 2012.
  * Suk et al. [2021] Youmi Suk, Hyunseung Kang, and Jee-Seon Kim.  Random forests approach for causal inference with clustered observational data.  _Multivariate Behavioral Research_ , 56(6):829–852, 2021.
  * Sun and Saenko [2016] Baochen Sun and Kate Saenko.  Deep coral: Correlation alignment for deep domain adaptation.  In Gang Hua and Hervé Jégou, editors, _Computer Vision – ECCV 2016 Workshops_ , pages 443–450, Cham, 2016. Springer International Publishing.
  * Sun et al. [2016] Baochen Sun, Jiashi Feng, and Kate Saenko.  Return of frustratingly easy domain adaptation.  In _Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence_ , page 2058–2065. AAAI Press, 2016.
  * Tachet des Combes et al. [2020] Remi Tachet des Combes, Han Zhao, Yu-Xiang Wang, and Geoffrey J Gordon.  Domain adaptation with conditional distribution matching and generalized label shift.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 19276–19289. Curran Associates, Inc., 2020.
  * Tibshirani et al. [2020] Julie Tibshirani, Susan Athey, and Stefan Wager.  _grf: Generalized Random Forests_ , 2020.  URL https://CRAN.R-project.org/package=grf.  R package version 1.2.0.
  * Villani [2008] Cédric Villani.  _Optimal Transport: Old and New_ , volume 338.  Springer Science & Business Media, 2008.
  * Welch [1967] P. Welch.  The use of fast fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms.  _IEEE Transactions on Audio and Electroacoustics_ , 15(2):70–73, 1967.
  * Zhang and Lu [2012] Guoyi Zhang and Yan Lu.  Bias-corrected random forests in regression.  _Journal of Applied Statistics_ , 39(1):151–160, 2012.
  * Zhang et al. [2013] Kun Zhang, Bernhard Schölkopf, Krikamol Muandet, and Zhikun Wang.  Domain adaptation under target and conditional shift.  In _Proceedings of the International Conference on International Conference on Machine Learning - Volume 28_ , page 819–827, 2013.
  * Zhao et al. [2019] Han Zhao, Remi Tachet Des Combes, Kun Zhang, and Geoffrey Gordon.  On learning invariant representations for domain adaptation.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, _Proceedings of the International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pages 7523–7532. PMLR, 2019.

## Checklist

  1. 1.

For all authors…

    1. (a)

Do the main claims made in the abstract and introduction accurately reflect
the paper’s contributions and scope? [Yes]

    2. (b)

Did you describe the limitations of your work? [Yes] In Section 7.

    3. (c)

Did you discuss any potential negative societal impacts of your work? [Yes] In
Section 7.

    4. (d)

Have you read the ethics review guidelines and ensured that your paper
conforms to them? [Yes]

  2. 2.

If you are including theoretical results…

    1. (a)

Did you state the full set of assumptions of all theoretical results? [Yes] In
Sections 3, 4 and A.

    2. (b)

Did you include complete proofs of all theoretical results? [Yes] In Section
A.

  3. 3.

If you ran experiments…

    1. (a)

Did you include the code, data, and instructions needed to reproduce the main
experimental results (either in the supplemental material or as a URL)? [Yes]
In the supplemental material.

    2. (b)

Did you specify all the training details (e.g., data splits, hyperparameters,
how they were chosen)? [Yes] In Section B.1.

    3. (c)

Did you report error bars (e.g., with respect to the random seed after running
experiments multiple times)? [Yes] In Table 1 and 2.

    4. (d)

Did you include the total amount of compute and the type of resources used
(e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] In Section
B.1.

  4. 4.

If you are using existing assets (e.g., code, data, models) or
curating/releasing new assets…

    1. (a)

If your work uses existing assets, did you cite the creators? [Yes] In Section
B.1 and 6.2.

    2. (b)

Did you mention the license of the assets? [Yes] In Section B.1 and 6.2.

    3. (c)

Did you include any new assets either in the supplemental material or as a
URL? [Yes] In the supplemental material.

    4. (d)

Did you discuss whether and how consent was obtained from people whose data
you’re using/curating? [N/A] The assets were approved by Institutional Animal
Care and Use Committee (IACUC).

    5. (e)

Did you discuss whether the data you are using/curating contains personally
identifiable information or offensive content? [N/A]

  5. 5.

If you used crowdsourcing or conducted research with human subjects…

    1. (a)

Did you include the full text of instructions given to participants and
screenshots, if applicable? [N/A]

    2. (b)

Did you describe any potential participant risks, with links to Institutional
Review Board (IRB) approvals, if applicable? [N/A]

    3. (c)

Did you include the estimated hourly wage paid to participants and the total
amount spent on participant compensation? [N/A]

##  Appendix A Proofs for Theoretical Analysis in Section 5

We heuristically explain in the main article that the proposed method is
inspired by an error bound on the target domain. Here, we provide details in
its derivation. To make this section self-dependent, we repeat some
definitions, and further comment on their implied properties.

Definition 1 (Probabilistic Discrepancy). For two hypotheses
$h,h^{\prime}:\mathbb{R}^{P}\rightarrow\mathbb{R}$, their difference based on
a probabilistic distribution $D$ over $\mathcal{X}$ is defined as
$\gamma(h,h^{\prime}|D)=\mathbb{E}_{x\sim D}|h(x)-h^{\prime}(x)|$.

Definition 1 is used to quantify the distance between any two hypotheses based
on a distribution $D$. This discrepancy is symmetric as
$\gamma(h,h^{\prime}|D)=\gamma(h^{\prime},h|D)$. If we replace $h$ by the
ground truth outcome function, it represents the error of $h^{\prime}$.

Definition 2 (Lipschitz Continuity). A function
$f:\mathbb{R}^{P}\rightarrow\mathbb{R}$ is Lipschitz continuous with parameter
$\lambda$, if $|f(x_{1})-f(x_{2})|\leq\lambda||x_{1}-x_{2}||_{2}$ holds for
any vectors $x_{1},x_{2}\in\mathcal{X}$. We denote the family as
$\mathcal{F}_{\lambda}$.

In this family, for any $\lambda$, $f\in\mathcal{F}_{\lambda}$
$\implies-f\in\mathcal{F}_{\lambda}$. We assume that the proposed hypothesis
functions $\\{h_{0},h_{1}\\}\in F_{\lambda}$ which can be represented by
neural networks, and the true functions $\\{g_{s,0},g_{s,1}\\}_{s=1}^{S}\in
F_{\lambda^{*}}$.

The family of Lipschitz functions with parameter $\lambda$ can also measure
the Wasserstein distance between spaces. Its relationship to the Wasserstein-1
distance is described below.

$\sup\limits_{f:\mathcal{X}\rightarrow\mathcal{R},||f||_{L}<\lambda}\mathbb{E}_{x\sim
D}[f(x)]-\mathbb{E}_{x\sim D^{\prime}}[f(x)]=\lambda W_{1}(D,D^{\prime})$ (10)

######  Lemma S1 (Additivity of Lipschitz functions).

If $f_{1}\in F_{\lambda_{1}}$ and $f_{2}\in F_{\lambda_{2}}$, then
$f_{1}+f_{2}\in F_{\lambda_{1}+\lambda_{2}}$.

This can be shown trivially with the triangle inequality. By the symmetry of
the Lipschitz family, $f_{1}-f_{2},-f_{1}+f_{2},-f_{1}-f_{2}\in
F_{\lambda_{1}+\lambda_{2}}$.

######  Lemma S2 (Symmetry of the Wasserstein distance).

For any two probabilistic distributions $D$ and $D^{\prime}$ over
$\mathcal{X}$, $W_{1}(D,D^{\prime})=W_{1}(D^{\prime},D)$.

###### Proof.

Assume that $f^{\prime}\in F_{1}$ is the function that maximizes
$\mathbb{E}_{x\sim D}[f(x)]-\mathbb{E}_{x\sim D}[f(x)]$. We have
$\mathbb{E}_{x\sim D}[f^{\prime}(x)]-\mathbb{E}_{x\sim
D}[f^{\prime}(x)]=W_{1}(D,D^{\prime})$. Since
$f\in\mathcal{F}_{1}\implies-f^{\prime}\in\mathcal{F}_{1}$,

$W_{1}(D,D^{\prime})=\mathbb{E}_{x\sim
D^{\prime}}[-f^{\prime}(x)]-\mathbb{E}_{x\sim D}[-f^{\prime}(x)]\leq
W_{1}(D^{\prime},D).$

From the other direction, we can similarly show that $W_{1}(D^{\prime},D)\geq
W_{1}(D,D^{\prime})$, which leads to
$W_{1}(D^{\prime},D)=W_{1}(D,D^{\prime})$.

∎

To bound the error on the target domain, we take a two-step procedure. First,
we bound the error of the unobserved potential outcomes within each domain
given the observed data. Then we explore the cross-domain relationships to
bound the overall error on the target domain with labeled data from source
domains.

Between-treatment Bound

For any source domain $s$ with $D_{s}$, we want to minimize the distance
between our proposed hypotheses $h_{0},h_{1}$ and $g_{s,0},g_{s,1}$. Since we
are only able to observe $T=0$ and $T=1$ in $D_{s,0}$, and $D_{i,1}$,
respectively. The full potential outcome error on $D_{s}$ cannot be directly
calculated. Instead, Proposition S1 suggests that this error can be bounded.

######  Proposition S1.

For any source domain $s$ with the marginal probabilities of receiving
treatment and control as $p_{s}^{T=1}$ and $p_{s}^{T=0}$, the probabilistic
discrepancy between $\\{h_{0},h_{1}\\}$ and $\\{g_{s,0},g_{s,1}\\}$,
$\gamma(h_{0},g_{i,0}|D_{s})+\gamma(h_{1},g_{i,1}|D_{s})$ has an upper bound,

$\small\begin{split}\gamma(h_{0},g_{s,0}|D_{s})+\gamma(h_{1},g_{s,1}|D_{s})\leq\gamma(h_{0},g_{s,0}|D_{s,0})+\gamma(h_{1},g_{i,1}|D_{s,1})+{(\lambda+\lambda^{*})}W_{1}(D_{s,0},D_{s,1}).\end{split}$

###### Proof.

For ease of derivation, we only expand on the error bound of the control group
$T=0$, $\gamma(h_{0},g_{s,0}|D_{s})$. In T-learner, $h_{0}$ and $h_{1}$ are
constructed separately, and they do not influence each other. The bound on
$T=1$ can be derived with identical steps.

$\small\begin{split}\gamma(h_{0},g_{s,0}|D_{s})&=\mathbb{E}_{x\sim
D_{s}}|h_{0}(x)-g_{s,0}(x)|\\\ &=p_{s}^{T=0}\mathbb{E}_{x\sim
D_{s,0}}|h_{0}(x)-g_{s,0}(x)|+p_{s}^{T=1}\mathbb{E}_{x\sim
D_{s,1}}|h_{0}(x)-g_{s,0}(x)|\\\ &=p_{s}^{T=0}\mathbb{E}_{x\sim
D_{s,0}}|h_{0}(x)-g_{s,0}(x)|+p_{s}^{T=1}\mathbb{E}_{x\sim
D_{s,0}}|h_{0}(x)-g_{s,0}(x)|\\\ &\quad\quad\quad\quad-
p_{s}^{T=1}\mathbb{E}_{x\sim
D_{s,0}}|h_{0}(x)-g_{s,0}(x)|+p_{s}^{T=1}\mathbb{E}_{x\sim
D_{s,1}}|h_{0}(x)-g_{s,0}(x)|\\\
&=\gamma(h_{0},g_{s,0}|D_{s,0})+p_{s}^{T=1}\biggl{(}\mathbb{E}_{x\sim
D_{s,1}}|h_{0}(x)-g_{s,0}(x)|-\mathbb{E}_{x\sim
D_{s,0}}|h_{0}(x)-g_{s,0}(x)|\biggl{)}\\\
&\leq\gamma(h_{0},g_{s,0}|D_{s,0})+p_{s}^{T=1}\biggl{(}\sup\limits_{f:\mathcal{X}\rightarrow\mathcal{R},||f||_{L}<\lambda+\lambda^{*}}\mathbb{E}_{x\sim
D_{i,0}}[f(x)]-\mathbb{E}_{x\sim D_{s,1}}[f(x)]\biggl{)}\\\
&=\gamma(h_{0},g_{s,0}|D_{s,0})+p_{s}^{T=1}(\lambda+\lambda^{*})W_{1}(D_{s,0},D_{s,1})\end{split}$

In the fifth line, we use the additivity of the Lipschitz family,
$h_{0}-g_{s,0}\in\mathcal{F}_{\lambda+\lambda^{*}}$. Similarly, between
$h_{1}$ and $g_{s,1}$, we have

$\gamma(h_{1},g_{s,1}|D_{i})\leq\gamma(h_{1},g_{s,1}|D_{s,1})+p_{s}^{T=0}(\lambda+\lambda^{*})W_{1}(D_{s,1},D_{s,0}).$

With the symmetry of $W(\cdot,\cdot)$,
$(\lambda+\lambda^{*})W_{1}(D_{s,1},D_{s,0})=(\lambda+\lambda^{*})W_{1}(D_{s,0},D_{s,1})$.
Then we sum up the two groups, and the bound is proven.

∎

Cross-Domain Bound

In this step, we simplify the problem to a multiple domain setups but assuming
randomization within each domain, $D_{s,0}=D_{s,1}=D_{s}$. We keep the other
conditions unchanged such as the domain shift and the shift in outcome
functions. As the outcomes on the target domain are unobservable, the error
can still only be bounded rather than calculated. Proposition S2 suggests a
bound.

Before we state the bound, we define two quantities. We assume that
$w=\\{w_{s}\\}_{s=1}^{S-1}$ are weights on labeled source domains and they sum
up to 1. $\gamma^{*}_{0}$ and $\gamma^{*}_{0}$ define the minimum
discrepancies attainable for the weighted summation of all domains.

$\begin{split}&\gamma^{*}_{0}=\min\limits_{h_{0}}\biggl{[}\gamma(h_{0},g_{S,0}|D_{S})+\sum\limits_{s=1}^{S-1}w_{s}\gamma(h_{0},g_{s,0}|D_{s})\biggl{]};\\\
&\gamma^{*}_{1}=\min\limits_{h_{1}}\biggl{[}\gamma(h_{1},g_{S,1}|D_{S})+\sum\limits_{s=1}^{S-1}w_{s}\gamma(h_{1},g_{s,1}|D_{s})\biggl{]}.\end{split}$
(11)

In (11), the two quantities depict the fundamental difference in true outcome
functions across all domains, which is unoptimizable. The larger its value,
the harder the target domain error can be controlled. If $g_{s,0},g_{s,1}$
only differ slightly across domains, we are more likely to find some
universally defined $h_{0}^{*}$ and $h_{1}^{*}\in\mathcal{F}_{\lambda^{*}}$
that make $\gamma^{*}_{0}$ and $\gamma^{*}_{1}$ small.

######  Proposition S2.

Assume that we have full randomization within each domain, and $h_{0}^{*}$ and
$h_{1}^{*}\in\mathcal{F}_{\lambda^{*}}$ with minimum errors $\gamma^{*}_{0}$
and $\gamma^{*}_{1}$ in (11). For any positive weights on source domains
$w=\\{w_{s}\\}_{s=1}^{S-1}$ with $\sum_{s=1}^{S-1}w_{s}=1$. The target domain
error can be bounded by,

$\begin{split}\gamma(h_{0},g_{S,0}|D_{S})+\gamma(h_{1},g_{S,1}|D_{S})&\leq
2(\lambda+\lambda^{*})W_{1}(D_{S},\sum\limits_{s=1}^{S-1}w_{s}D_{s})+\gamma_{0}^{*}+\gamma_{1}^{*}\\\
&+\sum\limits_{s=1}^{S-1}w_{s}[\gamma(h_{0},g_{0,s}|D_{s})+\gamma(h_{1},g_{1,s}|D_{s})].\end{split}$
(12)

###### Proof.

Again, we prove it on the control group and extend it to both groups with the
property of “T-learner".

$\begin{aligned}
\gamma(h_{0},g_{S,0}|D_{S})&\leq\gamma(h_{0},h_{0}^{*}|D_{S})+\gamma(g_{S,0},h_{0}^{*}|D_{S})\text{~{}
( Triangle Inequality )}\\\
&=\gamma(h_{0},h_{0}^{*}|D_{S})+\gamma(g_{S,0},h_{0}^{*}|D_{S})+\sum\limits_{s=1}^{S-1}w_{s}\gamma(h_{0},h_{0}^{*}|D_{s})-\sum\limits_{s=1}^{S-1}w_{s}\gamma(h_{0},h_{0}^{*}|D_{s})\\\
&\leq\gamma(h_{0},h_{0}^{*}|D_{S})+[\gamma(g_{S,0},h_{0}^{*}|D_{S})+\sum\limits_{s=1}^{S-1}w_{s}\gamma(g_{s,0},h_{0}^{*}|D_{s})]\\\
&\quad\quad\quad\quad\quad\quad+\sum\limits_{s=1}^{S-1}w_{s}\gamma(h_{0},g_{s,0}|D_{s})-\sum\limits_{s=1}^{S-1}w_{s}\gamma(h_{0},h_{0}^{*}|D_{s})\\\
&=\gamma(h_{0},h_{0}^{*}|D_{S})-\sum\limits_{s=1}^{S-1}w_{s}\gamma(h_{0},h_{0}^{*}|D_{s})+\gamma_{0}^{*}+\sum\limits_{s=1}^{S-1}w_{s}\gamma(h_{0},g_{0,s}|D_{s})\\\
&\leq(\lambda+\lambda^{*})W_{1}(D_{S},\sum\limits_{s=1}^{S-1}w_{s}D_{s})+\gamma_{0}^{*}+\sum\limits_{s=1}^{S-1}w_{s}\gamma(h_{0},g_{0,s}|D_{s})~{}~{}(\text{as
}h_{0}-h_{0}^{*}\in\mathcal{F}_{\lambda+\lambda^{*}})\\\ \end{aligned}$

Likewise, we can obtain the following quantity for the treatment group.

$\gamma(h_{1},g_{S,1}|D_{S})\leq(\lambda+\lambda^{*})W_{1}(D_{S},\sum\limits_{s=1}^{S-1}w_{s}D_{s})+\gamma_{1}^{*}+\sum\limits_{s=1}^{S-1}w_{s}\gamma(h_{1},g_{1,s}|D_{s})$

Lastly, we sum over two groups to get the bound in (12). ∎

Overall Bound

Proposition S1 or Proposition S2 only partially address our big goal. However,
their combination gives us the overall bound on the target domain in a
multicenter observational study setup.

Theorem 2 For any positive weights on source domain
$w=\\{w_{s}\\}_{s=1}^{S-1}$ with $\sum_{s=1}^{S-1}w_{s}=1$. The overall target
domain error for the proposed hypothesis functions $\\{h_{0},h_{1}\\}$ can be
bounded by,

$\small\begin{split}\gamma(h_{0},g_{S,0}|D_{S})+\gamma(h_{1},g_{S,1}|D_{S})&\leq(\lambda+\lambda^{*})[2W_{1}(D_{S},\sum\limits_{s=1}^{S-1}w_{s}D_{s})+\sum\limits_{s=1}^{S-1}w_{s}W_{1}(D_{s,0},D_{s,1})]\\\
&+\sum\limits_{s=1}^{S-1}w_{s}[\gamma(h_{0},g_{s,0}|D_{s,0})+\gamma(h_{1},g_{i,1}|D_{s,1})]+\gamma_{0}^{*}+\gamma_{1}^{*}.\end{split}$

###### Proof.

We replace each $\gamma(h_{0},g_{0,s}|D_{s})$ and
$\gamma(h_{1},g_{1,s}|D_{s})$ in Proposition S2 with the bound derived in
Proposition S1, and the overall bound can be obtained. ∎

##  Appendix B Method Implementations and Data Preprocessing

###  B.1 Model Specifications

Algorithm 1 MDCN Algorithm

Data with treatment labels and observed outcomes,
$\\{\\{t_{i},x_{i},y_{i}\\}_{i=1}^{N_{s}}\\}_{s=1}^{S-1}$ from source domains
and unlabeled data from the target domain $\\{x_{i}\\}_{i=1}^{N_{S}}$.

Feature embedding function $\phi$, potential outcome models
$\\{h_{0},h_{1}\\}$, domain discriminator $f_{cd}$ and treatment discriminator
$f_{bt}$ .

for iter $=1$ to $n^{iter}$ do

Sample a mini-batch from source domains
$\\{\\{t_{i},x_{i},y_{i}\\}_{i=1}^{N_{s}}\\}_{s=1}^{S-1}$ and target domain
$\\{x\\}_{i=1}^{N_{S}}$.

Calculate $\\{\boldsymbol{l}_{s}\\}_{s=1}^{S}$ and update
$\\{\boldsymbol{w}_{s}\\}_{s=1}^{S}$ in (6) for cross-domain adjustment.

for iter $=1$ to $n_{1}$ do

Sample a mini-batch from source domains
$\\{\\{t_{i},x_{i},y_{i}\\}_{i=1}^{N_{s}}\\}_{s=1}^{S-1}$ and target domain
$\\{x\\}_{i=1}^{N_{S}}$.

Optimize full loss in (1) with respect to the parameters in
$\phi,h_{0},h_{1}$.

Update the parameters in $\phi,h_{0},h_{1}$.

end for

for iter $=1$ to $n_{2}$ do

Sample a mini-batch from source domains
$\\{\\{t_{i},x_{i},y_{i}\\}_{i=1}^{N_{s}}\\}_{s=1}^{S-1}$ and target domain
$\\{x\\}_{i=1}^{N_{S}}$.

Optimize the cross-domain $-L_{cd}(\phi,f_{cd})$ in (7) with respect to the
parameters in $f_{cd}$.

Update the parameters in $f_{cd}$.

Optimize the between-treatment $-L_{bt}(\phi,f_{bt})$ in (4) with respect to
the parameters in $f_{bt}$.

Update the parameters in $f_{bt}$.

end for

end for

For MLP, CFR, MDMN, MDNMCFR and MDCN, we use the same model architecture for a
fair comparison. The feature embedding network is parametrized through a
neural network with two hidden layers of 50 units each. The outcome networks
for both treatment and control groups are based on a neural network with two
hidden layers of 50 units. For the circular data example and the mice data
example, we set the dimensions of the embedded space to 10 and 20,
respectively. The Lipschitz constraints for the Wasserstein-1 distance is
maintained through the gradient penalty [19]. The tuning parameters have fixed
values, as $\alpha$=1e-3 and $\beta$=5e-4. In practice, these specifications
could also be tuned. We use ELU as the activation function [11] and ADAM as
the optimizer [25] with step size 1e-4 throughout both experiments. The full
algorithm of MDCN is sketched out in Algorithm 1. The code will be public on
Github with the MIT license when the manuscript is accepted.

The method CF is implemented through the R package grf with GPL-3 license
[50]. We let the model automatically tune its hyper-parameters by specifying
`tune.parameters="all"`.

All Python- based methods are run on a single NVIDIA P100 GPU; the R-based CF
is run on an Intel(R) Xeon(R) Gold 6154 CPU.

###  B.2 Preprocessing the mice data

The raw data have high variability, as the first 100 principle components only
explains 50% of the total variance. Additionally, the limited sample size
poses a challenge with regards to the variability assessment. To enlarge the
sample size and decrease the variability of the data. We use a recently
published auto-encoder for data interpolation [35]. The encoder
$f_{e}:\mathbb{R}^{9856}\rightarrow\mathbb{R}^{15}$ is a neural network that
has three hidden layers of 100 units each. The decoder
$f_{d}:\mathbb{R}^{15}\rightarrow\mathbb{R}^{9856}$ is a neural network with
one hidden layer of 50 units. Then, we construct a discriminative network:
$f_{i}:\mathbb{R}^{15}\rightarrow[0,1]$ with a single hidden layer of 50 units
to judge the qualify of interpolated data. Lastly, we add a network with no
hidden layer to predict mice identity
$f_{m}:\mathbb{R}^{15}\rightarrow[0,1]^{21}$.

For each domain $s$, we have the following loss for the control group (home
cage),

$\begin{aligned} L_{s,0}=&\mathbb{E}_{x_{1},x_{2}\sim{D_{s,0}},\alpha\sim
U(0,1)}\min_{f_{e},f_{d},f_{m}}\max_{f_{i}}\biggl{\\{}\\\
&\biggl{[}(f_{d}(f_{e}(x_{1}))-x_{1})^{2}+(f_{d}(f_{e}(x_{2}))-x_{2})^{2}\biggl{]}+~{}\hskip
199.16928pt\text{(I)}\\\ &.1\biggl{[}f_{e}(f_{d}(\alpha
f_{e}(x_{1})+(1-\alpha)f_{e}(x_{2})))-(\alpha
f_{e}(x_{1})+(1-\alpha)f_{e}(x_{2}))\biggl{]}^{2}+~{}\hskip
93.89418pt\text{(II)}\\\
&.05\biggl{[}{CE}(f_{m}(f_{e}(x_{1})),s))+{CE}(f_{m}(f_{e}(x_{2})),s))+{CE}(f_{m}(\alpha
f_{e}(x_{1})+(1-\alpha)f_{e}(x_{2}),s))\biggl{]}+~{}\hskip
5.69046pt\text{(III)}\\\
&.05\biggl{[}\log(f_{i}(f_{e}(x_{1})))+\log(f_{i}(f_{e}(x_{2})))+\log(1-f_{i}(\alpha
f_{e}(x_{1})+(1-\alpha)f_{e}(x_{2})))\biggl{]}\biggl{\\}}.~{}\hskip
42.67912pt\text{(IV)}\end{aligned}$

We have the data reconstruction in (I). (II) is the “cycle loss" [35], which
is to reconstruct the interpolated feature embedding. It also contributes to
make $f_{e}$ and $f_{d}$ as reciprocal of each other. (III) is the cross-
entropy loss for each domain to be predictive of its domain identity $s$. It
is helpful in maintaining domain-specific information and creating domain
shift. (IV) is the adversarial loss. It discriminates how similar interpolated
points are to the raw data. Enforcing it enables the interpolation to retain
the original data distributions. The full loss is calculated by summing over
all domains and all treatment groups, $\sum_{s=1}^{21}[L_{s,0}+L_{s,1}]$. We
train this framework for 12,000 iterations. In each iteration, we randomly
draw 20 samples from each treatment group and each domain, which makes a total
batch size of 840. After the training is finished, we observe that the
embedded space explains >30% of the total variance. we randomly interpolate
10,000 samples for each domain and each treatment group to make an enlarged
dataset. These interpolated points overall reaches >90% accuracy in predicting
the mouse identity, with $>60\%$ confidence indicated by the discriminative
network $f_{i}$ that they are from the true data.

For outcome models, we prepare 30 mini neural networks with input and output
size: $\mathbb{R}^{15}\rightarrow\mathbb{R}$,
$\\{f_{1}^{0},\cdots,f_{15}^{0}\\}$, $\\{f_{1}^{1},\cdots,f_{15}^{1}\\}$.
These neural networks all have a single hidden layer with 15 units, and
hyperbolic tangent activation functions. Their weights are initialized through
the Xavier initialization [17]. Then we use 15 evenly spaced points from
$[-1.5,1.5]$ to shift their weight parameters to make them distinct.
Explicitly, we shift the weight parameters in
$f_{s}^{0}\in\\{f_{1}^{0},\cdots,f_{15}^{0}\\}$ by $+(3s/15-1.5)$. For
$\\{f_{1}^{1},\cdots,f_{15}^{1}\\}$, we shift them by another permutation of
these 15 points. Lastly, they are combined with the the center of each domain
$s$, $\bar{x}^{s}=\\{\bar{x}^{s}_{1},\cdots,\bar{x}^{s}_{15}\\}$ to create the
overall shift in outcome functions. It encourages domains that are close in
center to have similar outcome functions,

$\displaystyle
Y(0)|X,s\sim\mathcal{N}(\sum\limits_{j=1}^{15}sigmoid(\bar{x}^{s}_{1})f_{j}^{0}(X),1);~{}~{}Y(1)|X,s\sim\mathcal{N}(\sum\limits_{j=1}^{15}sigmoid(\bar{x}^{s}_{1})f_{j}^{1}(X),1).$

##  Appendix C Additional Visualizations

(a) MDMN

(b) MDMNCFR

(c) DCFR

Figure S1:  Additional visualizations of domain-level similarities. DCFR in
1(c) better maintains the relative closeness across domains than CFR does in
Figure 2. However, it is still not appropriately capturing the original
domain-level similarity as shown in 1(a) and 1(b), both of which are backed by
the cross-domain adjustment.

(a) CFR

(b) MDMNCFR

Figure S2:  Additional t-SNE plots for learned $\phi(x)$. The domain 9 is
chosen as the target domain and is labeled as orange "T". For convenience, its
3 nearest neighboring domains are labeled as "0", "1" and "2". They are found
through the weight matrix in MDMN. Labels with solid shade and transparent
shade represent the control group and the treatment group, respectively.

Figure S1 visualizes the similarity at domain level for the other 3 methods.
For methods without the domain level adjustment (CFR included), we use a
surrogate to calculate $\boldsymbol{l}_{s}$ in (6), so that similarity could
be represented in a comparable manner. We first use 1-d t-SNE on the learned
feature embedding $\phi(x)$, as $\text{tsne}(\phi(x))$. Then we calculate the
distance and weight according to (13),

$\begin{split}&\boldsymbol{l}_{s}=\\{l_{s,i}\\}_{i\neq
s},~{}l_{s,i}=\mathbb{E}_{x\sim{D}_{s}}[\text{tsne}(\phi(x))]-\mathbb{E}_{x\sim{D}_{i}}[\text{tsne}(\phi(x))],\\\
&\boldsymbol{w}_{s}=\\{w_{s,i}\\}_{i\neq
s}=\text{softmax}(-\boldsymbol{l}_{s}).\end{split}$ (13)

Diagonal elements are set to 1 as reference values. DCFR visualized in Figure
S1 better maintains the original domain-level connectivity than CFR in Figure
3. It connect such as domains 3 and 4 or domains 6 and 7. MDMN and MDMNCFR
both have the cross-domain component to learn domain similarities. They also
clearly reflect the original data distributions by showing strong connectivity
to two neighboring domains.

Figure S2 includes two additional t-SNE plots of the learned feature $\phi(x)$
in mouse data. In 2(a), the target domain is distant to the neighboring
domains "0" and "2", which gives the model fewer relevant points to learn from
to predict on the target domain. In 2(b), the treatment group and control
group are still less balanced compared to that in Figure 4(c), this suggests
worse generalization from our bound in Theorem 1.