# Multiply Robust Federated Estimation of Targeted Average Treatment Effects

Larry Han
Department of Health Sciences
Northeastern University
Boston, MA 02115
<EMAIL_ADDRESS>
&Zhu Shen
Department of Biostatistics
Harvard University
Boston, MA 02115
<EMAIL_ADDRESS>
Jose Zubizarreta
Departments of Health Care Policy, Biostatistics, and Statistics
Harvard University
Boston, MA 02115
<EMAIL_ADDRESS>

###### Abstract

Federated or multi-site studies have distinct advantages over single-site
studies, including increased generalizability, the ability to study
underrepresented populations, and the opportunity to study rare exposures and
outcomes. However, these studies are challenging due to the need to preserve
the privacy of each individual’s data and the heterogeneity in their covariate
distributions. We propose a novel federated approach to derive valid causal
inferences for a target population using multi-site data. We adjust for
covariate shift and covariate mismatch between sites by developing multiply-
robust and privacy-preserving nuisance function estimation. Our methodology
incorporates transfer learning to estimate ensemble weights to combine
information from source sites. We show that these learned weights are
efficient and optimal under different scenarios. We showcase the finite sample
advantages of our approach in terms of efficiency and robustness compared to
existing approaches.

##  1 Introduction

Compared to single-site studies, federated or multi-site studies confer
distinct advantages, such as the potential for increased generalizability of
findings, the opportunity to learn about underrepresented populations, and the
ability to study rare exposures and outcomes. However, deriving valid causal
inferences using multi-site data is difficult due to numerous real-world
challenges, including heterogeneity of site populations, different data
structures, and privacy-preserving constraints stemming from policies such as
the General Data Protection Regulation (GDPR) and Health Insurance Portability
and Accountability Act (HIPAA) that prohibit direct data pooling.

Recent methodological developments have focused on privacy-preserving
estimation strategies. These strategies typically involve sharing summary-
level information from multiple data sources [26, 27, 12, 13, 18]. However,
they often require restrictive assumptions such as homogeneous data structures
and model specifications (e.g., a common set of observed covariates measured
using a common data model), which are not realistic in practice.

To address these methodological gaps, we propose a multiply robust and
privacy-preserving estimator that leverages multi-site information to estimate
causal effects in a target population of interest. Compared to existing
approaches, our method allows investigators from different sites to
incorporate site-specific covariate information and domain knowledge and
provides increased protection against model misspecification. Our method
allows for flexible identification under different settings, including
systematically missing covariates and different site-specific covariates
(termed covariate mismatch). Our proposed method adopts an adaptive ensembling
approach that optimally combines estimates from source sites and serves as a
data-driven metric for the transportability of source sites. Moreover, the
proposed method relaxes the assumption of homogeneous model specifications by
adopting a class of multiply robust estimators for estimating the nuisance
functions.

###  1.1 Related Work and Contributions

The current literature on multi-site causal inference typically assumes that a
common set of confounders is observed in all sites [7, 6, 12, 13]. However,
this assumption is rarely met due to variations in local practices, e.g.,
differing data collection standards and coding practices. In particular, the
target site often lacks data on certain covariates available in the source
sites, and ignoring them can result in biased and inefficient inference [31].
Recently, [28] proposed a method to address covariate mismatch by integrating
source samples with unmeasured confounders and a target sample containing
information about these confounders. However, they assume that the target
sample is either a simple random sample or is sampled according to a known
procedure from the source population. [11] extended the method to the setting
where the average treatment effect (ATE) is not identifiable in some sites by
constructing control variates. However, their approach is limited to
addressing selection biases in case-control studies and cannot be easily
extended to other outcome types. [31] extended the framework by [7, 6] to
handle covariate mismatch by regressing predicted conditional outcomes on
effect modifiers and then taking the means of these regression models
evaluated on target site samples. Our work extends [31] to the multi-site and
federated data setting by utilizing an adaptive weighting approach that
optimally combines estimates from source sites.

Most existing approaches in the generalizability and transportability
literature deal with heterogeneous covariate distributions by modeling site
selection processes [5, 1, 25, 7, 6]. For instance, covariate shift can be
accounted for via inverse probability of selection weights [5, 1],
stratification [25], or augmentation [7, 6], which require pooling individual-
level information across sites. Our work differs from those in that we
preserve individual data privacy, sharing only summary-level information about
the target site. Specifically, we adopt density ratio models [22, 8] that only
share covariate moments of the target samples. Under certain specifications,
these density ratio models are equivalent to logistic regression-based
selection models for adjusting heterogeneity between target and source
populations. Our approach shares similarities with calibration weighting
methods, but we employ semi-parametric efficiency theory to enable a closed-
form approximation of variance rather than relying on the bootstrap.

Further, when data sources are heterogeneous, it would be beneficial for
investigators at different sites to incorporate site-specific knowledge when
specifying candidate models. However, to the best of our knowledge, existing
methods require common models to be specified across sites, which may not be
realistic or flexible enough [26, 27, 12, 13, 18]. We relax this requirement
by adopting a multiply robust estimator, allowing investigators in each site
to propose multiple, different outcome and treatment models. The estimator is
consistent if any one of the multiple outcome or treatment models is correctly
specified. Our work builds on [20], which established an equivalence between
doubly robust and multiply robust estimators using mixing weights determined
by predictive risks of candidate models [16, 14, 15, 2, 3, 4]. Further, to
avoid negative transfer due to non-transportable source estimates, we adopt a
data-adaptive ensembling approach [12, 13] that guarantees that the federated
estimator achieves improved precision compared to an estimator using target
site data alone when at least one source estimate is sufficiently similar to
the target estimate.

##  2 Preliminaries

We consider data from $K$ sites, where each site has access to its individual-
level data but is prohibited from sharing this data with any other site. The
set of sites will be denoted by $\mathcal{K}=\\{1,2,...,K\\}$. Without loss of
generality, we define the target site to be the first site, i.e., $T=\\{1\\}$
and the source sites as the remaining sites, i.e.,
$\mathcal{S}=\mathcal{K}\setminus T=\\{2,...,K\\}$.

For individual $i$, let $Y_{i}$ denote an observed outcome, which can be
continuous or discrete. $X_{i}\in\mathbb{R}^{p}$ represents the
$p$-dimensional baseline covariates in source site $k\in\mathcal{S}$.
$V_{i}\in\mathbb{R}^{q}$ represents the (partial) baseline covariates in the
target site $T$ such that $V_{i}\subseteq X_{i}$. To simplify the
presentation, we assume an identical set of covariates across all source
sites, although our method can accommodate scenarios where distinct covariate
sets are present among the source sites. Let $A_{i}$ represent a binary
treatment indicator, with $A_{i}=1$ denoting treatment and $A_{i}=0$ denoting
control. $R_{i}$ is a site indicator with $R_{i}=k$ if patient $i$ is from the
site $k$. We observe $n_{T}$ target samples,
$D_{T}=\\{Y_{i},V_{i},A_{i},R_{i}=T,1\leq i\leq n_{T}\\}$ and $n_{k}$ source
samples, $D_{k}=\\{Y_{i},X_{i},A_{i},R_{i}=k,1\leq i\leq n_{k}\\}$ for each
$k\in\mathcal{S}$. The total sample size is $N=\sum_{k\in\mathcal{K}}n_{k}$.
Under the potential outcomes framework [21, 23], we denote the counterfactual
outcomes under treatment and control as $\left\\{Y_{i}(1),Y_{i}(0)\right\\}$,
and only one of them is observed: $Y_{i}=A_{i}Y_{i}(1)+(1-A_{i})Y_{i}(0)$
[24]. The data structure is illustrated in Figure 1.

Figure 1: Schematic of the data structure in the multi-site setting.

Our goal is to estimate the average treatment effect in the target population,

$\Delta_{T}=\mu_{1,T}-\mu_{0,T}\quad\text{ where
}\quad\mu_{a,T}={E}\left\\{Y_{i}{(a)}\mid R_{i}=T\right\\}\text{ for
}a\in\\{0,1\\},$ (1)

where $\mu_{a,T}$ is the mean potential outcome under treatment $a$ in the
target population. To identify this quantity, we consider the following
assumptions:

  1. (A1)

(Consistency): For every individual $i$, if $A_{i}=a$, then
$Y_{i}=Y_{i}{(a)}$.

  2. (A2)

(Mean exchangeability over treatment assignment in the target population):
$\\\
{E}\\{Y_{i}{(a)}\mid{V}_{i}=v,A_{i},R_{i}=T\\}={E}\left\\{Y_{i}{(a)}\mid{V}_{i}=v,R_{i}=T\right\\}$.

  3. (A3)

(Positivity of treatment assignment in the target population): $\\\
0<P(A_{i}=1\mid{V}_{i}=v,R_{i}=T)<1$ for any $v$ s.t. $P({V}_{i}=v\mid
R_{i}=T)>0$.

  4. (A4)

(Mean exchangeability over treatment assignment in the source populations):
$\\\
{E}\\{Y_{i}{(a)}\mid{X}_{i}=x,A_{i},R_{i}=k\\}={E}\left\\{Y_{i}{(a)}\mid{X}_{i}=x,R_{i}=k\right\\}$,
$k\in\mathcal{S}$.

  5. (A5)

(Positivity of treatment assignment in the source populations): $\\\
0<P(A_{i}=1\mid{X}_{i}=x,R_{i}=k)<1$ for any $x$ s.t. $P({X}_{i}=x\mid
R_{i}=k)>0$, $k\in\mathcal{S}$.

  6. (A6)

(Mean exchangeability over site selection): $\\\ {E}\\{Y_{i}{(a)}\mid
V_{i}=v,R_{i}=k\\}={E}\left\\{Y_{i}{(a)}\mid V_{i}=v\right\\}$,
$k\in\mathcal{K}$.

  7. (A7)

(Positivity of site selection): $\\\ 0<P(R_{i}=k\mid{V}_{i}=v)<1$ for
$k\in\mathcal{S}$ and any $v$ s.t. $P({V}_{i}=v)>0$.

Assumption (A1) is the stable unit treatment value assumption (SUTVA),
requiring no interference between individuals. Assumption (A2) (Assumption
(A4)) states that the mean counterfactual outcome under treatment $a$ is
independent of treatment assignment, conditional on baseline covariates in the
target (source) populations. For Assumption (A2) and (A4) to hold, we require
all effect modifiers to be measured in $V$. Assumption (A3) (Assumption (A5))
states that each individual in the target (source) populations has a positive
probability of receiving each treatment. Assumption (A6) states that the mean
counterfactual outcome is independent of site selection, conditional on
covariates in the target population. For Assumption (A6) to hold, we require
all covariates that are distributed differently between target and source
populations (shifted covariates) to be measured in $V$. Thus, if these effect
modifiers are measured in $V$, Assumption (A2), (A4) and (A6) automatically
hold. Assumption (A7) requires that in each stratum defined by $V$, the
probability of being in a source population for each individual is positive.
Theorem 1 shows that under Assumption (A1), (A4) \- (A7), the mean
counterfactual outcome for the target can be identified in the sources. Since
these assumptions may not hold in practice, we devise a data-adaptive
ensembling procedure in Section 4 to screen out sites that significantly
violate these assumptions.

######  Theorem 1.

If Assumptions (A1) \- (A3) hold, the mean counterfactual outcomes in the
target population can be identified using the target sample.

$\mu_{a,T}=E\left\\{Y_{i}{(a)}\mid
R_{i}=T\right\\}=E\left\\{E\left\\{Y_{i}\mid
V_{i}=v,A_{i}=a,R_{i}=T\right\\}\mid R_{i}=T\right\\}.$ (2)

If Assumptions (A1), (A4) \- (A7) hold, the mean counterfactual outcomes in
the target population can be identified using the source samples.

$\displaystyle\mu_{a,T}$ $\displaystyle=E\left\\{Y_{i}{(a)}\mid
R_{i}=T\right\\}$ (3)
$\displaystyle=E\left\\{E\left\\{E\left\\{Y_{i}\mid{X}_{i}=x,A_{i}=a,R_{i}=k\right\\}\mid
V_{i}=v,R_{i}=k\right\\}\mid R_{i}=T\right\\}.$

##  3 Site-specific Estimators

For the target site $k=\\{T\\}$, a standard AIPW estimator is used for
$\mu_{a,T}$ as follows

$\widehat{\mu}_{a,T}=\frac{1}{n_{T}}\sum_{i=1}^{n}\biggl{[}\frac{I(A_{i}=a,R_{i}=T)}{\widehat{\pi}_{a,T}(V_{i})}\Bigl{\\{}Y_{i}-\widehat{m}_{a,T}(V_{i})\Bigr{\\}}+\widehat{m}_{a,T}(V_{i})\biggr{]},$
(4)

where $\widehat{m}_{a,T}(V_{i})$ is an estimator for $E\left\\{Y_{i}\mid
V_{i}=v,A_{i}=a,R_{i}=T\right\\}$, the outcome model in the target population,
and $\widehat{\pi}_{a,T}(V_{i})$ is an estimator for $P(A_{i}=1\mid
V_{i}=v,R_{i}=T)$, the probability of receiving treatment $a$ in the target
population.

For each source site $k\in\mathcal{S}$, we propose an estimator for
$\mu_{a,T}$ as follows

$\displaystyle\widehat{\mu}_{a,k}=$
$\displaystyle\frac{1}{n_{k}}\sum_{i=1}^{n}\biggl{[}\frac{I(A_{i}=a,R_{i}=k)}{\widehat{\pi}_{a,k}(X_{i})}\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}Y_{i}-\widehat{m}_{a,k}(X_{i})\Bigr{\\}}\biggr{]}$
(5) $\displaystyle+$
$\displaystyle\frac{1}{n_{k}}\sum_{i=1}^{n}\biggl{[}I(R_{i}=k)\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}\widehat{m}_{a,k}(X_{i})-\widehat{\tau}_{a,k}(V_{i})\Bigr{\\}}\biggr{]}+\frac{1}{n_{T}}\sum_{i=1}^{n}I(R_{i}=T)\widehat{\tau}_{a,k}(V_{i}),$

where $\widehat{\tau}_{a,k}(V_{i})$ is an estimator for
$E\left\\{m_{a,k}(x)\mid V_{i}=v,R_{i}=k\right\\}$ and
$\widehat{m}_{a,k}(X_{i})$ is an estimator for $E\left\\{Y_{i}\mid
X_{i}=x,A_{i}=a,R_{i}=k\right\\}$. $\widehat{\zeta}_{k}(V_{i})$ estimates
$f(V_{i}\mid R_{i}=T)/f(V_{i}\mid R_{i}=k)$, the density ratios of covariate
distributions in the target population $T$ and source population
$k\in\mathcal{S}$. $\widehat{\pi}_{a,k}(X_{i})$ estimates $P(A_{i}=1\mid
X_{i}=x,R_{i}=k)$, the probability of receiving treatment $a$ in source
$k\in\mathcal{S}$.

