\section{Electronic Submission}
\label{submission}

Submission to ICML 2022 will be entirely electronic, via a web site
(not email). Information about the submission process and \LaTeX\ templates
are available on the conference web site at:
\begin{center}
\textbf{\texttt{http://icml.cc/}}
\end{center}

The guidelines below will be enforced for initial submissions and
camera-ready copies. Here is a brief summary:
\begin{itemize}
\item Submissions must be in PDF\@. 
\item \textbf{New to this year}: If your paper has appendices, submit the appendix together with the main body and the references \textbf{as a single file}. Reviewers will not look for appendices as a separate PDF file. So if you submit such an extra file, reviewers will very likely miss it.
\item Page limit: The main body of the paper has to be fitted to 8 pages, excluding references and appendices; the space for the latter two is not limited. For the final version of the paper, authors can add one extra page to the main body.
\item \textbf{Do not include author information or acknowledgements} in your
    initial submission.
\item Your paper should be in \textbf{10 point Times font}.
\item Make sure your PDF file only uses Type-1 fonts.
\item Place figure captions \emph{under} the figure (and omit titles from inside
    the graphic file itself). Place table captions \emph{over} the table.
\item References must include page numbers whenever possible and be as complete
    as possible. Place multiple citations in chronological order.
\item Do not alter the style template; in particular, do not compress the paper
    format by reducing the vertical spaces.
\item Keep your abstract brief and self-contained, one paragraph and roughly
    4--6 sentences. Gross violations will require correction at the
    camera-ready phase. The title should have content words capitalized.
\end{itemize}

\subsection{Submitting Papers}

\textbf{Paper Deadline:} The deadline for paper submission that is
advertised on the conference website is strict. If your full,
anonymized, submission does not reach us on time, it will not be
considered for publication. 

\textbf{Anonymous Submission:} ICML uses double-blind review: no identifying
author information may appear on the title page or in the paper
itself. \cref{author info} gives further details.

\textbf{Simultaneous Submission:} ICML will not accept any paper which,
at the time of submission, is under review for another conference or
has already been published. This policy also applies to papers that
overlap substantially in technical content with conference papers
under review or previously published. ICML submissions must not be
submitted to other conferences and journals during ICML's review
period.
Informal publications, such as technical
reports or papers in workshop proceedings which do not appear in
print, do not fall under these restrictions.

\medskip

Authors must provide their manuscripts in \textbf{PDF} format.
Furthermore, please make sure that files contain only embedded Type-1 fonts
(e.g.,~using the program \texttt{pdffonts} in linux or using
File/DocumentProperties/Fonts in Acrobat). Other fonts (like Type-3)
might come from graphics files imported into the document.

Authors using \textbf{Word} must convert their document to PDF\@. Most
of the latest versions of Word have the facility to do this
automatically. Submissions will not be accepted in Word format or any
format other than PDF\@. Really. We're not joking. Don't send Word.

Those who use \textbf{\LaTeX} should avoid including Type-3 fonts.
Those using \texttt{latex} and \texttt{dvips} may need the following
two commands:

{\footnotesize
\begin{verbatim}
dvips -Ppdf -tletter -G0 -o paper.ps paper.dvi
ps2pdf paper.ps
\end{verbatim}}
It is a zero following the ``-G'', which tells dvips to use
the config.pdf file. Newer \TeX\ distributions don't always need this
option.

Using \texttt{pdflatex} rather than \texttt{latex}, often gives better
results. This program avoids the Type-3 font problem, and supports more
advanced features in the \texttt{microtype} package.

\textbf{Graphics files} should be a reasonable size, and included from
an appropriate format. Use vector formats (.eps/.pdf) for plots,
lossless bitmap formats (.png) for raster graphics with sharp lines, and
jpeg for photo-like images.

The style file uses the \texttt{hyperref} package to make clickable
links in documents. If this causes problems for you, add
\texttt{nohyperref} as one of the options to the \texttt{icml2022}
usepackage statement.


\subsection{Submitting Final Camera-Ready Copy}

The final versions of papers accepted for publication should follow the
same format and naming convention as initial submissions, except that
author information (names and affiliations) should be given. See
\cref{final author} for formatting instructions.

The footnote, ``Preliminary work. Under review by the International
Conference on Machine Learning (ICML). Do not distribute.'' must be
modified to ``\textit{Proceedings of the
$\mathit{39}^{th}$ International Conference on Machine Learning},
Baltimore, Maryland, USA, PMLR 162, 2022.
Copyright 2022 by the author(s).''

For those using the \textbf{\LaTeX} style file, this change (and others) is
handled automatically by simply changing
$\mathtt{\backslash usepackage\{icml2022\}}$ to
$$\mathtt{\backslash usepackage[accepted]\{icml2022\}}$$
Authors using \textbf{Word} must edit the
footnote on the first page of the document themselves.

Camera-ready copies should have the title of the paper as running head
on each page except the first one. The running title consists of a
single line centered above a horizontal rule which is $1$~point thick.
The running head should be centered, bold and in $9$~point type. The
rule should be $10$~points above the main text. For those using the
\textbf{\LaTeX} style file, the original title is automatically set as running
head using the \texttt{fancyhdr} package which is included in the ICML
2022 style file package. In case that the original title exceeds the
size restrictions, a shorter form can be supplied by using

\verb|\icmltitlerunning{...}|

just before $\mathtt{\backslash begin\{document\}}$.
Authors using \textbf{Word} must edit the header of the document themselves.

\section{Format of the Paper}

All submissions must follow the specified format.

\subsection{Dimensions}




The text of the paper should be formatted in two columns, with an
overall width of 6.75~inches, height of 9.0~inches, and 0.25~inches
between the columns. The left margin should be 0.75~inches and the top
margin 1.0~inch (2.54~cm). The right and bottom margins will depend on
whether you print on US letter or A4 paper, but all final versions
must be produced for US letter size.
Do not write anything on the margins.

The paper body should be set in 10~point type with a vertical spacing
of 11~points. Please use Times typeface throughout the text.

\subsection{Title}

The paper title should be set in 14~point bold type and centered
between two horizontal rules that are 1~point thick, with 1.0~inch
between the top rule and the top edge of the page. Capitalize the
first letter of content words and put the rest of the title in lower
case.

\subsection{Author Information for Submission}
\label{author info}

ICML uses double-blind review, so author information must not appear. If
you are using \LaTeX\/ and the \texttt{icml2022.sty} file, use
\verb+\icmlauthor{...}+ to specify authors and \verb+\icmlaffiliation{...}+ to specify affiliations. (Read the TeX code used to produce this document for an example usage.) The author information
will not be printed unless \texttt{accepted} is passed as an argument to the
style file.
Submissions that include the author information will not
be reviewed.

\subsubsection{Self-Citations}

If you are citing published papers for which you are an author, refer
to yourself in the third person. In particular, do not use phrases
that reveal your identity (e.g., ``in previous work \cite{langley00}, we
have shown \ldots'').

Do not anonymize citations in the reference section. The only exception are manuscripts that are
not yet published (e.g., under submission). If you choose to refer to
such unpublished manuscripts \cite{anonymous}, anonymized copies have
to be submitted
as Supplementary Material via CMT\@. However, keep in mind that an ICML
paper should be self contained and should contain sufficient detail
for the reviewers to evaluate the work. In particular, reviewers are
not required to look at the Supplementary Material when writing their
review (they are not required to look at more than the first $8$ pages of the submitted document).

\subsubsection{Camera-Ready Author Information}
\label{final author}

If a paper is accepted, a final camera-ready copy must be prepared.
For camera-ready papers, author information should start 0.3~inches below the
bottom rule surrounding the title. The authors' names should appear in 10~point
bold type, in a row, separated by white space, and centered. Author names should
not be broken across lines. Unbolded superscripted numbers, starting 1, should
be used to refer to affiliations.

Affiliations should be numbered in the order of appearance. A single footnote
block of text should be used to list all the affiliations. (Academic
affiliations should list Department, University, City, State/Region, Country.
Similarly for industrial affiliations.)

Each distinct affiliations should be listed once. If an author has multiple
affiliations, multiple superscripts should be placed after the name, separated
by thin spaces. If the authors would like to highlight equal contribution by
multiple first authors, those authors should have an asterisk placed after their
name in superscript, and the term ``\textsuperscript{*}Equal contribution"
should be placed in the footnote block ahead of the list of affiliations. A
list of corresponding authors and their emails (in the format Full Name
\textless{}email@domain.com\textgreater{}) can follow the list of affiliations.
Ideally only one or two names should be listed.

A sample file with author names is included in the ICML2022 style file
package. Turn on the \texttt{[accepted]} option to the stylefile to
see the names rendered. All of the guidelines above are implemented
by the \LaTeX\ style file.

\subsection{Abstract}

The paper abstract should begin in the left column, 0.4~inches below the final
address. The heading `Abstract' should be centered, bold, and in 11~point type.
The abstract body should use 10~point type, with a vertical spacing of
11~points, and should be indented 0.25~inches more than normal on left-hand and
right-hand margins. Insert 0.4~inches of blank space after the body. Keep your
abstract brief and self-contained, limiting it to one paragraph and roughly 4--6
sentences. Gross violations will require correction at the camera-ready phase.

\subsection{Partitioning the Text}

You should organize your paper into sections and paragraphs to help
readers place a structure on the material and understand its
contributions.

\subsubsection{Sections and Subsections}

Section headings should be numbered, flush left, and set in 11~pt bold
type with the content words capitalized. Leave 0.25~inches of space
before the heading and 0.15~inches after the heading.

Similarly, subsection headings should be numbered, flush left, and set
in 10~pt bold type with the content words capitalized. Leave
0.2~inches of space before the heading and 0.13~inches afterward.

Finally, subsubsection headings should be numbered, flush left, and
set in 10~pt small caps with the content words capitalized. Leave
0.18~inches of space before the heading and 0.1~inches after the
heading.

Please use no more than three levels of headings.

\subsubsection{Paragraphs and Footnotes}

Within each section or subsection, you should further partition the
paper into paragraphs. Do not indent the first line of a given
paragraph, but insert a blank line between succeeding ones.

You can use footnotes\footnote{Footnotes
should be complete sentences.} to provide readers with additional
information about a topic without interrupting the flow of the paper.
Indicate footnotes with a number in the text where the point is most
relevant. Place the footnote in 9~point type at the bottom of the
column in which it appears. Precede the first footnote in a column
with a horizontal rule of 0.8~inches.\footnote{Multiple footnotes can
appear in each column, in the same order as they appear in the text,
but spread them across columns and pages if possible.}

\begin{figure}[ht]
\vskip 0.2in
\begin{center}
\centerline{\includegraphics[width=\columnwidth]{icml_numpapers}}
\caption{Historical locations and number of accepted papers for International
Machine Learning Conferences (ICML 1993 -- ICML 2008) and International
Workshops on Machine Learning (ML 1988 -- ML 1992). At the time this figure was
produced, the number of accepted papers for ICML 2008 was unknown and instead
estimated.}
\label{icml-historical}
\end{center}
\vskip -0.2in
\end{figure}

\subsection{Figures}

You may want to include figures in the paper to illustrate
your approach and results. Such artwork should be centered,
legible, and separated from the text. Lines should be dark and at
least 0.5~points thick for purposes of reproduction, and text should
not appear on a gray background.

Label all distinct components of each figure. If the figure takes the
form of a graph, then give a name for each axis and include a legend
that briefly describes each curve. Do not include a title inside the
figure; instead, the caption should serve this function.

Number figures sequentially, placing the figure number and caption
\emph{after} the graphics, with at least 0.1~inches of space before
the caption and 0.1~inches after it, as in
\cref{icml-historical}. The figure caption should be set in
9~point type and centered unless it runs two or more lines, in which
case it should be flush left. You may float figures to the top or
bottom of a column, and you may set wide figures across both columns
(use the environment \texttt{figure*} in \LaTeX). Always place
two-column figures at the top or bottom of the page.

\subsection{Algorithms}

If you are using \LaTeX, please use the ``algorithm'' and ``algorithmic''
environments to format pseudocode. These require
the corresponding stylefiles, algorithm.sty and
algorithmic.sty, which are supplied with this package.
\cref{alg:example} shows an example.

\begin{algorithm}[tb]
   \caption{Bubble Sort}
   \label{alg:example}
\begin{algorithmic}
   \STATE {\bfseries Input:} data $x_i$, size $m$
   \REPEAT
   \STATE Initialize $noChange = true$.
   \FOR{$i=1$ {\bfseries to} $m-1$}
   \IF{$x_i > x_{i+1}$}
   \STATE Swap $x_i$ and $x_{i+1}$
   \STATE $noChange = false$
   \ENDIF
   \ENDFOR
   \UNTIL{$noChange$ is $true$}
\end{algorithmic}
\end{algorithm}

\subsection{Tables}

You may also want to include tables that summarize material. Like
figures, these should be centered, legible, and numbered consecutively.
However, place the title \emph{above} the table with at least
0.1~inches of space before the title and the same after it, as in
\cref{sample-table}. The table title should be set in 9~point
type and centered unless it runs two or more lines, in which case it
should be flush left.


\begin{table}[t]
\caption{Classification accuracies for naive Bayes and flexible
Bayes on various data sets.}
\label{sample-table}
\vskip 0.15in
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{lcccr}
\toprule
Data set & Naive & Flexible & Better? \\
\midrule
Breast    & 95.9$\pm$ 0.2& 96.7$\pm$ 0.2& $\surd$ \\
Cleveland & 83.3$\pm$ 0.6& 80.0$\pm$ 0.6& $\times$\\
Glass2    & 61.9$\pm$ 1.4& 83.8$\pm$ 0.7& $\surd$ \\
Credit    & 74.8$\pm$ 0.5& 78.3$\pm$ 0.6&         \\
Horse     & 73.3$\pm$ 0.9& 69.7$\pm$ 1.0& $\times$\\
Meta      & 67.1$\pm$ 0.6& 76.5$\pm$ 0.5& $\surd$ \\
Pima      & 75.1$\pm$ 0.6& 73.9$\pm$ 0.5&         \\
Vehicle   & 44.9$\pm$ 0.6& 61.5$\pm$ 0.4& $\surd$ \\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

Tables contain textual material, whereas figures contain graphical material.
Specify the contents of each row and column in the table's topmost
row. Again, you may float tables to a column's top or bottom, and set
wide tables across both columns. Place two-column tables at the
top or bottom of the page.

\subsection{Theorems and such}
The preferred way is to number definitions, propositions, lemmas, etc. consecutively, within sections, as shown below.
\begin{definition}
\label{def:inj}
A function $f:X \to Y$ is injective if for any $x,y\in X$ different, $f(x)\ne f(y)$.
\end{definition}
Using \cref{def:inj} we immediate get the following result:
\begin{proposition}
If $f$ is injective mapping a set $X$ to another set $Y$, 
the cardinality of $Y$ is at least as large as that of $X$
\end{proposition}
\begin{proof} 
Left as an exercise to the reader. 
\end{proof}
\cref{lem:usefullemma} stated next will prove to be useful.
\begin{lemma}
\label{lem:usefullemma}
For any $f:X \to Y$ and $g:Y\to Z$ injective functions, $f \circ g$ is injective.
\end{lemma}
\begin{theorem}
\label{thm:bigtheorem}
If $f:X\to Y$ is bijective, the cardinality of $X$ and $Y$ are the same.
\end{theorem}
An easy corollary of \cref{thm:bigtheorem} is the following:
\begin{corollary}
If $f:X\to Y$ is bijective, 
the cardinality of $X$ is at least as large as that of $Y$.
\end{corollary}
\begin{assumption}
The set $X$ is finite.
\label{ass:xfinite}
\end{assumption}
\begin{remark}
According to some, it is only the finite case (cf. \cref{ass:xfinite}) that is interesting.
\end{remark}

\subsection{Citations and References}

Please use APA reference format regardless of your formatter
or word processor. If you rely on the \LaTeX\/ bibliographic
facility, use \texttt{natbib.sty} and \texttt{icml2022.bst}
included in the style-file package to obtain this format.

Citations within the text should include the authors' last names and
year. If the authors' names are included in the sentence, place only
the year in parentheses, for example when referencing Arthur Samuel's
pioneering work \yrcite{Samuel59}. Otherwise place the entire
reference in parentheses with the authors and year separated by a
comma \cite{Samuel59}. List multiple references separated by
semicolons \cite{kearns89,Samuel59,mitchell80}. Use the `et~al.'
construct only for citations with three or more authors or after
listing all authors to a publication in an earlier reference \cite{MachineLearningI}.

Authors should cite their own work in the third person
in the initial version of their paper submitted for blind review.
Please refer to \cref{author info} for detailed instructions on how to
cite your own papers.

Use an unnumbered first-level section heading for the references, and use a
hanging indent style, with the first line of the reference flush against the
left margin and subsequent lines indented by 10 points. The references at the
end of this document give examples for journal articles \cite{Samuel59},
conference publications \cite{langley00}, book chapters \cite{Newell81}, books
\cite{DudaHart2nd}, edited volumes \cite{MachineLearningI}, technical reports
\cite{mitchell80}, and dissertations \cite{kearns89}.

Alphabetize references by the surnames of the first authors, with
single author entries preceding multiple author entries. Order
references for the same authors by year of publication, with the
earliest first. Make sure that each reference includes all relevant
information (e.g., page numbers).

Please put some effort into making references complete, presentable, and
consistent, e.g. use the actual current name of authors.
If using bibtex, please protect capital letters of names and
abbreviations in titles, for example, use \{B\}ayesian or \{L\}ipschitz
in your .bib file.

\section*{Accessibility}
Authors are kindly asked to make their submissions as accessible as possible for everyone including people with disabilities and sensory or neurological differences.
Tips of how to achieve this and what to pay attention to will be provided on the conference website \url{http://icml.cc/}.

\section*{Software and Data}

If a paper is accepted, we strongly encourage the publication of software and data with the
camera-ready version of the paper whenever appropriate. This can be
done by including a URL in the camera-ready copy. However, \textbf{do not}
include URLs that reveal your institution or identity in your
submission for review. Instead, provide an anonymous URL or upload
the material as ``Supplementary Material'' into the CMT reviewing
system. Note that reviewers are not required to look at this material
when writing their review.

\section*{Acknowledgements}

\textbf{Do not} include acknowledgements in the initial version of
the paper submitted for blind review.

If a paper is accepted, the final camera-ready version can (and
probably should) include acknowledgements. In this case, please
place such acknowledgements in an unnumbered section at the
end of the paper. Typically, this will include thanks to reviewers
who gave useful comments, to colleagues who contributed to the ideas,
and to funding agencies and corporate sponsors that provided financial
support.


\nocite{langley00}


\section{Introduction} \label{sec:introduction}


Estimating individualized treatment effects has been a hot topic because of its wide applications, ranging from personalized medicine, policy research, to customized marketing advertisement. Treatment effects of certain subgroups within the population are often of interest. Recently, there has been an explosion of research devoted to improving estimation and inference of covariate-specific treatment effects, or conditional average treatment effects (CATE) at a target research site \citep{athey2016recursive,wager2018estimation,hahn2020bayesian,kunzel2019metalearners,nie2020quasioracle}. 
However, due to the limited sample size in a single study, improving the accuracy of the estimation of treatment effects remains challenging. 

Leveraging data and models from various research sites to conduct statistical analyses is becoming increasingly popular \citep{reynolds2020leveraging, cohen2020leveraging, berger2015optimizing}.
Distributed research networks have been established in many large scale studies \citep{fleurence2014launching,hripcsak2015observational,platt2018fda,donohue2021use}.  %
A question often being asked is whether additional data or models from other research sites could bring improvement to a local estimation task, especially when a single site does not have enough data to achieve a desired statistical precision. 
This concern is mostly noticeable in estimating treatment effects where sample size requirement is high yet observations are typically limited. 
Furthermore, information exchange between data sites is often highly restricted due to privacy, feasibility, or other concerns, prohibiting centralized analyses that pool data from multiple sources \citep{maro2009design,brown2010distributed,toh2011comparative,raghupathi2014big,deshazo2015comparison,donahue2018veterans,dayan2021federated}.
One way to tackle this challenge is through model averaging \citep{raftery1997bayesian}, where multiple research sites collectively contribute to the tasks of statistical modeling without sharing sensitive subject-level data.
Although this idea has existed in supervised learning problems  \citep{dai2011greedy,mcmahan2017communication}, 
to our best knowledge, there are no established model averaging approach and theoretical results on estimating CATE in a distributed environment. The extension is non-trivial because CATE is unobserved in nature, as opposed to prediction problems where labels are given. 






This paper focuses on improving the prediction accuracy of CATE concerning a target site by leveraging models derived from other sites 
where 
\textit{transportability}  \citep[to be formally defined in Section~\ref{sec:assumption},][]{pearl2011transportability,stuart2011use,pearl2014external,bareinboim2016causal,buchanan2018generalizing,dahabreh2019generalizing} may not hold.
Specifically, there may exist heterogeneity in treatment effects. In the context of our multi-hospital example, these are: 
1) \textbf{local heterogeneity}: within a hospital, patients with different characteristics may have different treatment effects.
This is the traditional notion of CATE;
and 2) \textbf{global heterogeneity}: where the same patient may experience different treatment effects at different hospitals. The second type of heterogeneity is driven by site-level confounding, and hampers the transportability of models across hospital sites. 
We also note that these two types of heterogeneity may interact with each other in the sense that transportability is dependent on patient characteristics, which we will address.

We propose a model averaging framework that uses a flexible tree-based weighting scheme to combine learned models from sites that takes into account heterogeneity. The contribution of each learned model to the target site depends on subject characteristics. This is achieved by applying tree splittings \citep{breiman1984classification} at both the site and the subject levels. 
For example, effects of a treatment in two hospitals may be similar for female patients but not for male, suggesting us to consider borrowing information across sites only on selective subgroups. 
Our approach extends the classic model averaging framework \citep{raftery1997bayesian,wasserman2000bayesian,hansen2007least,yang2001adaptive} by allowing data-adaptive weights, which are interpretable in a sense that they can be used to lend credibility to transportability. 
For example, in the case of extreme heterogeneity where other sites merely contribute to the target, the weights can be used as a diagnostic tool to inform the decision against borrowing information.

\textbf{Main contributions.}
{{1)}} We propose a model averaging scheme with interpretable weights that are adaptive to both local and global heterogeneity via tree-splitting dedicated to improving CATE estimation under distributed data networks. 
{{2)}} We generalize model averaging techniques to study the transportability of causal inference. Causal assumptions with practical implications are explored to warrant the use of our approach. 
{{3)}} We provide an extensive empirical evaluation of the proposed approach with a concrete real-data example on how to apply the method in practice. 
{{4)}} Compared to other distributed learning methods, the proposed framework enables causal analysis without sharing subject-level data, is easy to implement, offers ease of operations, and minimizes infrastructure, which facilitates practical collaboration within research networks. 