Compared to the transportation estimators in [7, 6], we introduce two
additional nuisance functions, $\zeta_{k}(V_{i})$ and ${\tau}_{a,k}(V_{i})$.
Specifically, $\zeta_{k}(V_{i})$ accounts for covariate shift across sites,
while ${\tau}_{a,k}(V_{i})$ is introduced to address covariate mismatch across
sites. We provide estimation procedures for these nuisance functions in the
following subsections, and the theoretical guarantees of the estimator are
presented in Section 5.

###  3.1 Density Ratio Weighting

Most existing methods for adjusting for heterogeneity of site populations rely
on inverse probability of selection weighting, which requires pooling target
and source samples. However, such pooling is often restricted to protect
individuals’ data privacy. We propose a density ratio weighting approach,
which offers equivalent estimation without the need for direct data pooling
(see Appendix A).

Formally, we model the density ratios of covariate distributions in the target
$T$ and source $k\in\mathcal{S}$ by specifying an exponential tilt model [22,
8]; ${\zeta}_{k}(V_{i};\gamma_{k})={f(V_{i}\mid R_{i}=T)}/{f(V_{i}\mid
R_{i}=k)}=\exp\left\\{-{\gamma}_{k}^{\top}\psi({V_{i}})\right\\}$ where
$f(V_{i}\mid R_{i}=T)$ and $f(V_{i}\mid R_{i}=k)$ are density functions of
covariates $V_{i}$ in the target $T$ and source $k\in\mathcal{S}$,
respectively, and $\psi({V_{i}})$ is some $d$-dimensional basis with $1$ as
its first element. With this formulation, ${\zeta}_{k}(V_{i};\gamma_{k})=1$
for $\gamma_{k}=0$ and $\int{\zeta}_{k}(V_{i};\gamma_{k}){f(V_{i}\mid
R_{i}=k)}dx=1$. If we choose $\psi({V_{i}})=V_{i}$, we can recover the entire
class of natural exponential family distributions. If we include higher-order
terms, the exponential tilt model has greater flexibility in characterizing
the heterogeneity between two populations [9]. We solve for
$\widehat{\gamma}_{k}$ with the following estimating equation

$\frac{1}{n_{T}}\sum_{i=1}^{N}I\left(R_{i}=T\right)\psi\left(V_{i}\right)=\frac{1}{n_{k}}\sum_{i=1}^{N}I\left(R_{i}=k\right)\psi\left(V_{i}\right)\exp\left\\{-{\gamma}_{k}^{\top}\psi({V_{i}})\right\\}.$
(6)

This procedure preserves individual privacy; choosing $\psi({V_{i}})=V_{i}$,
the target site only needs to share its covariate means with the source sites;
each source site then solves (6) with its own data to obtain the density
ratios.

###  3.2 Multiply Robust Estimator

We relax the assumption of homogeneous model specifications across sites and
allow each site to propose multiple models for nuisance functions. Our
proposal follows the construction of multiply robust estimators for nuisance
functions via a model-mixing approach [20].

Formally, for each site $k\in\mathcal{K}$, we consider a set of $J$ candidate
treatment models for the propensity scores
$\\{\pi^{j}_{a,k}\left(x\right):j\in\mathcal{J}=\\{1,...,J\\}\\}$. Let
$\widehat{\pi}_{a,k}^{j}(x)$ be the estimator of $\pi^{j}_{a,k}(x)$ obtained
by fitting the corresponding candidate models on the data, which can be
parametric, semiparametric, or nonparametric machine learning models.
$\widehat{\pi}_{a,k}(X_{i})=\sum_{j=1}^{J}\widehat{\Lambda}_{j}\widehat{\pi}^{j}_{a,k}(X_{i})$
denotes the weighted predictions of propensity scores, with weights
$\widehat{\Lambda}_{j}$ assigned to predictions by each candidate model $j$.
To calculate the weights $\widehat{\Lambda}_{j}$, we adapt a model-mixing
algorithm developed in [29] and [20] based on the cumulative predictive risks
of candidate models.

First, we randomly partition the data within each site into a training set
$D^{\text{train}}_{k}$ of units indexed by $\\{1,...,n_{k}^{\text{train}}\\}$
and a validation set $D^{\text{val}}_{k}$ of units indexed by
$\\{n_{k}^{\text{train}}+1,...,n_{k}\\}$. Then, each candidate treatment model
is fit on $D^{\text{train}}_{k}$ to obtain
$\widehat{\pi}_{a,n_{k}^{\text{train}}}^{j}$ for $j\in\mathcal{J}$. The model-
mixing weights are determined by the models’ predictive risks assessed on
$D^{\text{val}}_{k}$ according to the Bernoulli likelihood. Specifically,

$\displaystyle\widehat{\Lambda}_{j}$
$\displaystyle=\left(n_{k}-n_{k}^{\text{train}}\right)^{-1}\sum_{i=n_{k}^{\text{train}}+1}^{n_{k}}\widehat{\Lambda}_{j,i}\quad\text{and}$
(7) $\displaystyle\widehat{\Lambda}_{j,i}$
$\displaystyle=\frac{\Pi_{q=n_{k}^{\text{train}}+1}^{i-1}\widehat{\pi}_{a,n_{k}^{\text{train}}}^{j}\left(X_{q}\right)^{A_{q}}\left\\{1-\widehat{\pi}_{a,n_{k}^{\text{train}}}^{j}\left(X_{q}\right)\right\\}^{1-A_{q}}}{\sum_{j^{\prime}=1}^{J}\Pi_{q=n_{k}^{\text{train}}+1}^{i-1}\widehat{\pi}_{a,n_{k}^{\text{train}}}^{j^{\prime}}\left(X_{q}\right)^{A_{q}}\left\\{1-\widehat{\pi}_{a,n_{k}^{\text{train}}}^{j^{\prime}}\left(X_{q}\right)\right\\}^{1-A_{q}}}\quad\text{for}\quad
n_{k}^{\text{train}}+2\leq i\leq n_{k},$

where $\widehat{\Lambda}_{j,n_{k}^{\text{train}}+1}=1/J$. The model mixing
estimators are consistent if one of the $j\in\mathcal{J}$ candidate models is
correctly specified [20]. A similar strategy extends for conditional outcomes
$m_{a,k}(X_{i})$ by combining a set of $L$ candidate outcome models
$\\{m^{l}_{a,k}\left(x\right):l\in\mathcal{L}=\\{1,...,L\\}\\}$. We obtain
$\widehat{m}_{a,k}(X_{i})=\sum_{l=1}^{L}\widehat{\Omega}_{l}\widehat{m}_{a,k}^{l}(X_{i})$
as the predicted outcomes with weights $\widehat{\Omega}_{l}$ of candidate
outcomes models under treatment $a$ in site $k$. Further details are provided
in Appendix B.

###  3.3 Handling Covariate Mismatch

To account for covariate mismatch, we adapt the approach in [31], introducing
the nuisance function $\tau_{a,k}(V_{i})=E\\{m_{a,k}(x)\mid
V_{i}=v,R_{i}=k\\}$, where $m_{a,k}(x)$ is the outcome regression for
treatment $a$ in site $k$. First, we estimate $m_{a,k}(X_{i})$ by regressing
the outcome $Y_{i}$ on covariates $X_{i}$ among units receiving treatment $a$
in site $k$. We then regress $\widehat{m}_{a,k}(X_{i})$, the estimates from
the previous step, on $V_{i}$ in the source site $k$ to obtain
$\widehat{\tau}_{a,k}(x)$. By doing so, we project all site-specific estimates
of conditional outcomes to a common hyperplane defined by $V_{i}$. If all
effect modifiers that are distributed differently between target and source
populations are measured in $V_{i}$, then the information contained in the
projected site-specific estimates can be transported to the target site.
Finally, we take the mean of $\widehat{\tau}_{a,k}(x)$ over the target sample,
which gives us the transported estimate $\hat{\tau}_{a,k}(V_{i})$ for the mean
counterfactual outcomes under treatment $a$ in the target population.

##  4 Federated Global Estimator

Let $\widehat{\mu}_{a,T}$ denote the estimate of $\mu_{a,T}$ based on target
data only and $\widehat{\mu}_{a,k}$ be the estimates of $\mu_{a,T}$ using
source data $k\in\mathcal{S}$. We propose a general form of the federated
global estimator as follows

$\widehat{\mu}_{a,G}=\widehat{\mu}_{a,T}+\sum_{k\in\mathcal{K}}\widehat{\eta}_{k}\left\\{\widehat{\mu}_{a,k}-\widehat{\mu}_{a,T}\right\\},$
(8)

where $\widehat{\eta}_{k}\geq 0$ is a non-negative weight assigned to site-
specific estimates and $\sum_{k\in\mathcal{K}}\widehat{\eta}_{k}=1$. The role
of $\eta_{k}$ is to determine the ensemble weight given to the site-specific
estimates. We can employ diverse weighting methods by selecting appropriate
values of $\eta_{k}$. For example, if $\eta_{k}=0$, the global estimator is
simply the estimator based on target data only; if $\eta_{k}=n_{k}/N$, the
global estimator combines site-specific estimates by their sample sizes; if
$\eta_{k}=(1/\sigma_{k}^{2})/\sum_{j\in\mathcal{K}}(1/\sigma_{j}^{2})$ where
$\sigma_{k}^{2}=\text{Var}(\widehat{\mu}_{a,k})$, the global estimator is the
inverse variance weighting estimator, which is known to be appropriate when
working models are homogeneous across sites [27]. To control for bias due to
non-transportable site estimates while achieving optimal efficiency, we
estimate $\eta_{k}$ data-adaptively by a penalized regression of site-specific
influence functions [12, 13]. This strategy ensembles the site-specific
estimates for higher efficiency if they are sufficiently similar to the target
estimates; if source estimates are significantly different, their weights will
be shrunk toward zero with high probability.

We denote the data-adaptive weights as $\eta_{k,L_{1}}$, obtained as the
solutions to a penalized regression of the site-specific influence functions
as follows

$\widehat{\eta}_{k,L_{1}}=\arg\min_{\eta_{k}\geq
0}\sum_{i=1}^{N}\left[\widehat{\xi}_{T,i}{(a)}-\sum_{k\in\mathcal{K}}\eta_{k}\left(\widehat{\xi}_{T,i}{(a)}-\widehat{\xi}_{k,i}{(a)}-\widehat{\delta}_{k}\right)\right]^{2}+\lambda\sum_{k\in\mathcal{K}}\left|\eta_{k}\right|\widehat{\delta}_{k}^{2},$
(9)

where $\widehat{\xi}_{T,i}{(a)}$ and $\widehat{\xi}_{k,i}{(a)}$ are the
estimated influence functions for the target and source site estimators (see
Appendix D.3 for the exact form of the influence functions). The estimated
difference $\widehat{\delta}_{k}=\widehat{\mu}_{a,k}-\widehat{\mu}_{a,T}$
quantifies the bias between the estimate from source $k\in\mathcal{S}$ and the
estimate from the target $T$. The tuning parameter $\lambda$ determines the
penalty imposed on source site estimates and in practice, is chosen via cross-
validation. Specifically, we create a grid of values of $\lambda$ and
iteratively train and evaluate the model using different $\lambda$ values,
selecting the one with the lowest average validation error after multiple
sample splits.

We estimate the variance of $\widehat{\mu}_{a,G}$ using the estimated
influence functions for $\widehat{\mu}_{a,T}$ and $\widehat{\mu}_{a,k}$. By
the central limit theorem,
$\sqrt{N}(\widehat{\mu}_{a,G}-\bar{\mu}_{a,G})\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,\Sigma)$, where
$\Sigma=E\\{\sum_{k\in\mathcal{K}}\bar{\eta}_{k}\xi_{k,i}(a)\\}^{2}$ and
$\bar{\mu}_{a,G}$ and $\bar{\eta}_{k}$ denote the limiting values of
$\widehat{\mu}_{a,G}$ and $\hat{\eta}_{k}$ respectively. The standard error of
$\widehat{\mu}_{a,G}$ is estimated as $\sqrt{\widehat{\Sigma}/N}$ where
$\widehat{\Sigma}=N^{-1}\sum_{k\in\mathcal{K}}\sum_{i=1}^{n_{k}}\\{\widehat{\eta}_{k}\widehat{\xi}_{k,i}{(a)}\\}^{2}$.
A two-sided $(1-\alpha)\times$100% confidence interval for $\mu_{a,G}$ is

$\widehat{\mathcal{C}}_{\alpha}=\left[\widehat{\mu}_{a,G}-\sqrt{\widehat{\Sigma}/N}\mathcal{Z}_{\alpha/2},\quad\widehat{\mu}_{a,G}+\sqrt{\widehat{\Sigma}/N}\mathcal{Z}_{\alpha/2}\right],$
(10)

where $\mathcal{Z}_{\alpha/2}$ is the $1-\alpha/2$ quantile for a standard
normal distribution.

##  5 Theoretical Guarantees

###  5.1 Site-specific Estimator

We first establish the theoretical properties of the site-specific estimators
constructed with the multiply robust model-mixing approach. Define
$\overline{\pi}_{a,k}^{j}$, $\overline{m}_{a,k}^{l}$, $\overline{\tau}_{a,k}$
and $\overline{\zeta}_{k}$ as non-stochastic functionals that the
corresponding estimators $\widehat{\pi}_{a,k}^{j}$, $\widehat{m}_{a,k}^{l}$,
$\widehat{\tau}_{a,k}$ and $\widehat{\zeta}_{k}$ converge to. That is,

$\displaystyle\|\widehat{\pi}_{a,k}^{j}-\overline{\pi}_{a,k}^{j}\|$
$\displaystyle=o_{p}(1),\quad\|\widehat{m}_{a,k}^{l}-\overline{m}_{a,k}^{l}\|=o_{p}(1),\quad\|\widehat{\tau}_{a,k}-\overline{\tau}_{a,k}\|=o_{p}(1),\quad\|\widehat{\zeta}_{k}-\overline{\zeta}_{k}\|$
$\displaystyle=o_{p}(1).$

As shown in Lemmas 1 and 2 in Appendix D.2, the $L_{2}$ risks of the model
mixing estimators $\widehat{\pi}_{a,k}$ and $\widehat{m}_{a,k}$ are bounded by
the smallest risks of all candidate models plus a remainder term that vanishes
at a faster rate than the risks themselves.

######  Theorem 2.

Suppose that the conditions in Lemmas 1 and 2 hold, and that
$\widehat{\pi}^{j}_{a,k}$, $\widehat{m}^{l}_{a,k}$, $\widehat{\zeta}_{k}$,
$\widehat{\tau}_{a,k}$, $\bar{\pi}^{j}_{a,k}$, $\bar{m}^{l}_{a,k}$,
$\bar{\zeta}_{k}$ and $\bar{\tau}_{a,k}$ are uniformly bounded for all
treatment models $j\in\mathcal{J}$ and for all outcome models
$l\in\mathcal{L}$. Consider the following conditions:

  1.     1. (B1)

$\overline{\pi}_{a,k}^{j}=\pi_{a,k}$ for some $j\in\mathcal{J}$,

    2. (B2)

$\overline{m}_{a,k}^{l}=m_{a,k}$ for some $l\in\mathcal{L}$,

  2.     1. (C1)

$\overline{\zeta}_{k}=\zeta_{k}$,

    2. (C2)

$\overline{\tau}_{a,k}=\tau_{a,k}$.

Then, under Assumptions (A1) \- (A7), and if one of (B1) or (B2) and one of
(C1) or (C2) hold,

$\displaystyle\|\widehat{\mu}_{a,k}-\mu_{a,T}\|=O_{p}\left(n^{-1/2}+\|\widehat{\pi}_{a,k}-\pi_{a,k}\|\|\widehat{m}_{a,k}-{m}_{a,k}\|+\|\widehat{\zeta}_{k}-{\zeta}_{k}\|\|\widehat{\tau}_{a,k}-\tau_{a,k}\|\right).$
(11)

Further, if the nuisance estimators satisfy the following convergence rate

$\displaystyle\left\|\widehat{m}_{a,k}-m_{a,k}\right\|\left\|\widehat{\pi}_{a,k}-\pi_{a,k}\right\|$
$\displaystyle=o_{p}(1/\sqrt{n}),\quad\|\widehat{\zeta}_{k}-\zeta_{k}\|\left\|\widehat{\tau}_{a,k}-\tau_{a,k}\right\|=o_{p}(1/\sqrt{n}),$

then $\sqrt{n}(\widehat{\mu}_{a,k}-\mu_{a,T})$ asymptotically converges to a
normal distribution with mean zero and asymptotic variance equal to the
semiparametric efficiency bound. The derivation of the result is provided in
the Appendix.

###  5.2 Federated Global Estimator

######  Theorem 3.

Under Assumptions (A1) \- (A7) and the regularity conditions specified in the
Appendix, the federated global estimator of $\Delta_{T}$, given by
$\widehat{\Delta}_{G}=\widehat{\mu}_{1,G}-\widehat{\mu}_{0,G}$, is consistent
and asymptotically normal,

$\sqrt{N/\widehat{\mathcal{V}}}\left(\widehat{\Delta}_{G}-\Delta_{T}\right)\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,1),$ (12)

with the variance estimated consistently as $\widehat{\mathcal{V}}$. The
variance of $\widehat{\Delta}_{G}$ is no larger than that of the estimator
based on target data only,
$\widehat{\Delta}_{T}=\widehat{\mu}_{1,T}-\widehat{\mu}_{0,T}$. Further, if
there exist some source sites with consistent estimators of $\Delta_{T}$ and
satisfy conditions specified in the Appendix, the variance of
$\widehat{\Delta}_{G}$ is strictly smaller than $\widehat{\Delta}_{T}$.

##  6 Experiments

We evaluate the finite sample properties of five different estimators: (i) an
augmented inverse probability weighted (AIPW) estimator using data from the
target site only (Target), (ii) an AIPW estimator that weights each site
proportionally to its sample size (SS), (iii) an AIPW estimator that weights
each site inverse-proportionally to its variance (IVW), (iv) an AIPW estimator
that weights each site with the $L_{1}$ weights defined in (9) (AIPW-$L_{1}$),
and (v) a multiply robust estimator with the $L_{1}$ weights defined in (9)
(MR-$L_{1}$).

Across different settings, we examine the performance of each estimator in
terms of bias, root mean square error, and coverage and length of 95%
confidence intervals (CI) across $500$ simulations.

We consider a total of five sites and fix the first site as the target site
with a relatively small sample size of $300$. The source sites have larger
sample sizes of $\left\\{500,500,1000,1000\right\\}$. We model heterogeneity
in the covariate distributions across sites with skewed normal distributions
and varying levels of skewness in each site,
$X_{kp}\sim\mathcal{SN}\left(x;\Xi_{kp},\Omega_{kp}^{2},\mathrm{~{}A}_{kp}\right)$,
where $k\in\\{1,...,5\\}$ indexes each site and $p\in\\{1,...,4\\}$ indexes
the covariates; $\Xi_{kp}$, $\Omega_{kp}^{2}$ and $\mathrm{~{}A}_{kp}$ are the
location, scale, and skewness parameters, respectively. Following [17], we
also generate covariates $Z_{kp}$ as non-linear transformation of $X_{kp}$
such that $Z_{k1}=\exp(X_{k1}/2)$, $Z_{k2}=X_{k2}/\\{1+\exp(X_{k1})\\}+10$,
$Z_{k3}=(X_{k1}X_{k3}/25+0.6)^{3}$ and $Z_{k4}=(X_{k2}+X_{k4}+20)^{2}$.

For the MR-$L_{1}$, we adaptively mix two outcome models and two treatment
models. We specify the first model with the covariates $X_{kp}$, and the
second model with the covariates $Z_{kp}$. The AIPW-$L_{1}$ estimator requires
a common model to be specified across sites, so we specify the outcome and
treatment models using covariates $X_{kp}$.

The tuning parameter $\lambda$ is selected through cross-validation using a
grid of values $\\{0,10^{-3},10^{-2},0.1,0.5,1,2,5,10\\}$. To perform cross-
validation, the simulated datasets in each site are split into two equally
sized training and validation datasets.

###  6.1 No Covariate Mismatch

We first consider the setting where there is no covariate mismatch, i.e. $p=4$
for both target and source sites. For each unit, we generate potential
outcomes as

$Y_{k}(a)=210+X_{k}\beta_{x}+Z_{k}\beta_{z}+\varepsilon_{k}$ (13)

where $\beta_{x}=\beta_{z}=(27.4,13.7,13.7,13.7)$. For units in the target
site, we generate outcomes with $X_{k}$ only by setting $\beta_{z}=0$; for
units in the source sites, either $X_{k}$ or $Z_{k}$ is used to generate
outcomes. If $\beta_{x}\neq 0$, then $\beta_{z}=0$ and vice versa. Similarly,
the treatment is generated as

$A_{k}\sim\operatorname{Bernoulli}\left(\pi_{k}\right)\quad\pi_{k}=\operatorname{expit}(X_{k}\alpha_{x}+Z_{k}\alpha_{z})$
(14)

where $\alpha_{x}=\alpha_{z}=(-1,0.5,-0.25,-0.1)$. For units in the target
site, we generate treatments with $X_{k}$ only by setting $\alpha_{z}=0$; for
units in the source sites, either $X_{k}$ or $Z_{k}$ is used to generate
treatments. If $\alpha_{x}\neq 0$, then $\alpha_{z}=0$ and vice versa. With
this data generation scheme, the true ATE is $\Delta_{T}=0$.

We compare the performance of the five estimators described above under the
following settings:

Setting 1 ($C=0$): outcomes and treatments in all source sites are generated
with $Z_{k}$. However, all source sites misspecify both models with $X_{k}$.
The target site correctly specifies both models.

Setting 2 ($C=1/2$): outcomes and treatments are generated with $X_{k}$ in
Sites 2 and 4, but with $Z_{k}$ in Sites 3 and 5; thus, the outcome and
treatment models are misspecified in Sites 3 and 5, and only half of the
source sites correctly specify the models.

Setting 3 ($C=1$): outcomes and treatments in all source sites are generated
with $X_{k}$, so all source sites correctly specify outcome and treatment
models with $X_{k}$.

Table 1: Mean absolute error (MAE), root mean squared error (RMSE), coverage (Cov.), and length (Len.) of $95\%$ CIs based on 500 simulated data sets in three (mis)specification settings. | Target | SS | IVW |  AIPW-$L_{1}$ |  MR-$L_{1}$
---|---|---|---|---|---
$C=0$ |  |  |  |  |
MAE | 0.109 | 1.933 | 0.177 | 0.110 | 0.050
RMSE | 0.141 | 1.987 | 0.219 | 0.144 | 0.061
Cov. | 0.950 | 0.998 | 0.826 | 0.936 | 0.960
Len. | 0.551 | 7.035 | 0.567 | 0.547 | 0.234
$C=1/2$ |  |  |  |  |
MAE | 0.109 | 1.111 | 0.107 | 0.109 | 0.050
RMSE | 0.141 | 1.189 | 0.139 | 0.140 | 0.062
Cov. | 0.950 | 1.000 | 0.942 | 0.950 | 0.962
Len. | 0.551 | 6.010 | 0.540 | 0.547 | 0.242
$C=1$ |  |  |  |  |
MAE | 0.109 | 0.036 | 0.035 | 0.050 | 0.049
RMSE | 0.141 | 0.045 | 0.044 | 0.064 | 0.063
Cov. | 0.950 | 0.968 | 0.956 | 0.958 | 0.960
Len. | 0.551 | 0.195 | 0.191 | 0.260 | 0.253

The results in Table 1 indicate that the MR-$L_{1}$ estimator has lower RMSE
than the Target estimator when some source sites have correctly specified
models ($C=1/2$ and $C=1$). Relative to the MR-$L_{1}$ estimator, the SS and
IVW estimators demonstrate larger biases and RMSE, and lower coverage when
some source sites have misspecified models ($C=0$ and $C=1/2)$. The MR-$L_{1}$
estimator shows reduced biases and RMSE compared to the AIPW-$L_{1}$
estimator, while maintaining similar coverage; this improvement can be
attributed to the inclusion of an additional model that closely resembles the
true model. When all source sites correctly specify working models ($C=1$),
the IVW estimator performs optimally with the shortest confidence interval as
expected.

###  6.2 Covariate Mismatch

To demonstrate that our proposed MR-$L_{1}$ estimator can handle covariate
mismatch across sites, we modify the data-generating process in the following
way: only two covariates are used in the outcome and treatment generation
processes in the target site. Specifically, the generating models remain the
same as in (13) and (14), using covariates $X_{k}$. However, for units in the
target site, we set $\beta_{x}=(27.4,13.7,0,0)$ for outcome generation and
$\alpha_{x}=(-1,0.5,0,0)$ for treatment generation.

The AIPW-$L_{1}$ estimator, which requires common models across sites, only
uses the shared covariates ($X_{k1}$ and $X_{k2}$) to specify outcome and
treatment models for all sites. On the other hand, our MR-$L_{1}$ estimator
allows for different covariates in different sites, so we utilize both shared
covariates with the target site and unique covariates to specify the outcome
and treatment models in the source sites.

Table 2: Mean absolute error (MAE), root mean squared error (RMSE), coverage (Cov.), and length (Len.) of $95\%$ CIs based on 500 simulated data sets in covariate mismatch settings. | Target | SS | IVW |  AIPW-$L_{1}$ |  MR-$L_{1}$
---|---|---|---|---|---
MAE | 0.108 | 4.331 | 0.150 | 0.107 | 0.053
RMSE | 0.136 | 4.401 | 0.186 | 0.134 | 0.067
Cov. | 0.946 | 1.000 | 0.882 | 0.950 | 0.944
Len. | 0.538 | 26.024 | 0.553 | 0.536 | 0.253

In Table 2, we observe that the AIPW-$L_{1}$ estimator exhibits similar bias,
RMSE, coverage, and length of confidence intervals as the Target estimator
while outperforming the SS and IVW estimators. This is because relying solely
on shared covariates leads to significant biases in all source sites (Figure
2, left panel), and the AIPW-$L_{1}$ estimator assigns nearly all of the
ensemble weight to the target site so as to reduce bias.

In contrast, the MR-$L_{1}$ estimator outperforms the AIPW-$L_{1}$ estimator
by exhibiting substantially smaller bias, lower RMSE, and better coverage.
This improvement can be attributed to the inclusion of unique covariates from
the source sites, which allows for the recovery of true models in those sites
and contributes to a more accurate estimation of $\Delta_{T}$ (Figure 2, right
panel). These findings suggest that neglecting covariate mismatch by solely
relying on shared covariates can lead to highly biased results.

Figure 2: Estimates of the TATE based on $500$ simulated data sets with
covariate mismatch comparing the site-specific estimators with nuisance
functions estimated by AIPW (left) and by multiply robust model-mixing
(right).

##  7 Conclusion

We have proposed a novel federated approach for privacy-preserving, multiply
robust, and flexible estimation of causal effects. Compared to existing
federated methods, our proposed approach accommodates covariate shift and
covariate mismatch across sites, while guaranteeing efficient estimation and
preserving privacy in the sense that only covariate means of the target
samples are shared in a single round of communication. Our proposal allows
investigators in each site to have greater flexibility in specifying candidate
models by utilizing site-specific information. Moreover, our method utilizes
adaptive ensemble weights to avoid negative transfer in the federation
process. The limitations of the current proposal provide opportunities for
further research. To handle high-dimensional covariates, future research can
explore ways to jointly model the propensity score and density ratio to reduce
the dimension of parameters for population balancing.