\section{Related Work} \label{sec:related}

There are two types of construct of a distributed database \citep{breitbart1986database}: \emph{homogeneous} versus \emph{heterogeneous}. 
For homogeneous data sources, data across sites are random samples of the global population. 
Recent modeling approaches \citep{lin2010relative, lee2017communication, mcmahan2017communication, battey2018distributed, jordan2018communication, tang2020distributed,wang2021tributarypca} all assume samples are randomly partitioned, which guarantees
identical data distribution across sites. 
The goal of these works is to improve overall prediction by averaging results from homogeneous sample divisions. 
The classic random effects meta-analysis (see, e.g.,  \citet{whitehead2002meta,sutton2000methods,borenstein2011introduction} describes heterogeneity using modeling assumptions, but its focus mostly is still on global patterns. 


\textbf{Heterogeneous models.} In practice, however, there is often too much global heterogeneity 
in a distributed data network 
to warrant direct aggregation of models obtained from local sites. The focus shifts to improving the estimation of a target site by selectively leveraging information from other data sources.
There are two main classes of approaches.
The first class is based on comparison of the learned model parameters $\{\widehat\bm{\theta}_1,\dots,\widehat\bm{\theta}_K\}$ from $K$ different sites where for site $k$ we adopt model $f_k(\bm{x}) = f(\bm{x}; \bm{\theta}_k)$ with subject features $\bm{x}$ to approximate the outcome of interest $Y$. 
Clustering and shrinkage approaches are then used
by merging data or models that are similar \citep{ke2015homogeneity,tang2016fused,smith2017federated,ma2017concave,wang2020sylvester,tang2020individualized,tang2021poststratification}. 
Most of these require the pooling of subject-level data.
The second class of approaches falls in the \textit{model averaging} framework \citep{raftery1997bayesian} with weights directly associated with the local prediction. 
Let site 1 be our target site, and the goal is to improve $f_{1}$ using a weighted estimator $f^*(\bm{x}) = \sum_{k=1}^K {\omega}_{k} f_k(\bm{x})$ with weights $\omega_k$ to balance the contribution of each model and $\sum_k \omega_{k} = 1$. 
It provides an immediate interpretation of usefulness of each data source. 
When the weights are proportional to the prediction performance of $f_k$ on site 1, for example, 
$${\omega}_{k} =\frac{\exp\{- \sum_{i \in \mathcal{I}_1}(f_k(\bm{x}_i) - y_i)^2\}}{ \sum_{\ell=1}^{K} \exp\{- \sum_{i \in \mathcal{I}_1}(f_\ell(\bm{x}_i) - y_i)^2\} },$$
with $y_i$ being the observed outcome of subject $i$ in site 1, indexed by $\mathcal{I}_1$, the method is termed as the exponential weighted model averaging (EWMA).
Several variations of ${\omega}_{k}$ can be found in
\citet{yang2001adaptive, dai2011greedy,yao2018using,dai2018bayesian}. 
In general, separate samples are used to obtain the estimates of $\omega_k$'s and $f_k$'s, respectively. 

Here we focus on the literature review of model averaging. 
We note that our framework is also related to federated learning \citep{mcmahan2017communication}. But the latter often
involves iterative updating  rather than a one-shot procedure, and could be hard to apply to nonautomated distributed
research networks. Besides, it has been developed mainly to estimate a global prediction model by leveraging distributed data, and is not designed to target any specific site. We further discuss these approaches and other related research topics and their distinctions with model averaging in Appendix~\ref{suppl-related}.


\textbf{Transportability.}
In causal inference, there is a lot of interest in identifying subgroups with enhanced treatment effects, targeting  at  the  feasibility  of  customizing  estimates  for  individuals \citep{athey2016recursive,wager2018estimation,hahn2020bayesian,kunzel2019metalearners,nie2020quasioracle}. These methods aim to estimate the CATE function $\tau(\bm{x})$, denoting the difference in potential outcomes between treatment and control, conditional on subject characteristics $\bm{x}$. 
To reduce uncertainty in estimation of personalized treatment effects, incorporating additional data or models are sought after. 
\citet{pearl2011transportability,pearl2014external,bareinboim2016causal} introduced the notion of transportability to warrant causal inference models be generalized to a new population. 
The issue of generalizability is common in practice due to the non-representative sampling of participants in randomized controlled trials \citep{cook2002experimental,druckman2011cambridge,allcott2015site,stuart2015assessing,egami2020elements}. 
Progress on bridging the findings from an experimental study with observational data can be found in, e.g., \citet{stuart2015assessing,kern2016assessing,stuart2018generalizability,ackerman2019implementing,yang2020elastic,harton2021combining}. 
See \citet{tipton2018review,colnet2020causal,degtiar2021review} and references therein for a comprehensive review. 
However, most methods require fully centralized data. In contrast, we leverage the distributed nature of model averaging to derive an integrative CATE estimator.













\section{A Tree-based Model Averaging Framework} \label{sec:method}


We first formally define the {conditional average treatment effect} (CATE). 
Let $Y$ denote the outcome of interest, $Z \in \{0,1\}$ denote a binary treatment indicator, and $\bm{X}$ denote subject features. Correspondingly, let $y$, $z$ and $\bm{x}$ denote their realizations. 
Using the potential outcome framework \citep{neyman1923applications,rubin1974estimating}, we define
CATE as $\tau(\bm{x})=E[Y^{(Z=1)}-Y^{(Z=0)} |\bm{X}=\bm{x}],$
where $Y^{(Z=1)}$ and $Y^{(Z=0)}$ are the potential outcomes under treatment arms $Z=1$ and $Z=0$, respectively. 
The expected difference of the potential outcomes is dependent on subject features $\bm{X}$. By the causal consistency assumption, the observed outcome is $Y = ZY^{(Z=1)} + (1 - Z)Y^{(Z=0)}$.


Now suppose the distributed data network 
consists of $K$ sites, each with sample size of $n_k$. 
Site $k$ contains data $\mathcal{D}_k = \{y_i, z_i, \bm{x}_i\}_{i \in \mathcal{I}_k}$, where $\mathcal{I}_k$ denotes its index set.
Its CATE function is given by $\tau_k (\bm{x}) = E_k[Y^{(Z=1)}-Y^{(Z=0)} |\bm{X}=\bm{x}],$ where the expectation is taken over the data distribution in site $k$. 
Without loss of generality, we assume the goal is to estimate the CATE function in site 1, $\tau_1$.





\subsection{Causal Assumptions} \label{sec:assumption}

To ensure information can be properly borrowed across sites, we first impose the following idealistic assumptions, and then present relaxed version of Assumption 2. 
Let $S$ be the site indicator taking values in $\mathcal{S} = \{1, \dots, K\}$ such that $S_i = k$ if $i \in \mathcal{I}_k$. 

\textit{Assumption 1:} $\{Y^{(Z=0)}, Y^{(Z=1)}\} \perp Z | \bm{X}, S$;

\textit{Assumption 2 (Transportability):} $$
\{Y^{(Z=0)}, Y^{(Z=1)}\} \perp S | \bm{X};$$
\textit{Assumption 3:} $	0< P(S = 1| \bm{X}) < 1 \mbox{ and } 0< P(Z = 1| \bm{X}, S) < 1   \quad \mbox{for all } \bm{X} \mbox{ and } S$.


Assumption 1
ensures treatment effects are unconfounded within sites so that $\tau_k(\bm{x})$ can be consistently identified. 
It holds by design when data are randomized controlled trials or when treatment assignment depends on $\bm{X}$. 
By this assumption, we have $\tau_{k}(\bm{x}) =
E[Y |\bm{X}=\bm{x}, S = k, Z = 1] - E[Y |\bm{X}=\bm{x}, S = k, Z = 0]$.
The equality directly results from the assumption. 
Assumption 2
essentially states that the CATE functions are transportable, i.e., $\tau_k(\bm{x}) = \tau_{k'}(\bm{x})$ for $k, k'\in \{1,\dots, K\}$. See also \citet{stuart2011use}, \citet{buchanan2018generalizing} and \citet{yang2020elastic} for similar consideration.
This assumption may not be satisfied due to heterogeneity across sites. In other words, site can be a confounder which prevents transporting of CATE functions across sites.
Our method allows Assumption 2
to be violated and use model averaging weights to determine transportability. Explicitly, we consider a relaxed Assumption 2a to hold for a subset of sites that contains site 1.

\textit{Assumption 2a (Partial Transportability):} $$
\{Y^{(Z=0)}, Y^{(Z=1)}\} \perp S_1 | \bm{X}.$$ 
Here, $S_1$ takes values in $\mathcal{S}_1 = \{k: \tau_k(\bm{x}) = \tau_1(\bm{x})\}$ and $\{1\} \subset \mathcal{S}_1 \subset \mathcal{S}$. 
We denote $\mathcal{S}_1$ as the set of transportable sites with regard to site 1. 
Hence, transportability holds across some sites and specific subjects.
In a special case in Section~\ref{sec:simulation} where $\mathcal{S}_1 = \{1\}$, bias may be introduced to by model averaging. However, our approach is still able to exploits the bias and variance trade off to improve estimation. 
Assumption 3
ensures that all subjects are possible to be observed in site 1 and all subjects in all sites are possible to receive either arm of treatment. 
The former ensures a balance of covariates between site 1 population and the population of other sites. 
Violation of either one may result in extrapolation and introduce unwanted bias to the ensemble estimates for site 1. This assumption is also used, e.g., in \citet{stuart2011use}. 





\subsection{Model Ensemble} \label{sec:adaptiveMA}

We consider an adaptive weighting of $\{\tau_1, \dots, \tau_K\}$ by
\begin{equation} \label{eq:agg}
   {\tau}^*(\bm{x}) = \sum_{k=1}^K \omega_{k}(\bm{x}) \tau_k(\bm{x})
\end{equation}
where ${\tau}^*$ is the weighted model averaging estimator. The weight functions $\omega_{k}(\bm{x})$'s are not only site-specific, but also depend on $\bm{x}$, and follow $\sum_{k=1}^{K} \omega_{k}(\bm{x}) = 1$.
It measures the importance of $\tau_k$ in assisting site 1 when subjects with characteristics $\bm{x}$ are of interest. 
We rely on each of the sites to derive their respective $\widehat\tau_k$ from $\mathcal{D}_k$ so that $\mathcal{D}_1, \dots, \mathcal{D}_K$ do not need to be pooled. Only the estimated functions $\{\widehat \tau_2, \dots, \widehat \tau_K\}$ are passed to site 1. 
We will describe the approaches to estimate $\widehat{\tau}_k$ in Section \ref{sec:local}.


A two-stage model averaging approach is proposed.
We first split $\mathcal{D}_1$, the data in the target site, into a training set and an estimation set indexed by $\{i \in \mathcal{I}_1^{(1)}\}$ and $\{i \in \mathcal{I}_1^{(2)}\}$, respectively.
\emph{1) Local stage:} 
Obtain $\widehat{\tau}_1$ from subjects in $\mathcal{I}_1^{(1)}$. 
Obtain $\widehat{\tau}_k$ from local subjects in $\mathcal{I}_k$, $k = 2, \dots, K$. These $\{\widehat{\tau}_k\}_{k=1}^K$ are then passed to site 1 to get $K$ predicted treatment effects for each subject in $\mathcal{I}_1^{(2)}$, resulting in an augmented data set as shown in Figure~\ref{fig:framework}(b).
\emph{2) Ensemble stage:} A tree-based ensemble model is trained on the augmented data by either an ensemble tree (ET) or an ensemble random forest (EF), with the predicted treatment effects from the previous stage, i.e., $\widehat{\tau}_k(\bm{x}_i)$ as the \emph{outcome}. The site indicator $S$ of which local model is used as well as the subject features $\bm{x}_i$ are fed into the ensemble model as \emph{predictors}. The resulting model will be used to compute our proposed model averaging estimator.
Figure~\ref{fig:framework}(a) illustrates a conceptual diagram of the proposed model averaging framework and structure of the augmented data. 
Note the idea of data augmentation has been used in, e.g., computer vision \citep{wang2017effectiveness}, statistical computing \citep{van2001art}, and imbalanced classification \citep{chawla2002smote}. 
Here the technique is being used to construct weights for model averaging, which will be discussed in the following paragraph. 
Algorithm~\ref{algo:code} provides an algorithmic overview. 
Our method has been implemented as an R package \texttt{ifedtree} available on GitHub (\url{https://github.com/ellenxtan/ifedtree}). 


\begin{figure}[tb
\centering
 \begin{subfigure}{0.5\textwidth}
  \centerline{\includegraphics[width=\linewidth]{figures/diagram_final_crop.pdf}}
  \vspace{-0.3cm}
  \caption{}
  \vspace{0.2pc}
 \end{subfigure}
 \begin{subfigure}{.25\textwidth}
  \centerline{\includegraphics[width=\linewidth]{figures/aug_df.pdf}}
  \caption{}
 \end{subfigure}
 \vspace{-0.3cm}
 \caption{(a) Schema of the proposed algorithm. (b) Illustration of the augmented data constructed from the estimation set of site 1.} 
 \label{fig:framework}
\end{figure}




\setlength{\textfloatsep}{0.4cm}
\begin{algorithm}[tb]
\smal
  \caption{Tree-based model averaging for heterogeneous data sources}
  \label{algo:code}
\begin{algorithmic
    \FOR{$k=1$ {\bfseries to} $K$}
    \STATE \hfill\COMMENT{Loop through $K$ sites. Can be run in parallel.}
    \STATE Build a local model using site $k$ data. Site 1 model uses its training set only.
    \ENDFOR
    \FOR{$i \in \mathcal{I}_1^{(2)}$} 
    \STATE \hfill\COMMENT{Loop through subjects in site 1 estimation set.}
        \FOR{$k=1$ {\bfseries to} $K$} 
        \STATE \hfill\COMMENT{Loop through $K$ local models.}
            \STATE Predict $\widehat{\tau}_k(\bm{x}_i)$ using local model $k$.
            \STATE $D_{i, k} = [\bm{x}_i, k, \widehat{\tau}_k(\bm{x}_i)]$.
        \ENDFOR
    \ENDFOR
    \STATE Create augmented site 1 data $\mathfrak{D}_{aug, 1}$ by concatenating $D_{i, k}$ vectors.
   
    \STATE $\widehat{\mathcal{T}}_{\text{EF}} (\bm{x}, s) =$ {\sc EnsembleForest}($\mathfrak{D}_{aug, 1}$)
   
    \STATE \hfill\COMMENT{Or {\sc EnsembleTree} when $B = 1$.}
\end{algorithmic}
\end{algorithm}
\setlength{\floatsep}{0.4cm}


\textbf{Construction of weights.} A tree-based ensemble is constructed to estimate the weighting functions $\{\omega_{k}\}_{k=1}^K$.
Heterogeneity across sites is explained by including the site index into an augmented training set when building trees.
An intuition of our approach is that sites that are split away from site 1 (by tree nodes) are ignored and the sites that fall into the same leaf node are considered homogeneous to site 1 hence contribute to the estimation of $\tau_1(\bm{x})$. A splitting by site may occur in any branches of a tree, resulting in an information sharing scheme across sites that is dependent on $\bm{x}$. 
We construct the ensemble by first creating an augmented data 
    $\mathfrak{D}_{aug, 1} = \{\bm{x}_i, k, \widehat{\tau}_k(\bm{x}_i)\}_{i \in \mathcal{I}_1^{(2)}, k \in \mathcal{S}}$,
for subjects in $\mathcal{I}_1^{(2)}$.
The illustration of this augmented site 1 data is given in Figure~\ref{fig:framework}(b). An ensemble is then trained on this data by either a tree or a random forest,
with the estimated treatment effects $\widehat{\tau}_k(\bm{x}_i)$ as the outcome, and a categorical site indicator of which local model is used along with all subject-level features as predictors, i.e., $(\bm{x}_i, k)$. 
We denote the resulting function as $\mathcal{T}(\bm{x}, s)$ which depends on both $\bm{x}$ and site $s$, specifically, 
$\mathcal{T}_{\text{ET}} (\bm{x}, s)$ and $ \mathcal{T}_{\text{EF}} (\bm{x}, s)$ for ensemble tree (ET) and ensemble forest (EF), respectively. Let $\mathcal{L}(\bm{x},s)$ denote the final partition of the feature space by the tree to which the pair $(\bm{x}, s)$ belongs. The ET estimate based on the augmented site 1 data can be derived by
\begin{align}\label{eq:ET}
    \widehat \mathcal{T}_{\text{ET}} (\bm{x}, s)
    &=  \{|\{(i, k): (\bm{x}_i, k) \in \mathcal{L}(\bm{x}, s)\}_{i \in \mathcal{I}_1^{(2)}, k \in \mathcal{S}}|\}^{-1}  \nonumber\\
    &\sum_{\{(i, k): (\bm{x}_i, k) \in \mathcal{L}(\bm{x}, s)\}_{i \in \mathcal{I}_1^{(2)}, k \in \mathcal{S}}} \widehat{\tau}_k(\bm{x}_i) \nonumber\\
    &=  \sum_{i\in\mathcal{I}_1^{(2)}} \sum_{k=1}^{K} \frac{\mathbbm{1}\{(\bm{x}_i, k)\in \mathcal{L}(\bm{x}, s)\}}{|\mathcal{L}(\bm{x}, s)|} \widehat{\tau}_k(\bm{x}_i).
\end{align}
Intuitively, observations with similar characteristics ($\bm{x}$ and $\bm{x}'$) and from similar sites ($s$ and $s'$) are more likely to fall in the same partition region in the ensemble tree, i.e., $(\bm{x}, s) \in \mathcal{L}(\bm{x}', s')$ or $(\bm{x}', s') \in \mathcal{L}(\bm{x}, s)$. This resembles a \emph{non-smooth kernel} where weights are $1/|\mathcal{L}(\bm{x}, s)|$ for observations that are within the neighborhood of $(\bm{x}, s)$, and 0 otherwise.
The estimator borrows information from neighbors in the space of $\bm{X}$ and $S$. 
The splits of the tree are based on minimizing in-sample MSE of $\widehat\tau$ within each leaf
and pruned by cross-validation over choices of the complexity parameter. 
Since a single tree is prone to be unstable,
in practice, we use random forest to reduce variance and smooth the partitioning boundaries.
By aggregating $B$ ET estimates each based on a subsample of the augmented data, $\{\widehat{\mathcal{T}}^{(b)}\}_{b=1}^B$, an EF estimate can be constructed by
\begin{align}
       \widehat \mathcal{T}_{\text{EF}} (\bm{x}, s) 
& = \frac{1}{B} \sum_{b=1}^{B} \widehat \mathcal{T}^{(b)} (\bm{x}, s) \nonumber \\
& = \sum_{i\in \mathcal{I}_1^{(2)}} \sum_{k=1}^{K} \lambda_{i,k}(\bm{x}, s) \widehat{\tau}_k(\bm{x}_i), \label{eq:EF}\\
    \text{where~} &  \lambda_{i,k}(\bm{x}, s) = \frac{1}{B} \sum_{b=1}^{B} \frac{\mathbbm{1}\{(\bm{x}_i, k)\in \mathcal{L}_b(\bm{x}, s)\}}{|\mathcal{L}_b(\bm{x}, s)|}.\nonumber
\end{align}
The form of $\widehat \mathcal{T}^{(b)} (\bm{x}, s)$ closely follows \eqref{eq:ET} but is based on a subsample of $\mathfrak{D}_{aug, 1}$. The weights, $\lambda_{i,k}(\bm{x}, s)$, are similar to that in \eqref{eq:ET}, and can be viewed as kernel weighting that defines an adaptive neighborhood of $\bm{x}$ and $s$. 
We then obtain the model averaging estimates defined in \eqref{eq:agg} by fixing $s = 1$ such that $\widehat \tau_{\text{ET}}^*(\bm{x}) = \widehat\mathcal{T}_{\text{ET}} (\bm{x}, s=1)$ or $\widehat \tau_{\text{EF}}^*(\bm{x}) = \widehat\mathcal{T}_{\text{EF}} (\bm{x}, s=1)$.
The weight functions $\{\omega_{k}(\bm{x})\}_{k=1}^K$ for $\widehat\tau^*(\bm{x})$ can be immediately obtained from the ET or EF by
\begin{align*}
    \widehat\tau_{\text{ET}}^*(\bm{x}) = \ &  \widehat \mathcal{T}_{\text{ET}} (\bm{x}, 1) = \sum_{k=1}^K \widehat \omega_{k}(\bm{x}) \widehat \tau_k(\bm{x}),  \\
    \text{where~~} & \widehat \omega_{k}(\bm{x}) = \sum_{i\in\mathcal{I}_1^{(2)}} \frac{\mathbbm{1}\{(\bm{x}_i, k)\in \mathcal{L}(\bm{x}, 1)\}}{|\mathcal{L}(\bm{x}, 1)|}; \nonumber \\
   \widehat\tau_{\text{EF}}^*(\bm{x}) = \ &  \widehat \mathcal{T}_{\text{EF}} (\bm{x}, 1) = \sum_{k=1}^K \widehat \omega_{k}(\bm{x}) \widehat \tau_k(\bm{x}),  \\
    \text{where~~} & \widehat \omega_{k}(\bm{x}) = \sum_{i\in\mathcal{I}_1^{(2)}} \lambda_{i,k}(\bm{x}, 1).
\end{align*}
It can be verified that $\sum_{k=1}^{K} \widehat \omega_k(\bm{x}) = 1$ for all $\bm{x}$.
As our simulations in Section \ref{sec:simulation} show, $\widehat \tau^*$ improves the local functional estimate $\widehat\tau_1$.
We set $B=2,000$ throughout the paper.
Tree and forest estimates are obtained by R packages \verb|rpart| and \verb|grf|, respectively. 


\textbf{Interpretability of weights.}
The choice of tree-based models naturally results in such kernel weighting $w_k(\bm{x})$ \citep{athey2019generalized}, which are not accessible by other ensemble techniques.
Such explicit and interpretable weight functions could deliver meaningful rationales for data integration.
For example, under scenarios where there exists extreme global heterogeneity (as shown in Section~\ref{sec:simulation} when $c$ is large), $w_k(\bm{x})$ can be used as a diagnostic tool to decide which external data sources should be co-used. 
Weights close to $0$ inform against model transportability, and they are adaptive to subject-level features $\bm{x}$ so that decisions can be made based on the subpopulations of interest.













\subsection{Local Models: Obtaining $\widehat\tau_k$} \label{sec:local}

Estimate of $\tau_k(\bm{x})$ at each local site must be obtained separately before the ensemble. 
Our proposed ensemble framework can be applied to a general estimator of $\tau_k(\bm{x})$. For each site, the local estimate could be obtained using different methods. 
Recently, there has been many work dedicated to the estimation of individualized treatment effects \citep{athey2016recursive,wager2018estimation,hahn2020bayesian,kunzel2019metalearners,nie2020quasioracle}. 
As an example, we consider using the causal tree (CT) \citep{athey2016recursive} to estimate the local model at each site. CT is a non-linear learner that \textit{(i)} allows different types of outcome such as discrete and continuous, and can be applied to a broad range of real data scenarios; \textit{(ii)} can manage hundreds of features and high order interactions by construction; \textit{(iii)} can be applied to both experimental studies and observational studies by propensity score weighting or doubly robust methods. 
CT is implemented in the R package \verb|causalTree|. We also explore another estimating option for local models in Appendix~\ref{suppl-sec:sim}. 





\subsection{Asymptotic Properties} \label{sec:consist}

We provide consistency guarantee of the proposed estimator $\widehat \mathcal{T}_{\text{EF}}$ for the true target $\tau_1$. 
Assuming point-wise consistent local estimators are used for $\{\tau_k\}_{k=1}^{K}$, EF with subsampling procedure described in Appendix~\ref{suppl-sec:consist} is consistent. 
\begin{theorem}\label{Theorem}
Suppose the subsample used to build each tree in an ensemble forest is drawn from different subjects of the augmented data and the following conditions hold:

  (i) Bounded covariates: Features $\bm{X}_i$ and the site indicator $S_i$ are independent and have a density that is bounded away from 0 and infinity.
  
  (ii) Lipschitz response: the conditional mean function $\mathbb{E}[ \mathcal{T}|\bm{X}=\bm{x},S=1]$ is Lipschitz-continuous.
  
  (iii) Honest trees: trees in the random forest use different data for placing splits and estimating leaf-wise responses.

Then $\widehat \mathcal{T}_{\text{EF}}(\bm{x},1) \overset{p}{\to} \tau_{1}(\bm{x})$, for all $\bm{x}$, as $\min_k n_k \to \infty$. Hence, $\widehat \tau_{\text{EF}}^*(\bm{x}) \overset{p}{\to} \tau_1(\bm{x}).$
\end{theorem}
The conditions and a proof of Theorem \ref{Theorem} is given in Appendix~\ref{suppl-sec:consist}.
To demonstrate the consistency properties of our methods, we add in Appendix~\ref{suppl-sec:sim} oracle versions of ET and EF estimators, denoted as ET-oracle and EF-oracle, which use the ground truth of local models $\{\tau_k\}_{k=1}^K$ in estimating $\{\widehat \omega_k\}_{k=1}^K$.
This removes the uncertainty in local models. 
The remaining uncertainty only results from the estimation of the ensemble weights, and we see both oracle estimators achieve minimal MSE.
Section~\ref{sec:simulation} gives a detailed evaluation of the finite sample performance. 









\begin{figure*}[tb
\centering
 \begin{subfigure}{0.49\textwidth}
  \centerline{\includegraphics[width=\linewidth]{figures/plt_ct500_tau_2grp_fair_mseRatio.pdf}}
  \caption{}
 \end{subfigure}
 \begin{subfigure}{0.49\textwidth}
  \centerline{\includegraphics[width=\linewidth]{figures/plt_ct500_tau_cont_fair_mseRatio.pdf}}
  \caption{}
 \end{subfigure}
 \vspace{-0.3cm}
 \caption{
 Box plots of MSE ratios of CATE estimators, respectively, over LOC, for \textbf{(a) discrete grouping} and \textbf{(b) continuous grouping} across site. 
  Different colors imply different estimators, and x-axis, i.e., the value of $c$, differentiates the scale of global heterogeneity. The red dotted line denotes an MSE ratio of 1. 
  MA performance is truncated due to large MSE ratios. 
The proposed ET and EF achieve smaller MSE ratios compared to standard model averaging or ensemble methods and are robust to heterogeneity across settings. 
}
\vspace{-0.3cm}
\label{fig:sim_box}
\end{figure*}


\section{Simulation Study}
\label{sec:simulation}


Monte Carlo simulations are conducted to assess the proposed methods. We specify $m(\bm{x}, k)$ as the conditional outcome surface and $\tau(\bm{x}, k)$ as the conditional treatment effect for individuals  with features $\bm{x}$ in site $k$. The treatment propensity is specified as $e(\bm{x}) = \Pr(Z=1|\bm{X}=\bm{x})$. 
The potential outcomes can be written as $Y_i = m(\bm{X}_i, S_i) + \{ Z_i - e(\bm{X}_i) \} \tau(\bm{X}_i, S_i) + \epsilon_i$, following notations in \citet{robinson1988root,athey2016recursive,wager2018estimation,nie2020quasioracle}. 
The mean function is $m(\bm{x}, k) = \frac{1}{2} x_1 + \sum_{d=2}^4 x_d +(x_1 - 3) \cdot c \cdot U_k$, and the treatment effect function is specified as 
$$ \tau(\bm{x}, k) = \mathbbm{1}\{x_1 > 0\} \cdot x_1 +(x_1 - 3) \cdot c \cdot U_k,$$
where $z=0,1$, $U_k$ denotes the global heterogeneity due to site-level confounding, controlled by a scaling factor $c$, and $\epsilon_i \sim N(0, 1)$. 
Features follow $\bm{X}_i \sim {N}(\bm{0},\bm{I}_D)$, where $D=5$, and are independent of $\epsilon_i$.
The simulation setting within each site (with $k$ fixed) is motivated by designs in \citet{athey2016recursive}. Features in $\tau$ are determinants of treatment effect while those in $m$ but not in $\tau$ are prognostic only. 
The data are generated under a distributed data networks. 
We assume there are $K=20$ sites in total, each with a sample size $n=500$. 
In our main exposition, we consider an experimental study design where treatment propensity is $e(\bm{x}) = 0.5$, i.e., individuals are randomly assigned to treatment and control.
Variations of the settings above are discussed, with results presented in Appendix~\ref{suppl-sec:sim}. 



\textbf{Global heterogeneity: discrete and continuous.} 
Two types for global heterogeneity are considered by the choice of $U_k$. 
For \textbf{\emph{discrete grouping}}, we assume there are two underlying groups among the $K$ sites $U_k \sim Bernoulli(0.5)$. Specifically, we assume odd-index sites and even-index sites form two distinct groups $\mathcal{G}_1 = 
\{1, 3, \dots, K-1\}$; $\mathcal{G}_2 = 
\{2, 4\dots, K\}$ such that $U_{k\in\mathcal{G}_1} = 0$ and $U_{k\in\mathcal{G}_2} = 1$. 
Sites from similar underlying groupings have similar treatment effects and mean effects, while sites from different underlying groupings have different treatment effects and mean effects. 
For \textbf{\emph{continuous grouping}}, we consider $U_k \sim Unif[0,1]$. 
We vary the scales of the global heterogeneity under the discrete and continuous cases, respectively, with $c$ taking values $c \in \{0,0.6,1,2\}$. A $c=0$ implies all data sources are homogeneous. 
In other words, Assumption~2 is satisfied when $c = 0$ but not when $c > 0$. 





\textbf{Compared estimators and evaluation.} The proposed approaches ET and EF are compared with several competing methods. 
\textbf{LOC:} A local CT estimator
that does not utilize external information. 
It is trained on $\mathcal{I}_1$ only, combining training and estimation sets.
\textbf{MA:} A naive model averaging method
with weights ${\omega}_{k}^{\text{MA}} = 1/k$. This approach assumes models are homogeneous. 
\textbf{EWMA:}
We consider a modified version of EWMA that can be used for CATE.
We obtain an approximation of $\tau_1(\bm{x})$ by fitting another local model using the estimation set of site 1, denoted by $\widetilde\tau_1(\bm{x})$. 
Its weights are given by $${\omega}_{k}^{\text{EWMA}} = \frac{\exp\{- \sum_{i \in \mathcal{I}_1^{(2)}}(\widehat\tau_k(\bm{x}_i) - \widetilde\tau_1(\bm{x}_i))^2\} }{ \sum_{\ell=1}^{K} \exp\{- \sum_{i \in \mathcal{I}_1^{(2)}}(\widehat\tau_\ell(\bm{x}_i) - \widetilde\tau_1(\bm{x}_i))^2\} }.$$ 
\textbf{STACK:} A stacking ensemble, which is a linear ensemble of predictions of several models \citep{breiman1996stacked}.
To our end, we regress $\widetilde\tau_1(\bm{x})$ on the predictions of the estimation set in site 1 from each local model, $\{\widehat\tau_1(\bm{x}), \dots, \widehat\tau_k(\bm{x}) \}$. The stacking weights are not probabilistic hence not directly interpretable.
We report the empirical mean squared error (MSE) of these methods over an independent testing set of sample size $n_{te}=2000$ from site 1.
$\mbox{MSE}(\widehat{\tau}) = n_{te}^{-1}\sum_{i=1}^{n_{te}} \{ \widehat{\tau}(\bm{x}_i) - \tau_1(\bm{x}_i) \}^2. $
Each simulation scenario is repeated for 1000 times. 
Experiments are performed on a 6-core Intel Xeon CPU E5-2620 v3 2.40GHz equipped with 64GB RAM.


\begin{table}[tb]
\centering
\caption{MSE ratios of EF over LOC. As $n$ increases, model averaging becomes more powerful due to better estimation of $\tau_k$, and is more pronounced when $c$ is small.}
\label{fig:ratio}
\resizebox{0.9\columnwidth}{!}{
\begin{tabular}{@{}cccccc@{}}
\toprule
\multicolumn{1}{l}{}  &  & $c=0$  & $c=0.6$ & $c=1$  & $c=2$  \\ \midrule
\multirow{3}{*}{\shortstack{Discrete\\ grouping}} & $n=100$  & 0.57 & 0.59  & 0.61 & 0.59 \\
& $n=500$  & 0.12 & 0.17  & 0.17 & 0.16 \\
& $n=1000$ & 0.07 & 0.12  & 0.12 & 0.13 \\
\multirow{3}{*}{\shortstack{Continuous\\ grouping}} & $n=100$  & 0.54 & 0.59  & 0.63 & 0.69 \\
& $n=500$  & 0.11 & 0.24  & 0.31 & 0.34 \\
& $n=1000$ & 0.08 & 0.17  & 0.21 & 0.26 \\ \bottomrule
\end{tabular}
}
\vspace{-0.1cm}
\end{table}


\begin{figure*}[tb]
  \begin{subfigure}{0.3\textwidth}
    \centering
    \includegraphics[width=\linewidth]{figures/plt_sim_cf_ni500_et_2grp.pdf}
    \caption{}
  \end{subfigure}
  \begin{subfigure}{0.325\textwidth}
    \centering
    \includegraphics[width=\linewidth]{figures/plt_sim_cf_ni500_ef_2grp_pred.pdf}
    \caption{}
  \end{subfigure}
  \begin{subfigure}{0.37\textwidth}
    \centering
    \includegraphics[width=\linewidth]{figures/plt_sim_wts_sim_2grp_updateLegend.pdf}
    \caption{}
  \end{subfigure}
  \medskip

  \begin{subfigure}{0.3\textwidth}
    \centering
    \includegraphics[width=\linewidth]{figures/plt_sim_cf_ni500_et_cont.pdf}
    \caption{}
  \end{subfigure}
  \begin{subfigure}{0.325\textwidth}
    \centering
    \includegraphics[width=\linewidth]{figures/plt_sim_cf_ni500_ef_cont_pred.pdf}
    \caption{}
  \end{subfigure}
  \begin{subfigure}{0.37\textwidth}
    \centering
    \includegraphics[width=\linewidth]{figures/plt_sim_wts_sim_cont_updateLegend.pdf}
    \caption{}
  \end{subfigure}
  \vspace{-0.3cm}
  \caption{
 Visualization of simulation
 results under \textbf{discrete grouping (a,b,c)} and \textbf{continuous grouping (d,e,f)} when  $c=1$. (a) and (d) visualize the proposed ETs where the site indicator and $X_1$ are selected as splitting variables, which is consistent with the underlying data generation process. 
 (b) and (e) show the predicted treatment effects of the proposed EFs varying $X_1$ in each site, marginalized over all other features. (b) is arranged according to the true grouping, odd sites versus even sites. 
 The plot recovers the pattern of local and global heterogeneity. 
 (c) and (f) plot the interpretable model averaging weights in EFs over $X_1$. 
 The weights of site 1 have a relatively large contribution to the weighted estimator while models from other sites have different contributions for different $X_1$ depending on their similarity in $\tau(\bm{x},k)$ to that in site 1. 
 Corresponding ET and EF show consistent patterns and recover the true grouping.
  }
  \label{fig:sim_visual}
  \vspace{-0.3cm}
\end{figure*}




\textbf{Estimation performance}. 
Figure~\ref{fig:sim_box} shows the performance of the proposed estimators and the competing estimators, using LOC as the benchmark.
The proposed ET and EF show the best performance in terms of the mean and variation of MSE among other estimators when $c > 0$, and comparable to equal weighting MA when $c = 0$.
Although, a forest is more stable than a tree in practice,
both ET and EF give similar results because the true model is relatively simple and can be accurately estimated by a single ensemble tree under the given sample size. 


Although asymptotically consistency, under finite sample, bias exists in local models and leads to biased model averaging estimates. 
While explicit quantification of bias and variance remains challenging due to extra uncertainty carried forward from the local estimates, we demonstrated that the proposed estimators can improve upon the local models under small sample size via Table~\ref{fig:ratio}.
It shows the MSE ratio of EF over LOC as a measure of gain resulting from model averaging by varying $n =100, 500, 1000$. 
The decrease in MSE ratio as $n$ increases, regardless of the choice of $c$, is consistent with our asymptotic results in Theorem~\ref{Theorem}.
This is due to a bias-and-variance trade-off in the ensemble that ensures a small MSE, which remains smaller than that in LOC despite varying $n$.
It also shows our method is robust to the existence of local uncertainty.







\textbf{Visualization of information borrowing.} Figure~\ref{fig:sim_visual} visualizes the proposed ET and EF. 
In (a) and (d), the site indicator and $X_1$ appear as splitting variables in the ETs, which is consistent with the data generation process. 
The estimated treatment effect (b) and (e) reveals the pattern of transportability across sites and with respect to $X_1$. 
Panels (c) and (f) plot the model averaging weights in EFs over $X_1$. Site 1 has a relatively large contribution to the weighted estimator while models from other sites have different contributions at different values of $X_1$ depending on their similarity in $\tau(\bm{x},k)$ to that in site 1. 
Corresponding ET and EF show consistent patterns. 


\begin{figure*}[htp
\centering
 \begin{subfigure}{0.2\textwidth}
  \centerline{\includegraphics[width=\linewidth]{figures/plt_real_reward.pdf}}
  \caption{}
 \end{subfigure}
 \begin{subfigure}{0.39\textwidth}
  \centerline{\includegraphics[width=\linewidth]{figures/plt_real_pdp_ef_obs.pdf}}
  \caption{}
 \end{subfigure}
 \begin{subfigure}{0.39\textwidth}
  \centerline{\includegraphics[width=\linewidth]{figures/plt_real_bmi_wt_updatedLegend.pdf}}
  \caption{}
 \end{subfigure}
 \vspace{-0.3cm}
 \caption{
 Application to estimating treatment effects of oxygen therapy on survival. 
 (a) Expected survival of treatment decision following different estimators. 
 The proposed EF shows the largest gain in improving survival rate, more promising than LOC and baseline. 
 (b) Estimated treatment effects varying duration and BMI, two important features in the fitted EF. 
 Patients with a BMI around 35, and a duration above 400 benefited the most. 
 (c) Visualization of data-adaptive weights in the estimated EF varying BMI. 

 Hospitals with a larger bed capacity tend to contribute more, the data of which might be more similar to hospital 1.
}
\vspace{-0.2cm}
\label{fig:real}
\end{figure*}


\textbf{Additional simulations.}
The detailed results of these additional simulations are included in Appendix~\ref{suppl-sec:sim}.
\emph{\textbf{1) Connection to supervised learning.}} The uniqueness of averaging $\tau_k(\bm{x})$ as opposed to supervised learning that averages prediction models $f_k(\bm{x})$ is that the outcome of $f_k(\bm{x})$ is immediately available. 
In our case, an additional estimation step is needed to construct the model averaging weights. 
We provide a comparison among estimators that utilize the ground truth $\{\tau_k(\bm{x})\}_{k=1}^K$ (denoted as ``-oracle'') when computing ensemble weights.
This mimics the case of supervised learning where weights are based on observed outcomes. 
Oracle methods achieve smaller MSE ratios; the pattern is consistent with Table~\ref{fig:ratio}.
\emph{\textbf{2) Simulation under observational studies.}} 
We also consider the treatment generation mechanism under an observational design. 
Specifically, the propensity is given as $e(\bm{x}) = \text{expit}(0.6x_1)$. We consider both a correctly specified propensity model using a logistic regression of $Z$ on $X_1$ and a misspecified propensity model with a logistic regression of $Z$ on all $\bm{X}$. In general, the proposed estimators obtain the best performance with similar results as in Figure~\ref{fig:sim_box}. With the correctly specified propensity score model, the local estimator is consistent in estimating $\tau_k(\bm{x})$, the proposed framework is valid. When the propensity model is misspecified, extra uncertainty is carried forward from the local estimates, but the proposed estimators can still improve upon LOC. This is due to a bias-and-variance trade-off that leads to small MSE, which remains smaller than the local models. 
\emph{\textbf{3) Covariate dimensions.}} 
Besides $D=5$, we consider other choices of covariate dimension including $D=20, 50$. With a higher dimension, the MSE ratio between the proposed estimates and LOC estimates increases but the same pattern across methods persists.
\emph{\textbf{4) Unequal sample size at each site.}} 
In the distributed date network, different sites may have a different sample size $n_k$. Those with a smaller sample size may not be representative of their population, leading to an uneven level of precision for local causal estimates. We consider a simulation setting where site 1 has a sample size of $n_1=500$ while other site $n_2,\ldots,n_K$ has a sample size of 200. Results show that the MSE ratio between the proposed estimates and LOC estimates increases compared to the scenario where the sample size in all sites are 500. However, the proposed estimators still enjoy the most robust performance. 
This also shows our method is robust to the existence of local uncertainty. 
\emph{\textbf{5) Different local estimators.}} We stress that other consistent estimators could be used as the local model. Options such as causal forest \citep{wager2018estimation} are explored varying the sample size at local sites. Similar performance is observed as in Figure~\ref{fig:sim_box}.
\emph{\textbf{6) Further comparisons to non-adaptive ensemble.}} Here we provide a brief discussion of the implications of the proposed method and how it differs from non-adaptive methods such as stacking. 
Although unrealistic, when the true weights are non-adaptive, the performance may be similar. Plus, our learned weights can be used to examine adaptivity, as shown in Figure~\ref{fig:sim_visual}(c,f) and Figure~\ref{fig:real}(c). 
Stacking is shown to be more robust than non-adaptive model averaging in case of model misspecification. See discussion in \citet{clarke2003comparing}. 
Our additional simulation results show that in case of a large global heterogeneity, as $c$ increases, the heterogeneity across sites gets larger, reducing the influence of important covariates on heterogeneity, hence the weights become more non-adaptive. However, the proposed methods still enjoy a comparable performance to STACK, which further indicates the robustness of the proposed methods. 






   