## References

  * [1] A. L. Buchanan, M. G. Hudgens, S. R. Cole, K. R. Mollan, P. E. Sax, E. S. Daar, A. A. Adimora, J. J. Eron, and M. J. Mugavero.  Generalizing evidence from randomized trials using inverse probability of sampling weights.  Journal of the Royal Statistical Society. Series A,(Statistics in Society), 181(4):1193, 2018.
  * [2] K. C. G. Chan.  A simple multiply robust estimator for missing response problem.  Stat, 2(1):143–149, 2013.
  * [3] K. C. G. Chan and S. C. P. Yam.  Oracle, multiple robust and multipurpose calibration in a missing response problem.  Statistical Science, 29(3):380–396, 2014.
  * [4] S. Chen and D. Haziza.  Multiply robust imputation procedures for the treatment of item nonresponse in surveys.  Biometrika, 104(2):439–453, 2017.
  * [5] S. R. Cole and E. A. Stuart.  Generalizing evidence from randomized clinical trials to target populations: the actg 320 trial.  American journal of epidemiology, 172(1):107–115, 2010.
  * [6] I. J. Dahabreh, S. E. Robertson, J. A. Steingrimsson, E. A. Stuart, and M. A. Hernan.  Extending inferences from a randomized trial to a new target population.  Statistics in Medicine, 39(14):1999–2014, 2020.
  * [7] I. J. Dahabreh, S. E. Robertson, E. J. Tchetgen, E. A. Stuart, and M. A. Hernán.  Generalizing causal inferences from individuals in randomized trials to all trial-eligible individuals.  Biometrics, 75(2):685–694, 2019.
  * [8] R. Duan, C. Luo, M. H. Schuemie, J. Tong, J. C. Liang, H. H. Chang, M. R. Boland, J. Bian, H. Xu, J. H. Holmes, et al.  Learning from local to global-an efficient distributed algorithm for modeling time-to-event data.  bioRxiv, 2020.
  * [9] R. Duan, Y. Ning, S. Wang, B. Lindsay, R. Carroll, and Y. Chen.  A fast score test for generalized mixture models.  Biometrics, 76:811–820, 2021.
  * [10] Y. Gu and H. Zou.  Aggregated expectile regression by exponential weighting.  Statistica Sinica, 29(2):671–692, 2019.
  * [11] W. Guo, S. Wang, P. Ding, Y. Wang, and M. I. Jordan.  Multi-source causal inference using control variates.  arXiv preprint arXiv:2103.16689, 2021.
  * [12] L. Han, J. Hou, K. Cho, R. Duan, and T. Cai.  Federated adaptive causal estimation (face) of target treatment effects.  arXiv preprint arXiv:2112.09313, 2021.
  * [13] L. Han, Y. Li, B. A. Niknam, and J. R. Zubizarreta.  Privacy-preserving and communication-efficient causal inference for hospital quality measurement.  arXiv preprint arXiv:2203.00768, 2022.
  * [14] P. Han.  A further study of the multiply robust estimator in missing data analysis.  Journal of Statistical Planning and Inference, 148:101–110, 2014\.
  * [15] P. Han.  Multiply robust estimation in regression analysis with missing data.  Journal of the American Statistical Association, 109(507):1159–1173, 2014.
  * [16] P. Han and L. Wang.  Estimation with missing data: beyond double robustness.  Biometrika, 100(2):417–430, 2013.
  * [17] J. D. Y. Kang and J. L. Schafer.  Demystifying double robustness: A comparison of alternative strategies for estimating a population mean from incomplete data.  Statistical Science, 22(4):523–539, 2007.
  * [18] Y. Kawamata, R. Motai, Y. Okada, A. Imakura, and T. Sakurai.  Collaborative causal inference on distributed data.  arXiv preprint arXiv:2208.07898, 2022.
  * [19] E. H. Kennedy, S. Balakrishnan, and L. Wasserman.  Minimax rates for heterogeneous causal effect estimation.  arXiv preprint arXiv:2203.00837, 2022.
  * [20] W. Li, Y. Gu, and L. Liu.  Demystifying a class of multiply robust estimators.  Biometrika, 107(4):919–933, 2020.
  * [21] J. Neyman.  On the application of probability theory to agricultural experiments.  Statistical Science, 5(5):463–480, 1923.
  * [22] J. Qin.  Inferences for case-control and semiparametric two-sample density ratio models.  Biometrika, 85(3):619–630, 1998.
  * [23] D. B. Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of Educational Psychology, 66(5):688, 1974.
  * [24] D. B. Rubin.  Randomization analysis of experimental data: the fisher randomization test comment.  Journal of the American Statistical Association, 75(371):591–593, 1980.
  * [25] E. Tipton.  Improving generalizations from experiments using propensity score subclassification: Assumptions, properties, and contexts.  Journal of Educational and Behavioral Statistics, 38(3):239–266, 2013.
  * [26] T. V. Vo, T. N. Hoang, Y. Lee, and T.-Y. Leong.  Federated estimation of causal effects from observational data.  arXiv preprint arXiv:2106.00456, 2021.
  * [27] R. Xiong, A. Koenecke, M. Powell, Z. Shen, J. T. Vogelstein, and S. Athey.  Federated causal inference in heterogeneous observational data.  arXiv preprint arXiv:2107.11732, 2021.
  * [28] S. Yang and P. Ding.  Combining multiple observational data sources to estimate causal effects.  Journal of the American Statistical Association, 115(531):1540–1554, 2020.
  * [29] Y. Yang.  Adaptive estimation in pattern recognition by combining different procedures.  Statistica Sinica, pages 1069–1089, 2000.
  * [30] Y. Yang.  Combining forecasting procedures: Some theoretical results.  Econometric Theory, 20(1):176–222, 2004.
  * [31] Z. Zeng, E. H. Kennedy, L. M. Bodnar, and A. I. Naimi.  Efficient generalization and transportation.  arXiv preprint arXiv:2302.00092, 2023.

APPENDIX

##  Appendix A Equivalence of density ratio weighting and inverse of selection
probability weighting

We show that the inverse probability of selection weighting (IPSW) to the site
$k$, $\rho_{k}(V_{i})=\frac{1-P(R_{i}=k\mid V_{i}=v)}{P(R_{i}=k\mid V_{i}=v)}$
is equivalent to the density ratio weighting
$\zeta_{k}(V_{i})=\frac{P(V_{i}=v\mid R_{i}=T)}{P(V_{i}=v\mid R_{i}=k)}$.

Under Assumptions (A1) \- (A7), the site-specific estimators based on the IPSW
for $\mu_{a,T}$ is

$\displaystyle\widehat{\mu}_{a,k}$
$\displaystyle=\frac{1}{n_{T}}\sum_{i=1}^{n}\biggl{[}\frac{I(A_{i}=a,R_{i}=k)}{\widehat{\pi}_{a,k}(X_{i})}\widehat{\rho}_{k}(V_{i})\\{Y_{i}-\widehat{m}_{a,k}(X_{i})\\}\biggr{]}$
(15)
$\displaystyle+\frac{1}{n_{T}}\sum_{i=1}^{n}\biggl{[}I(R_{i}=k)\widehat{\rho}_{k}(V_{i})\Bigl{\\{}\widehat{m}_{a,k}(X_{i})-\widehat{\tau}_{a,k}(V_{i})\Bigr{\\}}\biggr{]}$
$\displaystyle+\frac{1}{n_{T}}\sum_{i=1}^{n}I(R_{i}=T)\widehat{\tau}_{a,k}(V_{i})$

where $\widehat{\rho}_{k}(V_{i})=\frac{1-\widehat{P}(R_{i}=k\mid
V_{i}=v)}{\widehat{P}(R_{i}=k\mid V_{i}=v)}$ is an estimator for
$\rho_{k}(V_{i})$. Applying Baye’s rule, we show that the IPSW is equivalent
to the density ratio weighting up to a constant,

$\displaystyle\rho_{k}(V_{i})$ $\displaystyle=\frac{1-P(R_{i}=k\mid
V_{i}=v)}{P(R_{i}=k\mid V_{i}=v)}$ (16) $\displaystyle=\frac{P(R_{i}=T\mid
V_{i}=v)}{P(R_{i}=k\mid V_{i}=v)}$ $\displaystyle=\frac{P(V_{i}=v\mid
R_{i}=T)P(R_{i}=T)}{P(V_{i}=v\mid R_{i}=k)P(R_{i}=k)}$
$\displaystyle={\zeta}_{k}(V_{i})\frac{P(R_{i}=T)}{P(R_{i}=k)}.$

We re-write (15) by substituting $\widehat{\rho}_{k}(V_{i})$ with
$\widehat{\zeta}_{k}(V_{i})$,

$\displaystyle\widehat{\mu}_{a,k}$
$\displaystyle=\frac{1}{n_{T}}\sum_{i=1}^{n}\biggl{[}\frac{I(A_{i}=a,R_{i}=k)}{\widehat{\pi}_{a,k}(X_{i})}{\widehat{\zeta}}_{k}(V_{i})\frac{\widehat{P}(R_{i}=T)}{\widehat{P}(R_{i}=k)}\\{Y_{i}-\widehat{m}_{a,k}(X_{i})\\}\biggr{]}$
(17)
$\displaystyle+\frac{1}{n_{T}}\sum_{i=1}^{n}\biggl{[}I(R_{i}=k){\widehat{\zeta}}_{k}(V_{i})\frac{\widehat{P}(R_{i}=T)}{\widehat{P}(R_{i}=k)}\\{\widehat{m}_{a,k}(X_{i})-\widehat{\tau}_{a,k}(V_{i})\\}\biggr{]}$
$\displaystyle+\frac{1}{n_{T}}\sum_{i=1}^{n}I(R_{i}=T)\widehat{\tau}_{a,k}(V_{i}).$

A reasonable estimator $\frac{\widehat{P}(R_{i}=T)}{\widehat{P}(R_{i}=k)}$ is
$\frac{n_{T}}{n_{k}}$. Therefore, we recover our proposed site-specific
estimator for $\mu_{a,T}$,

$\displaystyle\widehat{\mu}_{a,k}$
$\displaystyle=\frac{1}{n_{k}}\sum_{i=1}^{n}\biggl{[}\frac{I(A_{i}=a,R_{i}=k)}{\widehat{\pi}_{a,k}(X_{i})}{\widehat{\zeta}}_{k}(V_{i})\\{Y_{i}-\widehat{m}_{a,k}(X_{i})\\}\biggr{]}$
(18)
$\displaystyle+\frac{1}{n_{k}}\sum_{i=1}^{n}\biggl{[}I(R_{i}=k){\widehat{\zeta}}_{k}(V_{i})\\{\widehat{m}_{a,k}(X_{i})-\widehat{\tau}_{a,k}(V_{i})\\}\biggr{]}$
$\displaystyle+\frac{1}{n_{T}}\sum_{i=1}^{n}I(R_{i}=T)\widehat{\tau}_{a,k}(V_{i}).$

##  Appendix B Multiply robust estimation for ${m}_{a,k}$

For multiply robust outcome estimation within each site $k$, we consider a set
of $L$ candidate models for conditional outcomes
$\\{m^{l}_{a,k}\left(x\right):l\in\mathcal{L}=\\{1,...,L\\}\\}$. Let
$\hat{m}_{a,k}^{l}(x)$ be the estimates of $m^{l}_{a,k}(x)$ obtained by
fitting the corresponding candidate models, which can be parametric,
semiparametric, or nonparametric machine learning models. Let
$\widehat{m}_{a,k}(X_{i})=\sum_{l=1}^{L}\widehat{\Omega}_{l}\widehat{m}_{a,k}^{l}(X_{i})$
be the predictions with ensemble weights $\widehat{\Omega}_{l}$ of candidate
outcome models under treatment $a$ in site $k$. We derive the ensemble weights
$\widehat{\Omega}_{l}$ based on the cumulative predictive risk of candidate
models. In particular, we denote the data corresponding to treated and control
units as $D_{k,1}$ and $D_{k,0}$ respectively, and the sample sizes of
$D_{k,1}$ and $D_{k,0}$ are denoted as $n_{k,1}$ and $n_{k,0}$. Consider the
treated samples first; we randomly partition $D_{k,1}$ into a training set
$D^{\text{train}}_{k,1}$ of units $i\in\\{1,...,n_{k,1}^{\text{train}}\\}$ and
a validation set $D_{k,1}^{\text{val}}$ of units
$i\in\\{n_{k,1}^{\text{train}}+1,...,n_{k,1}\\}$. Then, each candidate outcome
model is fit on $D^{\text{train}}_{k,1}$ to obtain
$\widehat{m}_{a,n_{k,1}^{\text{train}}}^{l}$ for $l\in\mathcal{L}$.

If the outcome is binary, the ensemble weights are determined by the fitted
models’ predictive risks evaluated on $D^{\text{val}}_{k,1}$ according to the
Bernoulli likelihood as shown in (7). If the outcome is continuous, the
ensemble weights are alternatively determined by the mean squared errors of
the fitted models on $D^{\text{val}}_{k,1}$ [30, 10]. Specifically,

$\displaystyle\widehat{\Omega}_{l}$
$\displaystyle=\left(n_{k,1}-{n_{k,1}^{\text{train}}}\right)^{-1}\sum_{i=n_{k,1}^{\text{train}}+1}^{n_{k,1}}\widehat{\Omega}_{l,i}\text{
and }$ (19) $\displaystyle\widehat{\Omega}_{l,i}$
$\displaystyle=\frac{\exp\left[-\kappa\sum_{q=n_{k,1}^{\text{train}}+1}^{i-1}\left\\{Y_{q}-\widehat{m}_{n_{k,1}^{\text{train}}}^{l}\left(X_{q}\right)\right\\}^{2}\right]}{\sum_{l^{\prime}=1}^{L}\exp\left[-\kappa\sum_{q=n_{k,1}^{\text{train}}+1}^{i-1}\left\\{Y_{q}-\widehat{m}_{n_{k,1}^{\text{train}}}^{l^{\prime}}\left(X_{q}\right)\right\\}^{2}\right]}\quad\text{
for }n_{k,1}^{\text{train}}+2\leq i\leq n_{k,1},$

where $\widehat{\Omega}_{l,n_{k,1}^{\text{train}}+1}=1/L$ and the ensemble
predictions for the conditional outcomes under treatment is
$\widehat{m}_{1,k}(X_{i})=\sum_{l=1}^{L}\widehat{\Omega}_{l}\widehat{m}_{1,k}^{l}(X_{i})$.
The above procedure is then repeated in the control group $D_{k,0}$ and $a=0$
to obtain the ensemble predictions for the conditional outcomes under control.

The tuning parameter $\kappa$ in (19) can be selected via cross-validation and
[10] showed that the performance of model mixture estimators is generally
robust across different choices of $\kappa$; they recommended choosing
$\kappa=\max(1,\lfloor\log(L)\rfloor)$.

##  Appendix C Additional experimental details

We consider five different estimators: (i) an augmented inverse probability
weighted (AIPW) estimator using data from the target site only (Target), (ii)
an AIPW estimator that weights each site proportionally to its sample size
(SS), (iii) an AIPW estimator that weights each site inverse-proportionally to
its variance (IVW), (iv) an AIPW estimator that weights each site with the
$L_{1}$ weights defined in (9) (AIPW-$L_{1}$), and (v) a multiply robust
estimator with the $L_{1}$ weights defined in (9) (MR-$L_{1}$).

###  C.1 Distribution of covariates

We consider a total of five sites and fix the first site as the target site
with a relatively small sample size of $300$. The source sites have larger
sample sizes of $\left\\{500,500,1000,1000\right\\}$. We model heterogeneity
in the covariate distributions across sites with skewed normal distributions
and varying levels of skewness in each site,
$X_{kp}\sim\mathcal{SN}\left(x;\Xi_{kp},\Omega_{kp}^{2},\mathrm{~{}A}_{kp}\right)$,
where $k\in\\{1,...,5\\}$ indexes each site and $p\in\\{1,...,4\\}$ indexes
the covariates. Specifically, for any given site $k\in\mathcal{K}$, we set the
location parameters $\Xi_{kp}=0$ and the scale parameter $\Omega_{kp}=1$.
Moreover, in the case of the target site, we specifically assign the skewness
parameter $\mathrm{~{}A}_{kp}$ a value of zero, denoting a symmetrical
distribution. However, for source site $k\in\mathcal{S}$, we adopt different
values for $\mathrm{~{}A}_{kp}$ based on the sample size. If the sample size
is equal to $500$, we assign $\mathrm{~{}A}_{kp}$ the value of $(1/2)^{p}$,
reflecting a positively skewed distribution. On the other hand, if the sample
size is $1000$, we assign $\mathrm{~{}A}_{kp}$ the value of $-(1/2)^{p}$,
indicating a negatively skewed distribution. Following [17], we generate
covariates $Z_{kp}$ as non-linear transformation of $X_{kp}$ such that
$Z_{k1}=\exp(X_{k1}/2)$, $Z_{k2}=X_{k2}/\\{1+\exp(X_{k1})\\}+10$,
$Z_{k3}=(X_{k1}X_{k3}/25+0.6)^{3}$ and $Z_{k4}=(X_{k2}+X_{k4}+20)^{2}$.

###  C.2 Site-specific estimates

The main text presented simulation results in Table 1 based on $500$
simulations, assuming no covariate mismatch. Here, to provide further insights
into the site-specific estimates, we provide simulation results in Figure 3.
The results show that in cases where some sites fail to properly specify their
working models, AIPW estimators exhibit significant bias, while the multiply
robust estimators are capable of accurately recovering the true TATE. This can
be attributed to the fact that the additional candidate model closely
approximates the true underlying models. These findings underscore the
enhanced safeguard against model misspecification provided by the multiply
robust estimators.

(a) No Site Correctly Specify Models ($C=0$)

(b) Site 2 and Site 4 Correctly Specify Models ($C=1/2$)

(c) All Sites Correctly Specify Models ($C=1$)

Figure 3: Results of site-specific estimation based on $500$ simulated data
sets in three (mis)specification settings.

###  C.3 Reproducibility and computing resources

All experiments in this study were performed using the statistical programming
language R (version 4.2.2). The package sn (version 2.1.0) was employed for
generating covariates that follow a skewed-normal distribution within each
site. To solve the density ratio estimating equations, we utilized the
rootSolve package (version 1.8.2.3). For the estimation of adaptive ensemble
weights based on penalized regression of site-specific influence functions, we
employed the glmnet package (version 4.1-4).

To enhance computational efficiency, parallel computing packages were
employed. Specifically, the packages foreach (version 1.5.2) and doParallel
were employed to facilitate the replication of experiments. For the purpose of
model mixing in multiply robust estimation, we employed the parallel package
(version 4.2.2). The replication of experiments was carried out using ten CPU
cores, while the implementation of the model mixing algorithm utilized five
CPU cores.

##  Appendix D Proofs

###  D.1 Proof of Theorem 1

If Assumptions (A1) \- (A3) hold, the mean counterfactual outcomes in the
target population can be identified using data from the target site;

$\displaystyle\mu_{a,T}$ $\displaystyle=E\left\\{Y_{i}{(a)}\mid
R_{i}=T\right\\}$ (20) $\displaystyle=E[E\left\\{Y_{i}{(a)}\mid
V_{i}=v,R_{i}=T\right\\}\mid R_{i}=T]$
$\displaystyle=E[E\left\\{Y_{i}{(a)}\mid V_{i}=v,A_{i}=a,R_{i}=T\right\\}\mid
R_{i}=T]$ $\displaystyle=E[E\left\\{Y_{i}\mid
V_{i}=v,A_{i}=a,R_{i}=T\right\\}\mid R_{i}=T].$

The second line follows the law of total expectation; the third line follows
Assumption (A2); the last line follows Assumption (A1).

If Assumptions (A1), (A4) \- (A7) hold, the mean counterfactual outcomes in
the target population can be identified using data from source sites;

$\displaystyle\mu_{a,T}$ $\displaystyle=E\left\\{Y_{i}{(a)}\mid
R_{i}=T\right\\}$ (21) $\displaystyle=E[E\left\\{Y_{i}{(a)}\mid
V_{i}=v,R_{i}=T\right\\}\mid R_{i}=T]$
$\displaystyle=E[E\left\\{Y_{i}{(a)}\mid V_{i}=v,R_{i}=k\right\\}\mid
R_{i}=T]$
$\displaystyle=E\left(E[E\left\\{Y_{i}{(a)}\mid{X}_{i}=x,V_{i}=v,R_{i}=k\right\\}\mid
V_{i}=v,R_{i}=k]\mid R_{i}=T\right)$
$\displaystyle=E\left(E[E\left\\{Y_{i}{(a)}\mid{X}_{i}=x,R_{i}=k\right\\}\mid
V_{i}=v,R_{i}=k]\mid R_{i}=T\right)$
$\displaystyle=E\left(E[E\left\\{Y_{i}(a)\mid{X}_{i}=x,R_{i}=k,A_{i}=a\right\\}\mid
V_{i}=v,R_{i}=k]\mid R_{i}=T\right)$
$\displaystyle=E\left(E[E\left\\{Y_{i}\mid{X}_{i}=x,R_{i}=k,A_{i}=a\right\\}\mid
V_{i}=v,R_{i}=k]\mid R_{i}=T\right).$

The second line follows the law of total expectation; the third line follows
Assumption (A6); the fourth line follows the law of total expectation; the
fifth line follows by our setup that $V_{i}\subseteq X_{i}$; the sixth line
follows Assumption (A4) the last line follows Assumption (A1).

###  D.2 Lemmas in [20]

We restate Lemmas 1 and 2 that were proved in [20]. These two lemmas together
show that the $L_{2}$ risks of the multiply robust estimators for $\pi_{a,k}$
and $m_{a,k}$ are bounded by the $L_{2}$ risks of the model with the smallest
risks, plus a negligible remainder term.

######  Lemma 1.

Suppose that for each $j\in\mathcal{J}$, there exists a constant
$0<\epsilon_{j}<1/2$ such that
$\epsilon_{j}<\widehat{\pi}_{a,k}^{j}(x)<1-\epsilon_{j}$ for all $x$. Then

$E\left(\left\|\widehat{\pi}_{a,k}-\pi_{a,k}\right\|^{2}\right)\leq\inf_{j\in\mathcal{J}}\frac{2}{\epsilon_{j}^{2}}E\left(\|\widehat{\pi}_{a,k}^{j}-\pi_{a,k}\|^{2}\right)+\frac{2\log(J)}{n_{k}-n_{k}^{\text{train}}}$
(22)

where
$\widehat{\pi}_{a,k}(x)=\sum_{j=1}^{J}\widehat{\Lambda}_{j}\widehat{\pi}^{j}_{a,k}(x)$.

If the candidate propensity models are parametric and one of them is correctly
specified for $\pi_{a,k}$, then
$\inf_{j\in\mathcal{J}}\frac{2}{\epsilon_{j}^{2}}E(\|\widehat{\pi}_{a,k}^{j}-\pi_{a,k}\|^{2})$
converges at a rate of $1/n$. If the candidate propensity models are non-
parametric, then
$\inf_{j\in\mathcal{J}}\frac{2}{\epsilon_{j}^{2}}E(\|\widehat{\pi}_{a,k}^{j}-\pi_{a,k}\|^{2})$
converges at a rate slower than $1/n$. The remainder term
$\frac{2\log(J)}{n_{k}-n_{k}^{\text{train}}}$ converges at the rate
$1/n_{k}-n_{k}^{\text{train}}$, so it vanishes at a faster rate than the
statistical risks of candidate models themselves.

######  Lemma 2.

Suppose we have a continuous outcome and follow the continuous outcome model
mixing algorithm presented in Appendix B. Suppose there exist constants
$C_{1},C_{2}>0$ such that
$\sup_{l\in\mathcal{L}}|\hat{m}_{a,k}^{l}(x)-m_{a,k}(x)|\leq C_{1}$ for all
$x$ and the subexponential norm of $Y-m_{a,k}(x)$ given $X=x$ is bounded above
by $C_{2}$ for all $x$. Then for $a\in\\{0,1\\}$

$E\left(\left\|\widehat{m}_{a,k}-m_{a,k}\right\|^{2}\right)\leq\inf_{l\in\mathcal{L}}E\left(\|\widehat{m}_{a,k}^{l}-m_{a,k}\|^{2}\right)+\frac{\log(L)}{\kappa\left(n_{k,a}-n_{k,a}^{\text{train}}\right)}$

for

$0<\kappa\leq\max\left[\frac{1}{16{e}C_{1}C_{2}},\frac{\exp\left\\{C_{1}\left(8{e}C_{2}\right)^{-1}\right\\}}{4\mathcal{M}_{2}\left\\{\left(4{e}C_{2}\right)^{-1}\right\\}+16C_{1}^{2}\mathcal{M}_{0}\left\\{\left(4{e}C_{2}\right)^{-1}\right\\}}\right]$
(23)

where
$\mathcal{M}_{0}(t)=2\exp\left(2{e}^{2}C_{2}^{2}t^{2}\right),\mathcal{M}_{2}(t)=16\sqrt{}2C_{2}^{2}\exp\left(8{e}^{4}C_{2}^{2}t^{2}\right)\text{
and }{e}=\exp(1)$.

###  D.3 Influence function of site-specific estimator

In this section, we derive the general form of the influence functions for the
site-specific estimators, following a similar derivation in [31]. The primary
difference is that we use the density ratio weights $\zeta_{k}(V)$ instead of
the inverse probability of selection weights. Define
$Z_{i}=(Y_{i},A_{i},X_{i},R_{i})$ with the (partial) baseline covariates in
the target site as $V_{i}\subseteq X_{i}$. The general form for our efficient
influence function is

$\displaystyle\xi_{a,k}(Z_{i})$
$\displaystyle=\frac{1}{P(R_{i}=k)}\Biggl{[}\frac{I(A_{i}=a,R_{i}=k)}{\pi_{a,k}(X_{i})}\zeta_{k}(V_{i})\\{Y_{i}-m_{a,k}(X_{i})\\}\Biggr{]}$
(24)
$\displaystyle+\frac{1}{P(R_{i}=k)}\Biggl{[}I(R_{i}=k)\zeta_{k}(V_{i})\\{m_{a,k}(X_{i})-\tau_{a,k}(V_{i})\\}\Biggr{]}$
$\displaystyle+\frac{1}{P(R_{i}=T)}\Biggl{[}I(R_{i}=T)\tau_{a,k}(V_{i})\Biggr{]}-\mu_{a,T}.$

The target site estimator has the following form,

$\widehat{\mu}_{a,T}=\frac{1}{n_{T}}\sum_{i=1}^{n}\Biggl{[}\frac{I(A_{i}=a,R_{i}=T)}{\pi_{a,T}(V_{i})}\left\\{Y_{i}-m_{a,T}(V_{i})\right\\}+I(R_{i}=T)m_{a,T}(V_{i})\Biggr{]}.$
(25)

This is a standard AIPW estimator whose influence function is derived in [7,
13] as

$\xi_{a,T}(Z_{i})=\frac{1}{P(R_{i}=T)}\Biggl{[}\frac{I(A_{i}=a,R_{i}=T)}{\pi_{a,T}(V_{i})}\\{Y_{i}-m_{a,T}(V_{i})\\}+I(R_{i}=T)m_{a,T}(V_{i})\Biggr{]}-\mu_{a,T}.$
(26)

###  D.4 Proof of Theorem 2

We first prove that given the conditions in Theorem 2, the site-specific
estimator is a consistent estimator of $\mu_{a,T}$, formalized in the
following Lemma:

######  Lemma 3.

Given Assumptions (A1) \- (A7), and 1 and 2, if one of (B1) or (B2) and one of
(C1) or (C2) hold, then $\widehat{\mu}_{a,k}$ is a consistent estimator for
$\mu_{a,T}$.

###### Proof of Lemma 3.

We divide the proof into four cases and show that $\widehat{\mu}_{a,k}$
achieves consistency.

  1. Case 1:

When $\overline{\pi}_{a,k}^{j}=\pi_{a,k}$ and
$\overline{\tau}_{a,k}=\tau_{a,k}$
By assumption, we have $\|\widehat{\pi}_{a,k}^{j}-{\pi}_{a,k}\|=o_{p}(1)$ and
$\|\widehat{\tau}_{a,k}-{\tau}_{a,k}\|=o_{p}(1)$. Given Lemma 1 for
$\widehat{\pi}_{a,k}$, $\|\widehat{\pi}_{a,k}-{\pi}_{a,k}\|=o_{p}(1)$. Define
$\overline{m}_{a,k}(x)=\frac{1}{L}\sum_{l=1}^{L}\overline{m}_{a,k}^{l}(x)$. By
definition, $\|\widehat{m}_{a,k}^{l}-\overline{m}_{a,k}^{l}\|=o_{p}(1)$, so
$\|\widehat{m}_{a,k}-\overline{m}_{a,k}\|=o_{p}(1)$. Together with
$\|\widehat{\zeta}_{k}-\overline{\zeta}_{k}\|=o_{p}(1)$, we can re-write the
site-specific estimator as

$\displaystyle\widehat{\mu}_{a,k}=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\biggl{[}\frac{n}{n_{k}}\frac{I(A_{i}=a,R_{i}=k)}{\pi_{a,k}(X_{i})}{\overline{\zeta}}_{k}(V_{i})\Bigl{\\{}Y_{i}-\overline{m}_{a,k}(X_{i})\Bigr{\\}}\biggr{]}$
(27)
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}\biggl{[}\frac{n}{n_{k}}I(R_{i}=k){\overline{\zeta}}_{k}(V_{i})\Bigl{\\{}\overline{m}_{a,k}(X_{i})-\tau_{a,k}(V_{i})\Bigr{\\}}\biggr{]}$
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}\frac{n}{n_{T}}I(R_{i}=T)\tau_{a,k}(V_{i})$
$\displaystyle+o_{p}(1)$

By assumption of i.i.d. units within each site and the law of large numbers,
$\widehat{\mu}_{a,k}$ converges in probability to

$\displaystyle
E\biggl{[}\frac{n}{n_{k}}\frac{I(A_{i}=a,R_{i}=k)}{\pi_{a,k}(X_{i})}{\overline{\zeta}}_{k}(V_{i})\Bigl{\\{}Y_{i}-\overline{m}_{a,k}(X_{i})\Bigr{\\}}\biggr{]}$
(28) $\displaystyle+$ $\displaystyle
E\biggl{[}\frac{n}{n_{k}}I(R_{i}=k){\overline{\zeta}}_{k}(V_{i})\Bigl{\\{}\overline{m}_{a,k}(X_{i})-\tau_{a,k}(V_{i})\Bigr{\\}}\biggr{]}$
$\displaystyle+$ $\displaystyle
E\biggl{[}\frac{n}{n_{T}}I(R_{i}=T)\tau_{a,k}(V_{i})\biggr{]}$
$\displaystyle=$
$\displaystyle\underbrace{E\biggl{[}\frac{n}{n_{k}}I(R_{i}=k){\overline{\zeta}}_{k}(V_{i})\Bigl{\\{}\frac{I(A_{i}=a)}{\pi_{a,k}(X_{i})}Y_{i}-\tau_{a,k}(V_{i})\Bigr{\\}}\biggr{]}}_{T_{1}}$
$\displaystyle-$
$\displaystyle\underbrace{E\biggl{[}\frac{n}{n_{k}}I(R_{i}=k){\overline{\zeta}}_{k}(V_{i})\Bigl{\\{}\frac{I(A_{i}=a)}{\pi_{a,k}(X_{i})}-1\Bigr{\\}}\overline{m}_{a,k}(X_{i})\biggr{]}}_{T_{2}}$
$\displaystyle+$
$\displaystyle\underbrace{E\biggl{[}\frac{n}{n_{T}}I(R_{i}=T)\tau_{a,k}(V_{i})\biggr{]}}_{T_{3}}.$