\section{Example: a Multi-Hospital Data Network}
\label{sec:application}












Application with contextual insights is provided based on an analysis of the eICU Collaborative Research Database, a multi-hospital database published by Philips Healthcare \citep{pollard2018eicu}. 
The analysis is motivated by a recent retrospective study that there is a higher survival rate when SpO$_2$ is maintained at 94-98\% among patients requiring oxygen therapy \citep{van2020search}, not ``the higher the better''. 
We use the same data extraction code to create our data. 
We consider SpO$_2$ within this range as treatment ($Z=1$) and outside of this range as control ($Z=0$). 
A total of 7,022 patients from 20 hospitals, each with at least 50 patients in each treatment arm, are included with a randomly selected target (hospital 1).
Hospital-level summary information is provided in Appendix~\ref{suppl-sec:real}.
Patient-level features include age, BMI, sex, Sequential Organ Failure Assessment (SOFA) score, and duration of oxygen therapy.
The outcome is hospital survival ($Y=1$) or death ($Y=0$).










Figure~\ref{fig:real} visualizes the performance of EF-based estimated effect of oxygen therapy setting on in-hospital survival. 
CT is used as the local model with propensity score modeled by a logistic regression. 
Figure~\ref{fig:real}(a) shows the propensity score-weighted \textbf{average survival} for those whose received treatment is consistent with the estimated decision. 
Specifically, the expected reward is given by $$\frac{\sum_i Y_i 1(Z_i = Z^{est}_i) / \pi(Z_i,\bm{X}_i)}{\sum_i 1(Z_i = Z^{est}_i) / \pi(Z_i,\bm{X}_i)},$$ where $Z^{est}_i = 1(\widehat\tau>0)$ denotes the estimated treatment rule and $\pi(Z_i,\bm{X}_i)$ is the probability of receiving the actual treatment.
We provide expected reward for the 1) observed treatment assignment (baseline), 2) LOC-based rule, and 3) EF-based rule. 
The treatment rule based on our method can increase mean survival by 3\% points compared to baseline, and is more promising than LOC. 