Given that $\pi_{a,k}$ is the true propensity score model, $E(T_{2})=0$ since

$E\Bigl{\\{}\frac{I(A_{i}=a)}{\pi_{a,k}(X_{i})}-1\mid
X_{i},R_{i}=k\Bigr{\\}}=\frac{P(A_{i}=a\mid
X_{i},R_{i}=k)}{\pi_{a,k}(X_{i})}-1=0.$ (29)

Similarly, given that $\tau_{a,k}=E\\{m_{a,k}(X_{i})\mid V_{i},R_{i}=k\\}$ is
the true model,

$\displaystyle E(T_{3})$ $\displaystyle=E\\{E(T_{3}\mid R_{i}=T)\\}$ (30)
$\displaystyle=\frac{n}{n_{T}}P(R_{i}=T)E(\tau_{a,k}(V_{i})\mid R_{i}=T)$
$\displaystyle=E(\tau_{a,k}(V_{i})\mid R_{i}=T)$ $\displaystyle=\mu_{a,T}.$

Finally, we consider $T_{1}$;

$\displaystyle E(T_{1})$
$\displaystyle=E\biggl{[}\frac{n}{n_{k}}I(R_{i}=k){\overline{\zeta}}_{k}(V_{i})\Bigl{\\{}\frac{I(A_{i}=a)}{\pi_{a,k}(X_{i})}Y_{i}-\tau_{a,k}(V_{i})\Bigr{\\}}\biggr{]}$
$\displaystyle=E\Biggl{(}E\Bigl{[}\frac{n}{n_{k}}I(R_{i}=k)\overline{\zeta}_{k}(V_{i})\Bigl{\\{}\frac{I(A_{i}=a)}{\pi_{a,k}(X_{i})}Y_{i}-\tau_{a,k}(V_{i})\Bigr{\\}}\mid
V_{i}\Bigr{]}\Biggr{)}$
$\displaystyle=E\Biggl{(}\frac{n}{n_{k}}I(R_{i}=k){\overline{\zeta}}_{k}(V_{i})\Bigl{[}E\Bigl{\\{}\frac{I(A_{i}=a)}{\pi_{a,k}(X_{i})}Y_{i}\mid
V_{i}\Bigr{\\}}-\tau_{a,k}(V_{i})\Bigr{]}\Biggr{)}$
$\displaystyle=E\Biggl{(}\frac{n}{n_{k}}I(R_{i}=k){\overline{\zeta}}_{k}(V_{i})\Bigl{[}E\Bigl{\\{}\frac{I(A_{i}=a)}{\pi_{a,k}(X_{i})}Y_{i}\mid
V_{i},R_{i}=k\Bigr{\\}}-\tau_{a,k}(V_{i})\Bigr{]}\Biggr{)}$
$\displaystyle=E\Biggl{\\{}\frac{n}{n_{k}}I(R_{i}=k){\overline{\zeta}}_{k}(V_{i})\Biggl{(}E\Bigl{[}E\Bigl{\\{}\frac{I(A_{i}=a)}{\pi_{a,k}(X_{i})}Y_{i}\mid
X_{i},R_{i}=k\Bigr{\\}}\mid
V_{i},R_{i}=k\Bigl{]}-\tau_{a,k}(V_{i})\Biggr{)}\Biggr{\\}}$
$\displaystyle=E\Biggl{(}\frac{n}{n_{k}}I(R_{i}=k){\overline{\zeta}}_{k}(V_{i})\Bigl{[}E\Bigl{\\{}E(Y_{i}\mid
A_{i}=a,X_{i},R_{i}=k)\mid
V_{i},R_{i}=k\Bigl{\\}}-\tau_{a,k}(V_{i})\Bigr{]}\Biggr{)}$
$\displaystyle=E\Biggl{(}\frac{n}{n_{k}}I(R_{i}=k){\overline{\zeta}}_{k}(V_{i})\Bigl{[}E\Bigl{\\{}m_{a,k}(X_{i})\mid
V_{i},R_{i}=k\Bigl{\\}}-\tau_{a,k}(V_{i})\Bigr{]}\Biggr{)}$
$\displaystyle=0,\text{ which completes the proof for Case 1.}$ (31)

  2. Case 2:

When $\overline{\pi}_{a,k}^{j}=\pi_{a,k}$ and $\overline{\zeta}_{k}=\zeta_{k}$
By assumption, we have $\|\widehat{\pi}_{a,k}^{j}-{\pi}_{a,k}\|=o_{p}(1)$ and
$\|\widehat{\zeta}_{k}-{\zeta}_{k}\|=o_{p}(1)$. Given Lemma 1 for
$\widehat{\pi}_{a,k}$, $\|\widehat{\pi}_{a,k}-{\pi}_{a,k}\|=o_{p}(1)$. Define
$\overline{m}_{a,k}(x)=\frac{1}{L}\sum_{l=1}^{L}\overline{m}_{a,k}^{l}(x)$. By
definition, $\|\widehat{m}_{a,k}^{l}-\overline{m}_{a,k}^{l}\|=o_{p}(1)$, so
$\|\widehat{m}_{a,k}-\overline{m}_{a,k}\|=o_{p}(1)$. Together with
$\|\widehat{\tau}_{a,k}-\overline{\tau}_{a,k}\|=o_{p}(1)$, we can write the
site-specific estimator as

$\displaystyle\widehat{\mu}_{a,k}=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\biggl{[}\frac{n}{n_{k}}\frac{I(A_{i}=a,R_{i}=k)}{\pi_{a,k}(X_{i})}\zeta_{k}(V_{i})\Bigl{\\{}Y_{i}-\overline{m}_{a,k}(X_{i})\Bigr{\\}}\biggr{]}$
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}\biggl{[}\frac{n}{n_{k}}I(R_{i}=k)\zeta_{k}(V_{i})\Bigl{\\{}\overline{m}_{a,k}(X_{i})-\overline{\tau}_{a,k}(V_{i})\Bigr{\\}}\biggr{]}$
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}\frac{n}{n_{T}}I(R_{i}=T)\overline{\tau}_{a,k}(V_{i})$
$\displaystyle+o_{p}(1)$ (32)

By assumption of i.i.d. units within each site and the law of large numbers,
$\widehat{\mu}_{a,k}$ converges in probability to

$\displaystyle
E\biggl{[}\frac{n}{n_{k}}\frac{I(A_{i}=a,R_{i}=k)}{\pi_{a,k}(X_{i})}\zeta_{k}(V_{i})\Bigl{\\{}Y_{i}-\overline{m}_{a,k}(X_{i})\Bigr{\\}}\biggr{]}$
$\displaystyle+$ $\displaystyle
E\biggl{[}\frac{n}{n_{k}}I(R_{i}=k)\zeta_{k}(V_{i})\Bigl{\\{}\overline{m}_{a,k}(X_{i})-\overline{\tau}_{a,k}(V_{i})\Bigr{\\}}\biggr{]}$
$\displaystyle+$ $\displaystyle
E\biggl{[}\frac{n}{n_{T}}I(R_{i}=T)\overline{\tau}_{a,k}(V_{i})\biggr{]}$
$\displaystyle=$
$\displaystyle\underbrace{E\Bigl{\\{}\frac{n}{n_{k}}\frac{I(A_{i}=a,R_{i}=k)}{\pi_{a,k}(X_{i})}\zeta_{k}(V_{i})Y_{i}\Bigr{\\}}}_{T_{1}}$
$\displaystyle+$
$\displaystyle\underbrace{E\biggl{[}\frac{n}{n_{k}}I(R_{i}=k)\zeta_{k}(V_{i})\overline{m}_{a,k}(X_{i})\Bigl{\\{}1-\frac{I(A_{i}=a)}{\pi_{a,k}(X_{i})}\Bigr{\\}}\biggr{]}}_{T_{2}}$
$\displaystyle+$
$\displaystyle\underbrace{E\biggl{[}\overline{\tau}_{a,k}(V_{i})\Bigl{\\{}\frac{n}{n_{T}}I(R_{i}=T)-\frac{n}{n_{k}}I(R_{i}=k)\zeta_{k}({V_{i}})\Bigr{\\}}\biggr{]}}_{T_{3}}$
(33)

Given that $\pi_{a,k}$ is the true propensity score model, $E(T_{2})=0$ since

$E\Bigl{\\{}\frac{I(A_{i}=a)}{\pi_{a,k}(X_{i})}-1\mid
X_{i},R_{i}=k\Bigr{\\}}=\frac{P(A_{i}=a\mid
X_{i},R_{i}=k)}{\pi_{a,k}(X_{i})}-1=0.$ (34)

Also since $\zeta_{k}={f(V_{i}|R_{i}=T)}/{f(V_{i}|R_{i}=k)}$ is the true
density ratio model,

$\displaystyle E(T_{3})$
$\displaystyle=E\biggl{(}E\biggl{[}\overline{\tau}_{a,k}(V_{i})\Bigl{\\{}\frac{n}{n_{T}}I(R_{i}=T)-\frac{n}{n_{k}}I(R_{i}=k)\zeta_{k}({V_{i}})\Bigr{\\}}\mid
V_{i}\biggr{]}\biggr{)}$ (35)
$\displaystyle=E\biggl{(}\overline{\tau}_{a,k}(V_{i})\biggl{[}\frac{n}{n_{T}}E\Bigl{\\{}I(R_{i}=T)\Bigr{\\}}-\frac{n}{n_{k}}E\Bigl{\\{}I(R_{i}=k)\Bigr{\\}}\zeta_{k}({V_{i}})\mid
V_{i}\biggr{]}\biggr{)}$ $\displaystyle=0\text{ following the results from
Appendix \ref{appendix_equivalence_density_ratio}. }$

Thus, we are left with $T_{1}$, which equals to $\mu_{a,T}$ by iterative
conditioning and Assumptions (A1), (A2), (A4) and (A6).

Proofs for Case 3, when $\overline{m}_{a,k}^{l}=m_{a,k}$ and
$\overline{\tau}_{a,k}=\tau_{a,k}$ and Case 4, when
$\overline{m}_{a,k}^{l}=m_{a,k}$ and $\overline{\zeta}_{k}=\zeta_{k}$, follow
closely to the proofs of Cases 1 and 2 and are omitted here for brevity. ∎

###### Proof of Theorem 2.

Given the general form of influence function (24), we can decompose the
estimation risk of the site-specific estimator as

$\mathbb{P}_{n}(\widehat{\mu}_{a,k})-{E}(\mu_{a,T})=\left(\mathbb{P}_{n}-{E}\right)(\widehat{\mu}_{a,k}-\mu_{a,T})+\left(\mathbb{P}_{n}-{E}\right)\mu_{a,T}+{E}(\widehat{\mu}_{a,k}-\mu_{a,T})$
(36)

Assuming proper sample splitting methods are employed, as stated in Lemma 2 of
[19], we have:

$\left(\mathbb{P}_{n}-E\right)\left(\widehat{\mu}_{a,k}-\mu_{a,T}\right)=o_{p}(n^{-1/2}).$
(37)

Then, we consider $\left(\mathbb{P}_{n}-E\right)\mu_{a,T}$. By the Central
Limit Theorem, we have
$\sqrt{n}\left\\{\mathbb{P}_{n}(\mu_{a,T})-E(\mu_{a,T})\right\\}\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,\sigma^{2})$, where $\sigma^{2}$ is given by:

$\displaystyle\sigma^{2}$
$\displaystyle=E\left\\{\xi_{a,k}(Z_{i})^{2}\right\\}$
$\displaystyle=E\left\\{\frac{1}{P^{2}(R_{i}=k)}\Biggl{[}\frac{I(A_{i}=a,R_{i}=k)}{\pi_{a,k}(X_{i})}\zeta_{k}(V_{i})\\{Y_{i}-m_{a,k}(X_{i})\\}\Biggr{]}^{2}\right\\}$
$\displaystyle+E\left\\{\frac{1}{P^{2}(R_{i}=k)}\Biggl{[}I(R_{i}=k)\zeta_{k}(V_{i})\\{m_{a,k}(X_{i})-\tau_{a,k}(V_{i})\\}\Biggr{]}^{2}\right\\}$
$\displaystyle+E\left\\{\frac{1}{P^{2}(R_{i}=T)}\Biggl{[}I(R_{i}=T)\left\\{\tau_{a,k}(V_{i})-\mu_{a,T}\right\\}\Biggr{]}^{2}\right\\}$
$\displaystyle+\text{Remaining cross-terms}.$ (38)

It can be verified that all remaining cross-terms have an expected value of
zero. Hence,