In the fitted EF, 
the hospital indicator is the most important, explaining about 50\% of the decrease in training error. 
Figure~\ref{fig:real}(b) shows the estimated CATE varying two important features, BMI and oxygen therapy duration.
Patients with BMI around 36 and duration above 400 show the most benefit from oxygen therapy in the target SpO$_2$ range. 
Patients with BMI between 20 and 30 and duration around 200 may not benefit from such alteration. 
Figure~\ref{fig:real}(c) visualizes 
the
data-adaptive weights $\omega_{k}(\bm{x})$ in the fitted EF with respect to BMI for different models, while holding other variables constant. The weights of hospital 1 are quite stable while models from other sites may have different contribution to the weighted estimator for different values of BMI. 
Judging from hospital information 
in Appendix~\ref{suppl-sec:real}, hospitals with a larger bed capacity tend to be similar to hospital 1, and are shown to provide larger contributions.





In this distributed research network, different hospitals have a different sample size.
For sensitivity analysis, we consider a weighting strategy to adjust for the sample size of site $k$. 
Results show similar patterns as in Figure \ref{fig:real}. Detailed results 
are provided in Appendix~\ref{suppl-sec:real}. 
The real-data access is provided in Appendix~\ref{suppl-sec:code}. 















\section{Discussion} 
\label{sec:disc}


We have proposed an efficient and interpretable tree-based model averaging framework for enhancing treatment effect estimation at a target site by borrowing information from potentially heterogeneous data sources. We generalize standard model averaging scheme in a data-adaptive way such that the generated weights depend on subject-level features. 
This work makes multi-site collaborations and especially treatment effect estimation more practical by avoiding the need to share subject-level data. 
Our approach extends beyond causal inference to estimating a general $f(\bm{x})$ from heterogeneous data. 



Unlike in classic model averaging where prediction performance can be assessed against observed outcomes or labels, treatment effects are not directly observed. 
While our approach is guaranteed to be consistent under randomized studies, 
the weights are estimated based on expected treatment effects, hence relying on  Assumption 1 (unconfoundedness) to hold. It may be a strong assumption in observational studies with unmeasured confounding. 





\section*{Acknowledgments}
This research was supported in part by the Competitive Medical Research Fund of the UPMC Health System and the University of Pittsburgh Center for Research Computing through the resources provided. 
The authors thank Gong Tang and Yu Wang for insightful discussion.