$\displaystyle\sigma^{2}$
$\displaystyle=E\left\\{\frac{1}{P^{2}(R_{i}=k)}\left[\frac{I(A_{i}=a,R_{i}=k)}{\pi^{2}_{a,k}(X_{i})}\zeta^{2}_{k}(V_{i})\\{Y_{i}-m_{a,k}(X_{i})\\}^{2}\right]\right\\}$
$\displaystyle+E\left\\{\frac{1}{P^{2}(R_{i}=k)}I(R_{i}=k)\left[\zeta_{k}(V_{i})\\{m_{a,k}(X_{i})-\tau_{a,k}(V_{i})\\}\right]^{2}\right\\}$
$\displaystyle+E\left[\frac{1}{P^{2}(R_{i}=T)}I(R_{i}=T)\left\\{\tau_{a,k}(V_{i})-\mu_{a,T}\right\\}^{2}\right]$
$\displaystyle=E\left[E\left\\{\frac{1}{P^{2}(R_{i}=k)}\left[\frac{I(A_{i}=a,R_{i}=k)}{\pi^{2}_{a,k}(X_{i})}\zeta^{2}_{k}(V_{i})\\{Y_{i}-m_{a,k}(X_{i})\\}^{2}\right]\mid
X_{i}\right\\}\right]$
$\displaystyle+E\left[E\left\\{\frac{1}{P^{2}(R_{i}=k)}I(R_{i}=k)\zeta^{2}_{k}(V_{i})\\{m_{a,k}(X_{i})-\tau_{a,k}(V_{i})\\}^{2}\mid
V_{i}\right\\}\right]$
$\displaystyle+E\left[\frac{1}{P^{2}(R_{i}=T)}I(R_{i}=T)\left\\{\tau_{a,k}(V_{i})-\mu_{a,T}\right\\}^{2}\right]$
$\displaystyle=E\left[\frac{P(R_{i}=k\mid
X_{i})}{P^{2}(R_{i}=k)}\zeta^{2}_{k}(V_{i})E\left\\{\frac{I(A_{i}=a)}{\pi^{2}_{a,k}(X_{i})}\\{Y_{i}-m_{a,k}(X_{i})\\}^{2}\mid
X_{i},A_{i}=a,R_{i}=k\right\\}\right]$
$\displaystyle+E\left\\{\frac{P(R_{i}=k\mid
X_{i})}{P^{2}(R_{i}=k)}\zeta^{2}_{k}(V_{i})E\left[\\{m_{a,k}(X_{i})-\tau_{a,k}(V_{i})\\}^{2}\mid
X_{i},R_{i}=k\right]\right\\}$
$\displaystyle+E\left(\frac{P(R_{i}=T)}{P^{2}(R_{i}=T)}E\left[\left\\{\tau_{a,k}(V_{i})-\mu_{a,T}\right\\}^{2}\mid
R_{i}=T\right]\right)$ $\displaystyle=E\left[\frac{P(R_{i}=k\mid
V_{i})}{P^{2}(R_{i}=k)}\zeta^{2}_{k}(V_{i})\frac{\text{Var}\\{Y_{i}\mid
X_{i},A_{i}=a,R_{i}=k\\}}{\pi_{a,k}(X_{i})}\right]$
$\displaystyle+E\left[\frac{P(R_{i}=k\mid
V_{i})}{P^{2}(R_{i}=k)}\zeta^{2}_{k}(V_{i})\text{Var}\\{m_{a,k}(X_{i})\mid
V_{i},R_{i}=k\\}\right]$
$\displaystyle+E\left[\frac{\text{Var}\left\\{\tau_{a,k}(V_{i})\mid
R_{i}=T\right\\}}{P(R_{i}=T)}\right].$ (39)

This expression represents the semiparametric efficiency bound for estimating
$\mu_{a,T}$ and is finite given that all nuisance functions are uniformly
bounded. Consequently, we can conclude that
$(\mathbb{P}_{n}-E)(\mu_{a,T})=O_{p}(n^{-1/2})$.

Moving forward, we analyze the conditional bias term. Firstly, let us define
$E\left\\{\widehat{m}_{a,k}(X_{i})\mid
V_{i},R_{i}=k\right\\}=\tilde{\tau}_{a,k}(V_{i})$. Given conditions in Theorem
2 and Lemma 3, we can rewrite the conditional bias as follows:

$\displaystyle E\left(\widehat{\mu}_{a,k}-\mu_{a,T}\right)$
$\displaystyle=E\left(\frac{n}{n_{k}}\biggl{[}\frac{I(A_{i}=a,R_{i}=k)}{\widehat{\pi}_{a,k}(X_{i})}\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}m_{a,k}(X_{i})-\widehat{m}_{a,k}(X_{i})\Bigr{\\}}\biggr{]}\right)$
$\displaystyle+E\left(\frac{n}{n_{k}}\biggl{[}I(R_{i}=k)\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}\tilde{\tau}_{a,k}(V_{i})-\widehat{\tau}_{a,k}(V_{i})\Bigr{\\}}\biggr{]}\right)$
$\displaystyle+E\left(\frac{n}{n_{T}}\biggl{[}I(R_{i}=T)\Bigl{\\{}\widehat{\tau}_{a,k}(V_{i})-\mu_{a,T}(V_{i})\Bigr{\\}}\biggr{]}\right)$
$\displaystyle=E\left(\frac{n}{n_{k}}\biggl{[}\frac{I(A_{i}=a,R_{i}=k)}{\widehat{\pi}_{a,k}(X_{i})}\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}m_{a,k}(X_{i})-\widehat{m}_{a,k}(X_{i})\Bigr{\\}}\biggr{]}\right)$
$\displaystyle+E\left(\frac{n}{n_{k}}\biggl{[}I(R_{i}=k)\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}\tilde{\tau}_{a,k}(V_{i})-\tau_{a,k}(V_{i})\Bigr{\\}}\biggr{]}\right)$
$\displaystyle-E\left(\frac{n}{n_{k}}\biggl{[}I(R_{i}=k)\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}\widehat{\tau}_{a,k}(V_{i})-{\tau}_{a,k}(V_{i})\Bigr{\\}}\biggr{]}\right)$
$\displaystyle+E\left(\frac{n}{n_{T}}\biggl{[}I(R_{i}=T)\Bigl{\\{}\widehat{\tau}_{a,k}(V_{i})-\mu_{a,T}(V_{i})\Bigr{\\}}\biggr{]}\right)$
$\displaystyle=E\left(\frac{n}{n_{k}}\biggl{[}\frac{I(A_{i}=a,R_{i}=k)}{\widehat{\pi}_{a,k}(X_{i})}\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}m_{a,k}(X_{i})-\widehat{m}_{a,k}(X_{i})\Bigr{\\}}\biggr{]}\right)$
$\displaystyle+E\left(\frac{n}{n_{k}}\biggl{[}I(R_{i}=k)\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}\widehat{m}_{a,k}(X_{i})-m_{a,k}(X_{i})\Bigr{\\}}\biggr{]}\right)$
$\displaystyle-E\left(\frac{n}{n_{k}}\biggl{[}I(R_{i}=k)\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}\widehat{\tau}_{a,k}(V_{i})-{\tau}_{a,k}(V_{i})\Bigr{\\}}\biggr{]}\right)$
$\displaystyle+E\left(\frac{n}{n_{T}}\biggl{[}I(R_{i}=T)\Bigl{\\{}\widehat{\tau}_{a,k}(V_{i})-\mu_{a,T}(V_{i})\Bigr{\\}}\biggr{]}\right)$
$\displaystyle=\underbrace{E\left[\frac{n}{n_{k}}\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}\frac{I(A_{i}=a,R_{i}=k)}{\widehat{\pi}_{a,k}(X_{i})}-I(R_{i}=k)\Bigr{\\}}\Bigl{\\{}m_{a,k}(X_{i})-\widehat{m}_{a,k}(X_{i})\Bigr{\\}}\right]}_{T_{1}}$
$\displaystyle+\underbrace{E\left[\Bigl{\\{}\frac{n}{n_{T}}I(R_{i}=T)-\frac{n}{n_{k}}I(R_{i}=k)\widehat{\zeta}_{k}(V_{i})\Bigr{\\}}\Bigl{\\{}\widehat{\tau}_{a,k}(V_{i})-\tau_{a,T}(V_{i})\Bigr{\\}}\right]}_{T_{2}}.$
(40)

Step 3 is derived from Step 2 through conditioning on the variables
$(V_{i},R_{i}=k)$, and we obtain the following equality:

$E\left(\frac{n}{n_{k}}\biggl{[}I(R_{i}=k)\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}\tilde{\tau}_{a,k}(V_{i})-\tau_{a,k}(V_{i})\Bigr{\\}}\biggr{]}\right)=E\left(\frac{n}{n_{k}}\biggl{[}I(R_{i}=k)\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}\widehat{m}_{a,k}(V_{i})-m_{a,k}(V_{i})\Bigr{\\}}\biggr{]}\right).$
(41)

Step 4 can be derived from Step 3 by definition; in the target site, the
conditional outcome is represented by $\mu_{a,T}={m}_{a,T}(V_{i})$ and
$\tau_{a,T}=E\left\\{{m}_{a,T}(V_{i})\mid
V_{i},R_{i}=T\right\\}={m}_{a,T}(V_{i})=\mu_{a,T}$. Conditioning on $X_{i}$,
we have the following expression:

$\displaystyle T_{1}=$ $\displaystyle
E\left[\frac{n}{n_{k}}\widehat{\zeta}_{k}(V_{i})\Bigl{\\{}\frac{I(A_{i}=a,R_{i}=k)}{\widehat{\pi}_{a,k}(X_{i})}-I(R_{i}=k)\Bigr{\\}}\Bigl{\\{}m_{a,k}(X_{i})-\widehat{m}_{a,k}(X_{i})\Bigr{\\}}\right]$
$\displaystyle=$ $\displaystyle
E\left[\frac{n}{n_{k}}\widehat{\zeta}_{k}(V_{i})E\Bigl{\\{}\frac{I(A_{i}=a,R_{i}=k)}{\widehat{\pi}_{a,k}(X_{i})}-I(R_{i}=k)\mid
X_{i}\Bigr{\\}}\Bigl{\\{}m_{a,k}(X_{i})-\widehat{m}_{a,k}(X_{i})\Bigr{\\}}\right]$
$\displaystyle=$ $\displaystyle
E\left[\frac{n}{n_{k}}\widehat{\zeta}_{k}(V_{i})P(R_{i}=k\mid
X_{i})\Bigl{\\{}\frac{{\pi}_{a,k}(X_{i})}{\widehat{\pi}_{a,k}(X_{i})}-1\Bigr{\\}}\Bigl{\\{}m_{a,k}(X_{i})-\widehat{m}_{a,k}(X_{i})\Bigr{\\}}\right]$
$\displaystyle=$ $\displaystyle
E\left[\frac{\widehat{\zeta}_{k}(V_{i})}{\widehat{\pi}_{a,k}(X_{i})}\Bigl{\\{}{\pi}_{a,k}(X_{i})-\widehat{\pi}_{a,k}(X_{i})\Bigr{\\}}\Bigl{\\{}m_{a,k}(X_{i})-\widehat{m}_{a,k}(X_{i})\Bigr{\\}}\right]$
(42)

Given the boundedness of $\widehat{\zeta}_{k}$ and $\widehat{\pi}_{a,k}$, it
can be concluded that
$T_{1}=O_{p}(\|\widehat{\pi}_{a,k}-\pi_{a,k}\|\|\widehat{m}_{a,k}-{m}_{a,k}\|)$.
Now let us consider $T_{2}$. Conditioning on $V_{i}$, we have the following
expression:

$\displaystyle T_{2}$
$\displaystyle=E\left[\Bigl{\\{}\frac{n}{n_{T}}I(R_{i}=T)-\frac{n}{n_{k}}I(R_{i}=k)\widehat{\zeta}_{k}(V_{i})\Bigr{\\}}\Bigl{\\{}\widehat{\tau}_{a,k}(V_{i})-\tau_{a,T}(V_{i})\Bigr{\\}}\right]$
$\displaystyle=E\left[E\Bigl{\\{}\frac{n}{n_{T}}I(R_{i}=T)-\frac{n}{n_{k}}I(R_{i}=k)\widehat{\zeta}_{k}(V_{i})\mid
V_{i}\Bigr{\\}}\Bigl{\\{}\widehat{\tau}_{a,k}(V_{i})-\tau_{a,T}(V_{i})\Bigr{\\}}\right]$
$\displaystyle=E\left[\Bigl{\\{}\frac{n}{n_{T}}P(R_{i}=T\mid
V_{i})-\frac{n}{n_{k}}P(R_{i}=k\mid
V_{i})\widehat{\zeta}_{k}(V_{i})\Bigr{\\}}\Bigl{\\{}\widehat{\tau}_{a,k}(V_{i})-\tau_{a,T}(V_{i})\Bigr{\\}}\right]$
$\displaystyle=E\left[\Bigl{\\{}\frac{1}{P(R_{i}=T)}P(R_{i}=T\mid
V_{i})-\frac{1}{P(R_{i}=k)}P(R_{i}=k\mid
V_{i})\widehat{\zeta}_{k}(V_{i})\Bigr{\\}}\Bigl{\\{}\widehat{\tau}_{a,k}(V_{i})-\tau_{a,T}(V_{i})\Bigr{\\}}\right]$
$\displaystyle=E\left[\frac{P(R_{i}=k\mid
V_{i})}{P(R_{i}=k)}\Bigl{\\{}{\zeta}_{k}(V_{i})-\widehat{\zeta}_{k}(V_{i})\Bigr{\\}}\Bigl{\\{}\widehat{\tau}_{a,k}(V_{i})-\tau_{a,T}(V_{i})\Bigr{\\}}\right].$
(43)

Therefore,
$T_{2}=O_{p}(\|\widehat{\zeta}_{k}-\zeta_{k}\|\left\|\widehat{\tau}_{a,k}-\tau_{a,k}\right\|)$.

Combining all results from the above three steps yields

$\displaystyle\|\widehat{\mu}_{a,k}-\mu_{a,T}\|=O_{p}\left(n^{-1/2}+\|\widehat{\pi}_{a,k}-\pi_{a,k}\|\|\widehat{m}_{a,k}-{m}_{a,k}\|+\|\widehat{\zeta}_{k}-{\zeta}_{k}\|\|\widehat{\tau}_{a,k}-\tau_{a,k}\|\right).$
(44)

Further, if the nuisance estimators satisfy the following convergence rate

$\displaystyle\left\|\widehat{m}_{a,k}-m_{a,k}\right\|\left\|\widehat{\pi}_{a,k}-\pi_{a,k}\right\|$
$\displaystyle=o_{p}(1/\sqrt{n}),\quad\|\widehat{\zeta}_{k}-\zeta_{k}\|\left\|\widehat{\tau}_{a,k}-\tau_{a,k}\right\|=o_{p}(1/\sqrt{n}),$
(45)

then the conditional bias,
$E\left(\widehat{\mu}_{a,k}-\mu_{a,T}\right)=o_{p}(1/\sqrt{n})$ since
$O_{p}(o_{p}(1/\sqrt{n}))=o_{p}(1/\sqrt{n})$.

Consequently, by Slutsky’s theorem,
$\sqrt{n}(\widehat{\mu}_{a,k}-\mu_{a,T})\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,\sigma^{2})$. Here, $\sigma^{2}$ represents
the semiparametric efficiency bound, which has been derived and defined in
(39). ∎

###  D.5 Proof of Theorem 3

###### Proof.

The proof presented herein closely follows the methodology outlined in [12];
we start by establishing the consistency and asymptotic normality of the
global estimator, assuming a fixed ${\eta}_{k}$. We then invoke Lemma 4 and 5
in [12], which state that the proposed adaptive estimation for $\eta_{k}$ as
shown in (9) allows for (i) the recovery of the optimal $\bar{\eta}_{k}$ by
the estimator $\widehat{\eta}_{k}$, and (ii) the uncertainty introduced by
$\widehat{\eta}_{k}$ is negligible when estimating $\Delta_{T}$.

Lemma 3 demonstrates the consistency of the site-specific estimators given
that the source site $k$ satisfies the conditions outlined in Theorem 2. We
denote the set of source sites that fulfill the conditions in Theorem 2 as
$\mathcal{S}^{*}$, and consider a fixed ${\eta}_{k}$ such that ${\eta}_{k}=0$
for $k\not\in\mathcal{S}^{*}$, then

$\displaystyle\widehat{\mu}_{a,G}(\eta_{k})=\widehat{\mu}_{a,T}+\sum_{k\in\mathcal{K}}{\eta}_{k}\\{\widehat{\mu}_{a,k}-\widehat{\mu}_{a,T}\\}$
(46)

is consistent for $\mu_{a,T}$ since $\widehat{\mu}_{a,k}$ for
$k\in\mathcal{S}^{*}$ are consistent estimators for ${\mu}_{a,T}$. Hence, we
can establish that
$\widehat{\Delta}_{G}(\eta_{k})=\widehat{\mu}_{1,G}(\eta_{k})-\widehat{\mu}_{0,G}(\eta_{k})$
consistently estimates $\Delta_{T}=\mu_{1,T}-\mu_{0,T}$. Moving forward, we
proceed to examine the asymptotic normality of the global estimator by
utilizing the influence functions for the site-specific estimators. First, we
rewrite the global estimator with the fixed $\eta_{k}$ as

$\displaystyle\widehat{\Delta}_{G}(\eta_{k})=\left(1-\sum_{k\in\mathcal{S}}{\eta}_{k}\right)\left(\widehat{\mu}_{1,T}-\widehat{\mu}_{0,T}\right)+\sum_{k\in\mathcal{S}}{\eta}_{k}\left(\widehat{\mu}_{1,k}-\widehat{\mu}_{0,k}\right)=\left(1-\sum_{k\in\mathcal{S}}{\eta}_{k}\right)\widehat{\Delta}_{T}+\sum_{k\in\mathcal{S}}{\eta}_{k}\widehat{\Delta}_{k}$
(47)

In Appendix D.3, the influence functions for the source site and target site
estimators have been derived. To facilitate representation, we decompose the
influence functions for the source site estimators into two parts, each
defined on the target sample and source sample, respectively,

$\displaystyle\xi^{(1)}_{a,k}(Z_{i})$
$\displaystyle=\frac{1}{P(R_{i}=k)}\Biggl{[}\frac{I(A_{i}=a,R_{i}=k)}{\pi_{a,k}(X_{i})}\zeta_{k}(V_{i})\\{Y_{i}-m_{a,k}(X_{i})\\}\Biggr{]}$
(48)
$\displaystyle+\frac{1}{P(R_{i}=k)}\Biggl{[}I(R_{i}=k)\zeta_{k}(V_{i})\\{m_{a,k}(X_{i})-\tau_{a,k}(V_{i})\\}\Biggr{]}$
$\displaystyle\xi^{(2)}_{a,k}(Z_{i})$
$\displaystyle=\frac{1}{P(R_{i}=T)}\Biggl{[}I(R_{i}=T)\tau_{a,k}(V_{i})\Biggr{]}-\mu_{a,T}.$

Here, it should be noted that the expectations of where
$\xi^{(1)}_{a,k}(Z_{i})$ and $\xi^{(2)}_{a,k}(Z_{i})$ are both equal to zero.
Furthermore, for the target site estimator, we define:

$\displaystyle\xi^{(2)}_{a,T}(Z_{i})$
$\displaystyle=\frac{1}{P(R_{i}=T)}\Biggl{[}\frac{I(A_{i}=a,R_{i}=T)}{\pi_{a,T}(V_{i})}\\{Y_{i}-m_{a,T}(V_{i})\\}+I(R_{i}=T)m_{a,T}(V_{i})\Biggr{]}-\mu_{a,T}$
(49)

and the expectation of $\xi^{(2)}_{a,T}(Z_{i})$ is also zero. We denote the
influence function for $\widehat{\Delta}_{k}$ as
$\xi^{(1)}_{1,k}(Z_{i})-\xi^{(1)}_{0,k}(Z_{i})+\xi^{(2)}_{1,k}(Z_{i})-\xi^{(2)}_{0,k}(Z_{i})$.
Similarly, the influence function for $\widehat{\Delta}_{T}$ is given by
$\xi^{(2)}_{1,T}(Z_{i})-\xi^{(2)}_{0,T}(Z_{i})$. Thus, the influence function
for the global estimator evaluated with target and source samples can be
expressed as:

$\displaystyle\widehat{\Delta}_{G}(\eta_{k})-\Delta_{T}$
$\displaystyle=\left(1-\sum_{k\in\mathcal{S}}{\eta}_{k}\right)\left(\widehat{\Delta}_{T}-\Delta_{T}\right)+\sum_{k\in\mathcal{S}}{\eta}_{k}\left(\widehat{\Delta}_{k}-\Delta_{T}\right)$
$\displaystyle=\left(1-\sum_{k\in\mathcal{S}}{\eta}_{k}\right)\frac{1}{n_{T}}\sum_{i=1}^{N}I(R_{i}=T)\left\\{\xi^{(2)}_{1,T}(Z_{i})-\xi^{(2)}_{0,T}(Z_{i})\right\\}$
$\displaystyle+\left(\sum_{k\in\mathcal{S}}{\eta}_{k}\right)\frac{1}{n_{T}}\sum_{i=1}^{N}I(R_{i}=T)\left\\{\xi^{(2)}_{1,k}(Z_{i})-\xi^{(2)}_{0,k}(Z_{i})\right\\}$
$\displaystyle+\sum_{k\in\mathcal{S}}{\eta}_{k}\frac{1}{n_{k}}\sum_{i=1}^{N}I(R_{i}=k)\left\\{\xi^{(1)}_{1,k}(Z_{i})-\xi^{(1)}_{0,k}(Z_{i})\right\\}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}I(R_{i}=T)\left(1-\sum_{k\in\mathcal{S}}{\eta}_{k}\right)\frac{\left\\{\xi^{(2)}_{1,T}(Z_{i})-\xi^{(2)}_{0,T}(Z_{i})\right\\}}{P(R_{i}=T)}$
$\displaystyle+\frac{1}{N}\sum_{i=1}^{N}I(R_{i}=T)\left(\sum_{k\in\mathcal{S}}{\eta}_{k}\right)\frac{\left\\{\xi^{(2)}_{1,k}(Z_{i})-\xi^{(2)}_{0,k}(Z_{i})\right\\}}{P(R_{i}=T)}$
$\displaystyle+\frac{1}{N}\sum_{k\in\mathcal{S}}\sum_{i=1}^{N}I(R_{i}=k){\eta}_{k}\frac{\left\\{\xi^{(1)}_{1,k}(Z_{i})-\xi^{(1)}_{0,k}(Z_{i})\right\\}}{P(R_{i}=k)}$
(50)

The asymptotic variance for $\widehat{\Delta}_{G}(\eta_{k})$ equals the
variance of the influence function (50). Its derivation is similar to that of
39. Let us denote this asymptotic variance as ${\mathcal{V}}(\eta_{k})$. Under
the assumption of i.i.d. units within each site, we have:

$\displaystyle{\mathcal{V}}(\eta_{k})$
$\displaystyle=\left(1-\sum_{k\in\mathcal{S}}{\eta}_{k}\right)^{2}\frac{\text{Var}\left\\{\xi^{(2)}_{1,T}(Z_{i})-\xi^{(2)}_{0,T}(Z_{i})\mid
R_{i}=T\right\\}}{P(R_{i}=T)}$
$\displaystyle+\left(\sum_{k\in\mathcal{S}}{\eta}_{k}\right)^{2}\frac{\text{Var}\left\\{\xi^{(2)}_{1,k}(Z_{i})-\xi^{(2)}_{0,k}(Z_{i})\mid
R_{i}=T\right\\}}{P(R_{i}=T)}$
$\displaystyle+2\left(1-\sum_{k\in\mathcal{S}}{\eta}_{k}\right)\left(\sum_{k\in\mathcal{S}}{\eta}_{k}\right)\frac{\text{Cov}\left\\{\xi^{(2)}_{1,T}(Z_{i})-\xi^{(2)}_{0,T}(Z_{i}),\xi^{(2)}_{1,k}(Z_{i})-\xi^{(2)}_{0,k}(Z_{i})\mid
R_{i}=T\right\\}}{P(R_{i}=T)}$
$\displaystyle+\sum_{k\in\mathcal{S}}{\eta}_{k}^{2}\frac{\text{Var}\left\\{\xi^{(1)}_{1,k}(Z_{i})-\xi^{(1)}_{0,k}(Z_{i})\mid
R_{i}=k\right\\}}{P(R_{i}=k)}$ (51)

Under the uniform boundedness conditions stated in Theorem 2, this variance is
finite. Consequently, we can express the asymptotic distribution of
$\sqrt{N}\left(\widehat{\Delta}{G}(\eta_{k})-\Delta{T}\right)$ as:

$\sqrt{N}\left(\widehat{\Delta}_{G}(\eta_{k})-\Delta_{T}\right)\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}\left(0,{\mathcal{V}}(\eta_{k})\right)$ (52)

We further define the optimal adaptive weight $\bar{\eta}_{k}$ as follows:

$\bar{\eta}_{k}=\arg\min_{{\eta}_{k}=0\forall
k\not\in\mathcal{S}^{*}}\mathcal{V}(\eta_{k})$ (53)

By leveraging Lemmas 4 and 5 from [12], we can recover the optimal
$\bar{\eta}_{k}$ with negligible uncertainty for estimating $\Delta_{T}$ if we
estimate $\widehat{\eta}_{k,L_{1}}$ using (9). The consistency of
$\widehat{\mathcal{V}}(\widehat{\eta}_{k,L_{1}})$ follows when we can
effectively approximate $\mathcal{V}(\bar{\eta}_{k})$ with
$\widehat{\mathcal{V}}(\widehat{\eta}_{k,L_{1}})$. Thus,

$\sqrt{N/\widehat{\mathcal{V}}(\widehat{\eta}_{k,L_{1}})}\left(\widehat{\Delta}_{G}(\widehat{\eta}_{k,L_{1}})-\Delta_{T}\right)\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}\left(0,1\right)$ (54)

We now proceed to analyze the efficiency gain resulting from the federation
process. The estimator that relies solely on the target data is denoted as
$\widehat{\Delta}_{T}=\widehat{\Delta}_{G}(\eta_{0})$, where $\eta_{0}$
assigns all weights to the target and none to the source. In contrast, the
estimator that leverages the proposed adaptive ensemble approach is denoted as
$\widehat{\Delta}_{G}(\widehat{\eta}_{k,L_{1}})=(1-\sum_{k\in\mathcal{S}}\widehat{\eta}_{k,L_{1}})\widehat{\Delta}_{T}+\sum_{k\in\mathcal{S}}\widehat{\eta}_{k,L_{1}}\widehat{\Delta}_{k}$.
Here, $\widehat{\eta}_{k,L_{1}}$ can recover the optimal weight
$\bar{\eta}_{k}$ that is associated with the minimum asymptotic variance.
Consequently, the variance of $\widehat{\Delta}_{G}(\widehat{\eta}_{k,L_{1}})$
is no larger than that of the estimator relying solely on the target data by
definition.

To establish that the asymptotic variance of
$\widehat{\Delta}_{G}(\widehat{\eta}_{k,L_{1}})$ is strictly smaller than that
of the estimator based solely on the target data $\widehat{\Delta}_{T}$, we
adapt Proposition 1 in [12] with a modified informative source condition.
Specifically, for each source site $s\in\mathcal{S}^{*}$, we define
$\widehat{\Delta}_{G}(\eta_{s})$ a global estimator where $\eta_{s}$ is the
optimal ensemble weight if we only consider target site and this source site
$s$. Then the modified informative source condition is
$\operatorname{Cov}\left\\{\sqrt{N}\widehat{\Delta}_{T},\sqrt{N}\left(\widehat{\Delta}_{G}({\eta}_{s})-\widehat{\Delta}_{T}\right)\right\\}$,
where $\widehat{\Delta}_{G}({\eta}_{s})-\widehat{\Delta}_{T}$ can be expressed
as

$\displaystyle\widehat{\Delta}_{G}({\eta}_{s})-\widehat{\Delta}_{T}$
$\displaystyle=\widehat{\Delta}_{G}({\eta}_{s})-\Delta_{T}-\left(\widehat{\Delta}_{T}-\Delta_{T}\right)$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}I(R_{i}=T)\left(1-{\eta}_{s}\right)\frac{\left\\{\xi^{(2)}_{1,T}(Z_{i})-\xi^{(2)}_{0,T}(Z_{i})\right\\}}{P(R_{i}=T)}$
$\displaystyle+\frac{1}{N}\sum_{i=1}^{N}I(R_{i}=T){\eta}_{s}\frac{\left\\{\xi^{(2)}_{1,s}(Z_{i})-\xi^{(2)}_{0,s}(Z_{i})\right\\}}{P(R_{i}=T)}$
$\displaystyle+\frac{1}{N}\sum_{i=1}^{N}I(R_{i}=s){\eta}_{s}\frac{\left\\{\xi^{(1)}_{1,s}(Z_{i})-\xi^{(1)}_{0,s}(Z_{i})\right\\}}{P(R_{i}=s)}$
$\displaystyle-\frac{1}{N}\sum_{i=1}^{N}I(R_{i}=T)\frac{\left\\{\xi^{(2)}_{1,T}(Z_{i})-\xi^{(2)}_{0,T}(Z_{i})\right\\}}{P(R_{i}=T)}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}I(R_{i}=T){\eta}_{s}\frac{\left\\{\xi^{(2)}_{1,s}(Z_{i})-\xi^{(2)}_{0,s}(Z_{i})-\xi^{(2)}_{1,T}(Z_{i})+\xi^{(2)}_{0,T}(Z_{i})\right\\}}{P(R_{i}=T)}$
$\displaystyle+\frac{1}{N}\sum_{i=1}^{N}I(R_{i}=s){\eta}_{s}\frac{\left\\{\xi^{(1)}_{1,s}(Z_{i})-\xi^{(1)}_{0,s}(Z_{i})\right\\}}{P(R_{i}=s)}.$
(55)

In summary, if there exist a source site $s$ with consistent estimator of
$\Delta_{T}$ and further satisfy
$\left|\operatorname{Cov}\left\\{\sqrt{N}\widehat{\Delta}_{T},\sqrt{N}\left(\widehat{\Delta}_{G}({\eta}_{s})-\widehat{\Delta}_{T}\right)\right\\}\right|\geq\varepsilon$
where $\varepsilon$ denotes a positive constant, the variance of
$\widehat{\Delta}_{G}$ is strictly smaller than $\widehat{\Delta}_{T}$. ∎