# Federated Adaptive Causal Estimation (FACE) of Target Treatment Effects

Larry Han1, Jue Hou1, Kelly Cho2, Rui Duan1†, Tianxi Cai1,3†
1 Department of Biostatistics, Harvard T.H. Chan School of Public Health
2 Massachusetts Veterans Epidemiology Research and Information Center,
US Department of Veteran Affairs
3 Department of Biomedical Informatics, Harvard Medical School
${\dagger}$ Co-corresponding authors

###### Abstract

Federated learning of causal estimands may greatly improve estimation
efficiency by aggregating estimates from multiple study sites, but robustness
to extreme estimates is vital for maintaining consistency. We develop a
federated adaptive causal estimation (FACE) framework to incorporate
heterogeneous data from multiple sites to provide treatment effect estimation
and inference for a target population of interest. Our strategy is
communication-efficient and privacy-preserving and allows for flexibility in
the specification of the target population. Our method accounts for site-level
heterogeneity in the distribution of covariates through density ratio
weighting. To safely aggregate estimates from all sites and avoid negative
transfer, we introduce an adaptive procedure of weighing the estimators
constructed using data from the target and source populations through a
penalized regression on the influence functions, which achieves 1) consistency
and 2) optimal efficiency. We illustrate FACE by conducting a comparative
effectiveness study of BNT162b2 (Pfizer) and mRNA-1273 (Moderna) vaccines on
COVID-19 outcomes in U.S. veterans using electronic health records from five
VA sites.

Keywords: Adaptive weighting, Doubly robust, Causal inference, Federated
learning, Influence function

##  1 Introduction

Multi-center, federated causal inference is of great interest, particularly
when studying novel treatments, rare diseases, or in times of urgent health
crises. For example, the COVID-19 pandemic has highlighted the need for novel
approaches to efficiently and safely evaluate the efficacy of novel therapies
and vaccines, leveraging data from multiple healthcare systems to ensure the
generalizability of findings. Over the past few years, many research networks
and data consortia have been built to facilitate multi-site studies and have
been actively contributing to COVID-19 studies, including the Observational
Health Data Sciences and Informatics (OHDSI) consortium (Hripcsak et al.,
2016) and the Consortium for Clinical Characterization of COVID-19 by EHR
(Brat et al., 2020).

Analyzing data collected from multiple healthcare systems, however, is highly
challenging for several reasons. First, various sources of heterogeneity exist
in terms of differences in the underlying population of each data set and
policy level variations of treatment assignment. Since the treatment effect
may differ across different patient populations, it would be of interest to
infer the causal treatment effect for a specific target population. However,
the presence of heterogeneity and potential model misspecification poses great
difficulty to ensure valid estimates for the target average treatment effect
(ATE). Furthermore, patient level data typically cannot be shared across
healthcare centers, posing an additional challenge. To overcome these
challenges, we propose a federated causal inference framework to estimate the
ATE for the target population, while leveraging heterogeneous data from other
source populations, and accounting for data sharing constraints though a
privacy-preserving and communication-efficient federated learning strategy.

Most existing literature on federated learning has been focused on regression
and classification models (Chen et al., 2006; Li et al., 2013; Chen & Xie,
2014; Lee et al., 2017; Lian & Fan, 2017; Wang et al., 2019; Duan et al.,
2019). Limited federated learning methods currently exist to make causal
inference with multiple heterogeneous studies. Recently, Xiong et al. (2021)
proposed federated inverse probability weighted (IPW) estimation of the ATE,
where the estimand is the ATE for the entire population rather than a target
population. The ATE considered in Xiong et al. (2021) is a special case of our
setting when the target population is set to be the entire population.
Although Xiong et al. (2021) provided multiple methods for point estimation
and variance estimation, the choice of the proper method depends on prior
knowledge about model homogeneity and model specification, which are usually
unknown in practice. In addition, no empirical study in Xiong et al. (2021)
was provided to test the robustness of the approach to the covariate shift
assumption. Vo et al. (2021) proposed a Bayesian approach that models
potential outcomes as random functions distributed by Gaussian processes.
Their focus is also on the population ATE rather than a particular target
population, and their approach requires specifying parameters and
hyperparameters of Gaussian processes and modeling correlations through kernel
functions, which may be numerically intensive. Compared to these approaches,
our approach estimates the ATE in a particular target population, and accounts
for the potential heterogeneity in the covariate distributions across
populations. Most importantly, our approach additionally safeguards against
the case where incorporating source data can lower the performance of the
method compared to using target data alone, known as negative transfer (Pan &
Yang, 2009).

Another related strand of literature concerns the generalizability and
transportability of randomized clinical trials to EHR studies. For example,
Stuart et al. (2011, 2015, 2018) assessed the generalizability of results from
randomized trials to target popoulations of interest. Dahabreh et al. (2020);
Dong et al. (2020), and Josey et al. (2020) all focused on extending
inferences about treatments from a randomized trial to a new target population
by using different weighting schemes. For a comprehensive review of
statistical methods for generalizability and transportability, see Degtiar &
Rose (2021). However, to date, no literature in generalizability and
transportability has sought to leverage observational data from a potentially
large number of source sites in a data-adaptive manner to obtain unbiased,
efficient, and robust estimation of target treatment effects.

To incorporate heterogeneous data from multiple sites to make inference about
the target average treatment effect (TATE), we develop a federated adaptive
causal estimation (FACE) framework. Our major contributions can be summarized
as follows. First, our framework allows for flexibility in the specification
of the target population (e.g., a single site, some sites, or all sites).
Second, our method accounts for site-level heterogeneity in the distribution
of covariates through a density ratio weighting approach, such that each
source site can obtain a doubly robust estimate of the TATE. Third, to safely
aggregate estimates from all sites and avoid negative transfer, we introduce
an adaptive method that anchors on the target population estimate and then
estimates data-adaptive weights to combine source site estimates. Compared to
the estimator using only target population data, our proposed global estimator
adaptively combines data from all populations to meaningfully increase
estimation precision without introducing bias. In addition, we propose a
communication-efficient and privacy-preserving federated algorithm to obtain
the doubly robust estimator, which only requires each participating site to
share information once.

The remainder of the paper is organized as follows. In Section 2, we introduce
the problem setting, notation, and assumptions required for identification of
the TATE. In Section 3, we declare the proposed FACE framework for estimating
the TATE. We introduce the in-site estimators based on the target population
and source population separately in Sections 3.1 and 3.2 and present the
adaptive and distributed integration in Section 3.3. In Section 4, we provide
the theoretical guarantees of the method, including double robustness and
asymptotic normality. In Section 5, we conduct extensive simulations for
various numbers of sites, sample sizes, and show robustness to
misspecification of different models. In Section 6, we apply FACE to conduct a
comparative effectiveness study of COVID-19 vaccines using the EHRs from five
federated VA sites. We conclude in Section 7 with key takeaways and directions
for future research.

##  2 Setting and Notation

For the $i$-th observation, we denote the outcome as $Y_{i}\in\mathbb{R}$, the
$p$-dimensional baseline covariate vector as
$\mathbf{X}_{i}=(X_{i1},...,X_{ip})^{\top}\in\mathcal{X}\subset\mathbb{R}^{p}$,
and the indicator for binary treatment as $A_{i}\in\\{0,1\\}$. Under the
federated learning setting, a total of $N$ observations are stored at $K$
study sites, where the $k$-th site has sample size $n_{k}$, and
$N=\sum_{k=1}^{K}n_{k}$. Let $R_{i}$ be a site indicator such that $R_{i}=k$
indicates the $i$-th patient in the $k$-th site. Indexing the site by a single
integer $R_{i}$, we assume that each observation may only belong to one site.
We summarize the observed data at each site $k$ as
$\mathscr{D}_{k}=\\{(Y_{i},\mathbf{X}_{i}^{{}^{\top}},A_{i},R_{i})^{\top},R_{i}=k\\},$
where each site has access to its own patient-level data but cannot share this
data with other sites. We denote the index set for each site as
$\mathcal{I}_{k}=\\{i:R_{i}=k\\}$ and the data in the target sites as
$\mathscr{D}_{\mathcal{T}}$. Let $\mathcal{T}$ be a subset of
$\\{1,\dots,K\\}$, indicating sites that are in the target population, and
$\mathcal{T}^{c}=\\{1,\dots,K\\}\setminus\mathcal{T}$ containing sites that
are in the source population. We denote the total sample size of the target
sites as $N_{\mathcal{T}}=\sum_{k\in\mathcal{T}}n_{k}$. For simplicity of
notation, we use $(Y,\mathbf{X},A,R)$ without subscripts to state general
assumptions and conclusions.

Under the potential outcomes framework (Neyman, 1923; Rubin, 1974), we denote
$Y^{(a)}$ as the potential outcome of patient under treatment $A=a$, $a=0,1$.
Our goal is to estimate the target average treatment effect (TATE) for a
specified target population,

$\Delta_{\mathcal{T}}=\mu^{(1)}_{\mathcal{T}}-\mu^{(0)}_{\mathcal{T}},\quad\mu^{(a)}_{\mathcal{T}}=\mathbb{E}(Y^{(a)}\mid
R\in\mathcal{T}),$ (1)

where the expectation is taken over the covariate distribution in the target
population. The target population can be specified at multiple levels
corresponding to different targets of real-world interest. When the target of
interest includes patients across all sites, i.e.,
$\mathcal{T}=\\{1,...,K\\}$, then the TATE is simply the ATE. This corresponds
to the special case that is studied by Xiong et al. (2021). When the target
population includes some of the sites, leveraging source sites is useful as it
can allow for greater efficiency in the estimation of the TATE.

To identify the TATE, we make the following standard assumptions throughout
the paper:

######  Assumption 1.

For an absolute positive constant $\varepsilon>0$, $a\in\\{0,1\\}$, and
$\mathbf{x}\in\mathcal{X}$,

  1. (a)

Consistency: $Y=Y^{(A)}$.

  2. (b)

Overlapping of treatment arms:
$\mathbb{P}(A=a\mid\mathbf{X}=\mathbf{x},R=k)\in(\varepsilon,1-\varepsilon)$,
$k\in\\{1,\dots,K\\}$.

  3. (c)

Overlapping of site populations:
$\mathbb{P}(R=k\mid\mathbf{X}=\mathbf{x})>\varepsilon$, $k\in\\{1,\dots,K\\}$.

  4. (d)

Ignorability: $\left(Y^{(1)},Y^{(0)}\right)\perp\\!\\!\\!\perp
A\mid\mathbf{X},R=k$, $k\in\mathcal{T}$.

Assumption 1(a) states that the observed outcome for patient $i$ under
treatment $a$ is the patient’s potential outcome under the same treatment.
Assumptions 1(b) and 1(c) are the standard overlap assumptions (Rosenbaum &
Rubin, 1983), which state that the probability of being assigned to each
treatment/site, conditional on baseline covariates, is positive in each site.
Here we allow each site to have a site-specific propensity score. Assumption
1(d) can be derived from the typical ignorability condition (Imbens & Rubin,
2015; Hernán & Robins, 2020). For further discussion on assumptions for
identification, see Dahabreh et al. (2020).

######  Remark 1.

Note that Assumption 1(d) is only made for target sites $k\in\mathcal{T}$. We
make no stringent requirement for $(Y^{(1)},Y^{(0)},A)\mid\mathbf{X}$ among
the source sites $k\in\mathcal{T}^{c}$. In the worst case scenario, data from
a source site can be totally irrelevant to $\Delta_{\mathcal{T}}$ such that
the estimator from the site is severely biased. To prevent negative transfer
from such cases, we devise an adaptive selection and aggregation for the
source sites.

We denote the models for the site-specific propensity scores and outcome
regressions:

$\displaystyle\text{Propensity score }\mathbb{P}(A=a\mid
R=k,\mathbf{X}):\pi_{k}(a,\mathbf{X};\boldsymbol{\alpha}_{k}),$
$\displaystyle\text{Outcome regression }\mathbb{P}(Y\mid
R=k,A=a,\mathbf{X}):m_{k}(a,\mathbf{X};\boldsymbol{\beta}_{a,k}).$ (2)

Since the distribution of the covariates $\mathbf{X}$ may be heterogeneous
across sites, we characterize the difference in covariate distributions
between the target site(s) and the $k$-th source site through the density
ratio

$\omega_{k}(\mathbf{x})=\frac{f(\mathbf{X}\mid
R\in\mathcal{T})}{f(\mathbf{X}\mid
R=k)}=\frac{\mathbb{P}(R\in\mathcal{T}\mid\mathbf{X}=\mathbf{x})\mathbb{P}(R=k)}{\mathbb{P}(R=k\mid\mathbf{X}=\mathbf{x})\mathbb{P}(R\in\mathcal{T})}.$

We choose flexible semiparametric models for the density ratio

$\omega_{k}(\mathbf{X};\boldsymbol{\gamma}_{k})=\exp(\boldsymbol{\gamma}_{k}^{\top}\boldsymbol{\psi}(\mathbf{X})),$
(3)

where $\boldsymbol{\psi}$ is a $\mathbb{R}^{p}\mapsto\mathbb{R}^{q}$ function
with an intercept term.

######  Remark 2.

The model (3) is known as the exponential tilt density ratio model, which has
wide application in areas such as mixture models to account for heterogeneity
between two distributions (Qin, 1998; Qin & Liang, 2011; Hong et al., 2017;
Duan et al., 2020b). For example, Qin & Liang (2011) studied the exponential
tilt model with $\boldsymbol{\psi}(\mathbf{x})=\mathbf{x}$, which recovers the
whole class of natural exponential family distributions, including the normal
distribution with mean shift, Bernoulli distribution for binary covariates,
etc. By including higher-order terms of $\mathbf{x}$ in
$\boldsymbol{\psi}(\mathbf{x})$, higher-order differences such as variance,
skewness, and covariance can be captured. An important contribution in our
paper is that we propose a privacy-preserving and communication-efficient
approach to estimate the heterogeneity in covariate distributions, which we
will describe in more detail in Section 3.

##  3 Method

Staring with an overview, we illustrate the workflow of the FACE framework in
Figure 1. In step 1, an estimate $\widehat{\Delta}_{\mathcal{T},k}$ will be
locally constructed within the $k$-th target site, and these estimate(s) will
be shared with a leading site. In the meantime, the $k$-th target site will
also calculate summary statistics
$\overline{\boldsymbol{\psi}}(\mathbf{x}_{k})$ and broadcast them to the
source sites, which will later be used for estimation of the density ratio
model. In Step 2, each source site will construct a site-specific estimate for
the TATE, $\widehat{\Delta}_{\mathcal{T},k}$, and pass it to the leading site.
In Step 3, a leading site processes the aggregation with estimators and
parameters from Steps 1 and 2 to obtain the final FACE estimator,
$\widehat{\Delta}_{\mathcal{T},\rm{FACE}}$. Overall, each site is only
required to share information one time with other sites. We provide detailed
explanations for each step of FACE, including the summary-level information
that must be shared across sites, in Sections 3.1-3.3 with generic models. An
example with logistic regression models is given in Section 3.4.

Figure 1: Workflow of FACE to construct a global estimator in a federated
setting

###  3.1 Step 1: Estimation Using Target Data

The initial doubly robust TATE estimator is obtained from the site-specific
ATE of the target site(s). For each target site $k\in\mathcal{T}$, we estimate
the in-site propensity score and outcome models
$\pi_{k}(a,\mathbf{X};\widehat{\boldsymbol{\alpha}}_{k})$ and
$m_{k}(a,\mathbf{X};\widehat{\boldsymbol{\beta}}_{a,k})$ using classical
regression models and methods such as maximum likelihood estimation or
estimating equations.

We compute the in-site contribution to the doubly robust TATE,
$\widehat{\Delta}_{\mathcal{T},k}$ from its plug-in component
$\widehat{M}_{\mathcal{T},k}$ and augmentation component
$\widehat{\delta}_{\mathcal{T},k}$, as well as summary characteristics for the
covariate distribution from the target sites,
$\overline{\boldsymbol{\psi}}_{\mathcal{T},k}$:

$\displaystyle\widehat{M}_{\mathcal{T},k}$
$\displaystyle=n_{k}^{-1}\sum_{i\in\mathcal{I}_{k}}\\{m_{k}(1,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{1,k})-m_{k}(0,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{0,k}),$
$\displaystyle\widehat{\delta}_{\mathcal{T},k}$
$\displaystyle=n_{k}^{-1}\sum_{i\in\mathcal{I}_{k}}\left[\frac{\operatorname*{\mathrm{I}}(A_{i}=1)}{\pi_{k}(1,\mathbf{X};\widehat{\boldsymbol{\alpha}}_{k})}\\{Y_{i}-m_{k}(1,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{1,k})\\}-\frac{\operatorname*{\mathrm{I}}(A_{i}=0)}{\pi_{k}(0,\mathbf{X};\widehat{\boldsymbol{\alpha}}_{k})}\\{Y_{i}-m_{k}(0,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{0,k})\\}\right],$
$\displaystyle\widehat{\Delta}_{\mathcal{T},k}$
$\displaystyle=\widehat{M}_{\mathcal{T},k}-\widehat{\delta}_{\mathcal{T},k},\qquad\overline{\boldsymbol{\psi}}_{\mathcal{T},k}=n_{k}^{-1}\sum_{i\in\mathcal{I}_{k}}\boldsymbol{\psi}(\mathbf{X}_{i}).$
(4)

Once $\widehat{M}_{\mathcal{T},k}$, $\widehat{\delta}_{\mathcal{T},k}$,
$\widehat{\Delta}_{\mathcal{T},k}$,
$\overline{\boldsymbol{\psi}}_{\mathcal{T},k}$, and $n_{k}$ for
$k\in\mathcal{T}$ are shared, each site can compute the doubly robust TATE,
its augmentation component, and the mean covariate vector from the target
sites:

$\widehat{M}_{\mathcal{T},\mathcal{T}}=\frac{1}{N_{\mathcal{T}}}\sum_{k\in\mathcal{T}}n_{k}\widehat{M}_{\mathcal{T},k},\;\widehat{\delta}_{\mathcal{T},\mathcal{T}}=\frac{1}{N_{\mathcal{T}}}\sum_{k\in\mathcal{T}}n_{k}\widehat{\delta}_{\mathcal{T},k},\;\widehat{\Delta}_{\mathcal{T},\mathcal{T}}=\frac{1}{N_{\mathcal{T}}}\sum_{k\in\mathcal{T}}n_{k}\widehat{\Delta}_{\mathcal{T},k},\;\overline{\boldsymbol{\psi}}_{\mathcal{T}}=\frac{1}{N_{\mathcal{T}}}\sum_{k\in\mathcal{T}}n_{k}\overline{\boldsymbol{\psi}}_{\mathcal{T},k}.$

When $\mathcal{T}$ contains multiple sites, we also share the estimators for
the variance-covariance of
$(\widehat{M}_{\mathcal{T},k},\widehat{\delta}_{\mathcal{T},k},\overline{\boldsymbol{\psi}}_{\mathcal{T},k})$,
which we denote as $\widehat{\Sigma}_{k}$. Variance estimation
$\widehat{\Sigma}_{k}$ can be conducted through classical influence function
or bootstrap within site. The role of the matrix in the aggregation will be
explained in (10).

######  Remark 3.

The initial TATE estimator $\widehat{\Delta}_{\mathcal{T},\mathcal{T}}$ is the
doubly robust estimator with site-stratified propensity score and outcome
regression:

$\displaystyle\widehat{\Delta}_{\mathcal{T},\mathcal{T}}=N_{\mathcal{T}}^{-1}\sum_{k\in\mathcal{T}}\sum_{i\in\mathcal{I}_{k}}$
$\displaystyle\bigg{[}\\{m_{k}(1,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{1,k})-m_{k}(0,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{0,k})+\frac{\operatorname*{\mathrm{I}}(A_{i}=1)}{\pi_{k}(1,\mathbf{X};\widehat{\boldsymbol{\alpha}}_{k})}\\{Y_{i}-m_{k}(1,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{1,k})\\}$
$\displaystyle-\frac{\operatorname*{\mathrm{I}}(A_{i}=0)}{\pi_{k}(0,\mathbf{X};\widehat{\boldsymbol{\alpha}}_{k})}\\{Y_{i}-m_{k}(0,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{0,k})\\}\bigg{]}.$

By Bang & Robins (2005), $\widehat{\Delta}_{\mathcal{T},\mathcal{T}}$ is
consistent if either the propensity score or outcome regression is
consistently estimated for each $k\in\mathcal{T}$. Mix-and-match is allowed,
i.e. the propensity score is consistent for some target sites while the
outcome regression is consistent for the rest. If a common outcome regression
model were assumed, then $\boldsymbol{\beta}_{a}$ could be estimated jointly
for additional efficiency gain.

###  3.2 Step 2: Estimation Using Source Data

To assist in estimating $\Delta_{\mathcal{T}}$ using source data, the
covariate shifts between the source sites and the target population need to be
adjusted for in order to eliminate potential biases. We achieve this goal by
tilting the source population to the target population through the density
ratio model introduced in (3). When all data are pooled together, estimating
$\widehat{\boldsymbol{\gamma}}_{k}$ can be achieved by constructing a pseudo-
likelihood function as described in Qin (1998). In a federated setting, we can
only use $\overline{\boldsymbol{\psi}}_{\mathcal{T}}$ obtained from the target
sites and estimate the density ratio
$\omega_{k}(\mathbf{X};\widehat{\boldsymbol{\gamma}}_{k})$ by solving the
estimating equation

$n_{k}^{-1}\sum_{i\in\mathcal{I}_{k}}\omega_{k}\left(\boldsymbol{\psi}(\mathbf{X}_{i});\boldsymbol{\gamma}_{k}\right)\boldsymbol{\psi}(\mathbf{X}_{i})=\overline{\boldsymbol{\psi}}_{\mathcal{T}}.$
(5)

The above estimation is done completely within the $k$-th source site.
Together with the in-site propensity score
$\pi_{k}(a,\mathbf{X};\widehat{\boldsymbol{\alpha}}_{k})$ and outcome
regression models $m_{k}(a,\mathbf{X};\widehat{\boldsymbol{\beta}}_{a,k})$, we
construct the site-specific augmentation,

$\displaystyle\widehat{\delta}_{\mathcal{T},k}=n_{k}^{-1}\sum_{i\in\mathcal{I}_{k}}\omega_{k}(\mathbf{X}_{i};\widehat{\boldsymbol{\gamma}}_{k})$
$\displaystyle\bigg{[}\frac{\operatorname*{\mathrm{I}}(A_{i}=1)}{\pi_{k}(1,\mathbf{X}_{i};\widehat{\boldsymbol{\alpha}}_{k})}\\{Y_{i}-m_{k}(1,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{a,k})\\}$
(6)
$\displaystyle-\frac{\operatorname*{\mathrm{I}}(A_{i}=0)}{\pi_{k}(0,\mathbf{X}_{i};\widehat{\boldsymbol{\alpha}}_{k})}\\{Y_{i}-m_{k}(0,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{a,k})\\}\bigg{]}.$
(7)

Then, the site-specific augmentations are shared back to the leading site. The
source sites also share $\widehat{\sigma}^{2}_{k}$, an estimate for the
conditional variance $\operatorname*{{\rm
Var}}\left(\widehat{\delta}_{\mathcal{T},k}\mid\mathscr{D}_{\mathcal{T}}\right)$
and $\widehat{\mathbf{d}}_{k}$, an estimate for the partial derivative of
$\widehat{\delta}_{\mathcal{T},k}$ with respect to
$\overline{\boldsymbol{\psi}}$, whose role in the aggregation will be
explained in (10). Both $\widehat{\sigma}^{2}_{k}$ and
$\widehat{\mathbf{d}}_{k}$ can be constructed from classical influence
functions, which we provide in the Supplement. Alternatively,
$\widehat{\sigma}^{2}_{k}$ may be estimated by bootstrapping within site and
$\widehat{\mathbf{d}}_{k}$ may be estimated by numerical derivatives.

######  Remark 4.

Combining the site-specific augmentation $\widehat{\delta}_{\mathcal{T},k}$
with the plug-in component $\widehat{M}_{\mathcal{T},\mathcal{T}}$ of the
initial TATE estimator, we obtain the site-specific doubly robust estimator

$\displaystyle\widehat{\Delta}_{\mathcal{T},k}=$
$\displaystyle\widehat{M}_{\mathcal{T},\mathcal{T}}+\widehat{\delta}_{\mathcal{T},k}$
$\displaystyle=$ $\displaystyle
N_{\mathcal{T}}^{-1}\sum_{k\in\mathcal{T}}\sum_{i\in\mathcal{I}_{k}}\\{m_{k}(1,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{1,k})-m_{k}(0,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{0,k})\\}$
$\displaystyle+n_{k}^{-1}\sum_{i\in\mathcal{I}_{k}}\omega_{k}(\mathbf{X}_{i};\widehat{\boldsymbol{\gamma}}_{k})\bigg{[}\frac{\operatorname*{\mathrm{I}}(A_{i}=1)}{\pi_{k}(1,\mathbf{X}_{i};\widehat{\boldsymbol{\alpha}}_{k})}\\{Y_{i}-m_{k}(1,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{a,k})\\}$
$\displaystyle\hskip
108.405pt-\frac{\operatorname*{\mathrm{I}}(A_{i}=0)}{\pi_{k}(0,\mathbf{X}_{i};\widehat{\boldsymbol{\alpha}}_{k})}\\{Y_{i}-m_{k}(0,\mathbf{X}_{i};\widehat{\boldsymbol{\beta}}_{a,k})\\}\bigg{]}.$

The estimator $\widehat{\Delta}_{\mathcal{T},k}$ is doubly robust in the
following sense: either (1) the outcome regression is consistent for the
target site(s) and source sites, or (2) the propensity score and density ratio
are consistent for source sites while the outcome regression estimator is
shared across all sites.

###  3.3 Step 3: Adaptive Aggregation

In the final step, we obtain our FACE estimator by adaptively aggregating the
initial TATE estimator and site-specific estimators. The leading site can
estimate $\Delta_{\mathcal{T}}$ as

$\widehat{\Delta}_{\mathcal{T},\sf
FACE}=\left(1-\sum_{k\in\mathcal{T}^{c}}\eta_{k}\right)\widehat{\Delta}_{\mathcal{T},\mathcal{T}}+\sum_{k\in\mathcal{T}^{c}}\eta_{k}\widehat{\Delta}_{\mathcal{T},k}=\widehat{M}_{\mathcal{T},\mathcal{T}}+\left(1-\sum_{k\in\mathcal{T}^{c}}\eta_{k}\right)\widehat{\delta}_{\mathcal{T},\mathcal{T}}+\sum_{k\in\mathcal{T}^{c}}\eta_{k}\widehat{\delta}_{\mathcal{T},k}.$
(8)

with the weights $\\{\eta_{k},k\in\mathcal{T}^{c}\\}$ to be estimated in a
data-adaptive fashion to filter out potentially biased site-specific
estimators. The second expression of $\widehat{\Delta}_{\mathcal{T},\sf FACE}$
in (8) shows how the shared parameters from Steps 1 and 2 are used in the
construction.

######  Remark 5.

The FACE estimate $\widehat{\Delta}_{\mathcal{T},\sf FACE}$ is a linear
combination of $\widehat{\Delta}_{\mathcal{T},\mathcal{T}}$ and
$\widehat{\Delta}_{\mathcal{T},k}$. It can also be interpreted as an “anchor
and augmentation” estimator. We may alternatively express (8) as

$\widehat{\Delta}_{\mathcal{T},\sf
FACE}=\widehat{\Delta}_{\mathcal{T},\mathcal{T}}+\sum_{k\in\mathcal{T}^{c}}\eta_{k}\left(\widehat{\Delta}_{\mathcal{T},k}-\widehat{\Delta}_{\mathcal{T},\mathcal{T}}\right).$

The alternative expression makes it clear that we “anchor” on the estimator
from the target site, $\widehat{\Delta}_{\mathcal{T},\mathcal{T}}$, and
augment the anchor estimator with a weighted difference between the target and
source site estimates. In general, since it is possible that source sites may
present empirical bias in estimating the TATE, $\eta_{k}$ should be estimated
in a data-adaptive fashion to safeguard against biased estimates, i.e., to
down weight or exclude likely biased source site estimates of the TATE.

Moreover, the aggregation of the remaining unbiased site-specific
augmentations should also minimize variance. Under the federated learning
setting, the key to evaluate the variance of (8) is to decompose it into
contributions from separate sites so that they can be estimated and shared
within each site. We consider the following decomposition

$\displaystyle\operatorname*{{\rm
Var}}\left(\left\\{1-\sum_{k\in\mathcal{T}^{c}}\eta_{k}\right\\}\widehat{\Delta}_{\mathcal{T},\mathcal{T}}+\sum_{k\in\mathcal{T}^{c}}\eta_{k}\widehat{\Delta}_{\mathcal{T},k}\right)$
$\displaystyle=\operatorname*{{\rm
Var}}\left(\widehat{M}_{\mathcal{T},\mathcal{T}}+\left\\{1-\sum_{k\in\mathcal{T}^{c}}\eta_{k}\right\\}\widehat{\delta}_{\mathcal{T},\mathcal{T}}+\left\\{\sum_{k\in\mathcal{T}^{c}}\eta_{k}\mathbf{d}_{k}^{\top}\right\\}\overline{\boldsymbol{\psi}}_{\mathcal{T}}+\sum_{k\in\mathcal{T}^{c}}\eta_{k}\left\\{\widehat{\delta}_{\mathcal{T},k}-\mathbf{d}_{k}^{\top}\overline{\boldsymbol{\psi}}_{\mathcal{T}}\right\\}\right)$
$\displaystyle\approx\sum_{k\in\mathcal{T}^{c}}\eta_{k}^{2}\operatorname*{{\rm
Var}}\left(\widehat{\delta}_{\mathcal{T},k}\mid\mathscr{D}_{\mathcal{T}}\right)+\sum_{k\in\mathcal{T}}\operatorname*{{\rm
Var}}\left(\frac{n_{k}}{N_{\mathcal{T}}}\left(1,1-\sum_{k\in\mathcal{T}^{c}}\eta_{k},\sum_{k\in\mathcal{T}^{c}}\eta_{k}\mathbf{d}_{k}^{\top}\right)(\widehat{M}_{\mathcal{T},k},\widehat{\delta}_{\mathcal{T},k},\overline{\boldsymbol{\psi}}_{k}^{\top})^{\top}\right),$
(9)

where $\mathbf{d}_{k}$ is the limit for $\widehat{\mathbf{d}}_{k}$, which is
the partial derivative of $\widehat{\delta}_{\mathcal{T},k}$ with respect to
$\overline{\boldsymbol{\psi}}$. We decouple the dependence of the source site
augmentation $\widehat{\delta}_{\mathcal{T},k}$ on the target site(s) through
$\overline{\boldsymbol{\psi}}$ by subtracting the first order approximation of
the dependence
$\mathbf{d}_{k}^{\top}\overline{\boldsymbol{\psi}}_{\mathcal{T}}$. The
resulting
$\widehat{\delta}_{\mathcal{T},k}-\mathbf{d}_{k}^{\top}\overline{\boldsymbol{\psi}}$
is asymptotically independent of the target site(s).

Motivated by (9), we propose the adaptive $L_{1}$ penalized optimal
aggregation

$\displaystyle\widehat{\boldsymbol{\eta}}=$ $\displaystyle\mathop{\rm
arg\min}_{\boldsymbol{\eta}\in\mathbb{R}^{K-|\mathcal{T}|}}N\left[\sum_{k\in\mathcal{T}^{c}}\eta_{k}^{2}\widehat{\sigma}_{k}^{2}+\sum_{k\in\mathcal{T}}\widehat{\mathbf{h}}_{k}(\boldsymbol{\eta})^{\top}\widehat{\Sigma}_{k}\widehat{\mathbf{h}}_{k}(\boldsymbol{\eta})\right]+\lambda\sum_{k\in\mathcal{T}^{c}}|\eta_{k}|\left(\widehat{\Delta}_{\mathcal{T},k}-\widehat{\Delta}_{\mathcal{T},\mathcal{T}}\right)^{2},$
$\displaystyle\widehat{\mathbf{h}}_{k}(\boldsymbol{\eta})=$
$\displaystyle\frac{n_{k}}{n_{\mathcal{T}}}\left(1,1-\sum_{k\in\mathcal{T}^{c}}\eta_{k},\sum_{k\in\mathcal{T}^{c}}\eta_{k}\widehat{\mathbf{d}}_{k}^{\top}\right)^{\top}$
(10)

with $\widehat{\Sigma}_{k}$ from Step 1 and $\widehat{\sigma}_{k}^{2}$ and
$\widehat{\mathbf{d}}_{k}$ from Step 2. The multiplicative $N$ factor is
required to stabilize the loss, as the variance components
$\widehat{\sigma}_{k}^{2}$ and $\widehat{\Sigma}_{k}$ decay to zero at $\sim
1/N$ rate. Choosing $\lambda\asymp N^{\nu}$ with $\nu\in(0,1/2)$, we achieve
an oracle property for the selection and aggregation: (1) the biased site-
specific augmentations have zero weights with large probability; (2) the
regularization on the weights for unbiased site-specific augmentations is
asymptotically negligible ($\ll N^{-1/2}$).

Using the variance estimator (stabilized by “$N$” factor likewise)

$\widehat{\mathcal{V}}=N\sum_{k\in\mathcal{T}^{c}}\widehat{\eta}_{k}\widehat{\sigma}_{k}^{2}+N\sum_{k\in\mathcal{T}}\widehat{\mathbf{h}}_{k}(\widehat{\boldsymbol{\eta}})^{\top}\widehat{\Sigma}_{k}\widehat{\mathbf{h}}_{k}(\widehat{\boldsymbol{\eta}})$
(11)

and the $1-\alpha/2$ quantile for standard normal distribution
$\mathcal{Z}_{\alpha/2}$ we construct the $(1-\alpha)\times 100\%$ confidence
interval

$\widehat{\mathcal{C}}_{\alpha}=\left[\widehat{\Delta}_{\mathcal{T},\sf
FACE}-\sqrt{\widehat{\mathcal{V}}/N}\mathcal{Z}_{\alpha/2},\widehat{\Delta}_{\mathcal{T},\sf
FACE}+\sqrt{\widehat{\mathcal{V}}/N}\mathcal{Z}_{\alpha/2}\right].$ (12)

The full FACE workflow is summarized in Algorithm 1.

Data: For $k=1,...,K$ sites,
$(Y_{ik},\mathbf{X}_{ik}^{\top},A_{i})^{\top},i=1,...,n_{k}$

for  _Target $k\in\mathcal{T}$_ do

2 Calculate $\overline{\boldsymbol{\psi}}_{k}$ and transfer to source sites.
Estimate $\boldsymbol{\alpha}_{k}$, $\boldsymbol{\beta}_{a,k}$ to calculate
the initial TATE $\widehat{\Delta}_{\mathcal{T},k}$ its augmentation
$\widehat{\delta}_{\mathcal{T},k}$ and the variance estimator
$\widehat{\Sigma}_{k}$.

3 end for

4for  _Source sites $k\in\mathcal{T}_{c}$_ do

5 Estimate $\boldsymbol{\gamma}_{k}$, $\boldsymbol{\alpha}_{k}$ and
$\boldsymbol{\beta}_{a,k}$ to calculate the site-specific augmentation
$\widehat{\delta}_{\mathcal{T},k}$ and transfer to target site. Calculate
$\widehat{\sigma}^{2}_{k}$, $\widehat{\mathbf{d}}_{k}$ and transfer to target
site.

6 end for

7for  _Leading Site_ do

8 Estimate $\boldsymbol{\eta}$ by solving the penalized regression in (10).
Construct the final global estimator as $\widehat{\Delta}_{\mathcal{T}}$ by
(8). Calculate the global estimator variance by (11) and construct 95% CI.

9 end for

Result: Global TATE estimate, $\widehat{\Delta}_{\mathcal{T}}$ and $95\%$ CI

Algorithm 1 Pseudocode to obtain TATE estimator leveraging all federated sites

######  Remark 6.

The solution of (10) only involves a small number of summary statistics that
grows proportionally to the dimension of $\boldsymbol{\psi}$ for each site.
Therefore, our aggregation is communication-efficient and privacy-protected.
Classically, procedures like (10) are usually constructed with influence
functions. Let $\widehat{\zeta}_{i}$, $\widehat{\xi}_{i,\mathcal{T}}$ and
$\widehat{\xi}_{i,k}$ be the estimated influence functions for
$\widehat{M}_{\mathcal{T},\mathcal{T}}$,
$\widehat{\delta}_{\mathcal{T},\mathcal{T}}$ and
$\widehat{\delta}_{\mathcal{T},k}$, respectively. The adaptive selection and
aggregation weight $\widehat{\boldsymbol{\eta}}$ can be obtained by solving
the following $L_{1}$ penalized regression

$Q({\boldsymbol{\eta}})=\min_{\boldsymbol{\eta}\in\mathbb{R}^{K-|\mathcal{T}|}}N^{-1}\sum_{i=1}^{N}\left\\{\widehat{\zeta}_{i}+\left(1-\sum_{k\in\mathcal{T}^{c}}\eta_{k}\right)\widehat{\xi}_{i,\mathcal{T}}+\sum_{k\in\mathcal{T}^{c}}\widehat{\xi}_{i,k}\right\\}^{2}+\lambda\sum_{k\in\mathcal{T}^{c}}|\eta_{k}|\left(\widehat{\Delta}_{\mathcal{T},k}-\widehat{\Delta}_{\mathcal{T},\mathcal{T}}\right)^{2}.$

######  Remark 7.

We propose a novel sample splitting approach for the optimal selection of
$\lambda$ without sharing individual-level data. We split the data into
training and validation datasets across all sites. In the training dataset, we
estimate our nuisance parameters $\boldsymbol{\alpha}_{k}$,
$\boldsymbol{\beta}_{a}$, and $\boldsymbol{\gamma}_{k}$ and influence
functions, and solve $Q(\boldsymbol{\eta})$ distributively for a grid of
$\lambda$ values. Using the associated $\boldsymbol{\eta}$ weights from each
value of $\lambda$, we estimate the MSE in the validation data. We set the
value of the optimal tuning parameter, $\lambda_{\mathrm{opt}}$, to be the
value that minimizes the MSE in the validation data.

###  3.4 FACE Under Logistic Regression Models

We showcase FACE with the following realization under logistic regression
models. Let $Y$ be a binary outcome, $K=5$ and $\mathcal{T}=\\{1\\}$. For
notational ease, let $\mathbf{X}$ be the vector of covariates with an
intercept term. We fit logistic regression models with link
$g(x)=1/(1+e^{-x})$ and loss $\ell(y,x)=\log(1+e^{x})-yx$ for all propensity
score and outcome regression models. The mapping $\boldsymbol{\psi}$ density
ratio model is set as the identity map,
$\boldsymbol{\psi}(\mathbf{X})=\mathbf{X}$, which corresponds to a multinomial
logistic regression model on the site indicator $R$.

In Step 1, we first calculate the mean covariate vector in the target site
$k=1$ as
$\overline{\boldsymbol{\psi}}_{\mathcal{T}}=\frac{1}{n_{1}}\sum_{i\in\mathcal{I}_{1}}\mathbf{X}_{i}$
and transfer it to sites 2 through 5. Then, we estimate the models for $k=1$

$\widehat{\boldsymbol{\alpha}}_{1}=\mathop{\rm
arg\min}_{\boldsymbol{\alpha}\in\mathbb{R}^{p+1}}\frac{1}{n_{1}}\sum_{i\in\mathcal{I}_{1}}\ell(A_{i},\boldsymbol{\alpha}^{\top}\mathbf{X}_{i}),\;\widehat{\boldsymbol{\beta}}_{a,1}=\mathop{\rm
arg\min}_{\boldsymbol{\beta}\in\mathbb{R}^{p+1}}\frac{1}{n_{1}}\sum_{i\in\mathcal{I}_{1}}I(A_{i}=a)\ell(Y_{i},\boldsymbol{\alpha}^{\top}\mathbf{X}_{i}).$

Using the estimated models, we obtain the initial estimator and its
augmentation

$\displaystyle\widehat{M}_{\mathcal{T},\mathcal{T}}=$
$\displaystyle\frac{1}{n_{1}}\sum_{i\in\mathcal{I}_{1}}\left\\{g\left(\widehat{\boldsymbol{\beta}}_{1,1}^{\top}\mathbf{X}_{i}\right)-g\left(\widehat{\boldsymbol{\beta}}_{0,1}^{\top}\mathbf{X}_{i}\right)\right\\},$
$\displaystyle\widehat{\delta}_{\mathcal{T},\mathcal{T}}=$
$\displaystyle\frac{1}{n_{1}}\sum_{i\in\mathcal{I}_{1}}\left[\frac{A_{i}}{g\left(\widehat{\boldsymbol{\alpha}}_{1}^{\top}\mathbf{X}_{i}\right)}\left\\{Y_{i}-g\left(\widehat{\boldsymbol{\beta}}_{1,1}^{\top}\mathbf{X}_{i}\right)\right\\}-\frac{1-A_{i}}{g\left(-\widehat{\boldsymbol{\alpha}}_{1}^{\top}\mathbf{X}_{i}\right)}\left\\{Y_{i}-g\left(\widehat{\boldsymbol{\beta}}_{0,1}^{\top}\mathbf{X}_{i}\right)\right\\}\right]$

and
$\widehat{\Delta}_{\mathcal{T},\mathcal{T}}=\widehat{M}_{\mathcal{T},\mathcal{T}}+\widehat{\delta}_{\mathcal{T},\mathcal{T}}$.
The variance covariance matrix estimator $\widehat{\Sigma}_{1}$ can be
calculated as
$\widehat{\Sigma}_{1}=\frac{1}{n_{1}^{2}}\sum_{i\in\mathcal{I}_{1}}\widehat{\mathbf{U}}_{i}\widehat{\mathbf{U}}_{i}^{\top}$
through the estimated influence functions, where
$\widehat{\mathbf{U}}_{i}=(\widehat{\zeta}_{i},\widehat{\xi}_{i},\boldsymbol{\psi}(\mathbf{X}_{i})^{\top})^{\top},$
and the exact form of $\widehat{\xi}_{i,1}$ and $\widehat{\zeta}_{i}$ are
given in the Supplement.

In Step 2, we estimate the models for $k=2,\dots,5$

$\displaystyle\widehat{\boldsymbol{\alpha}}_{k}=$ $\displaystyle\mathop{\rm
arg\min}_{\boldsymbol{\alpha}\in\mathbb{R}^{p+1}}\frac{1}{n_{k}}\sum_{i\in\mathcal{I}_{k}}\ell(A_{i},\boldsymbol{\alpha}^{\top}\mathbf{X}_{i}),\;\widehat{\boldsymbol{\beta}}_{a,k}=\mathop{\rm
arg\min}_{\boldsymbol{\beta}\in\mathbb{R}^{p+1}}\frac{1}{n_{k}}\sum_{i\in\mathcal{I}_{k}}I(A_{i}=a)\ell(Y_{i},\boldsymbol{\alpha}^{\top}\mathbf{X}_{i}).$
$\displaystyle\widehat{\boldsymbol{\gamma}}_{k}=$ $\displaystyle\mathop{\rm
arg\min}_{\boldsymbol{\gamma}\in\mathbb{R}^{p+1}}\frac{1}{n_{k}}\sum_{i\in\mathcal{I}_{k}}\exp(\boldsymbol{\gamma}^{\top}\mathbf{X}_{i})-\boldsymbol{\gamma}^{\top}\overline{\boldsymbol{\psi}}_{\mathcal{T}}.$

Using the estimated models, we obtain the the site-specific augmentations

$\widehat{\delta}_{\mathcal{T},k}=\frac{1}{n_{k}}\sum_{i\in\mathcal{I}_{k}}e^{\widehat{\boldsymbol{\gamma}}_{k}^{\top}\mathbf{X}_{i}}\left[\frac{A_{i}}{g\left(\widehat{\boldsymbol{\alpha}}_{k}^{\top}\mathbf{X}_{i}\right)}\left\\{Y_{i}-g\left(\widehat{\boldsymbol{\beta}}_{1,k}^{\top}\mathbf{X}_{i}\right)\right\\}-\frac{1-A_{i}}{g\left(-\widehat{\boldsymbol{\alpha}}_{k}^{\top}\mathbf{X}_{i}\right)}\left\\{Y_{i}-g\left(\widehat{\boldsymbol{\beta}}_{0,k}^{\top}\mathbf{X}_{i}\right)\right\\}\right].$

The partial derivative of $\widehat{\delta}_{\mathcal{T},k}$ with respect to
$\overline{\boldsymbol{\psi}}_{\mathcal{T}}$ is estimated by

$\displaystyle\widehat{\mathbf{d}}_{k}=$
$\displaystyle-\left(\frac{1}{n_{k}}\sum_{j\in\mathcal{I}_{k}}e^{\widehat{\boldsymbol{\gamma}}_{k}^{\top}\mathbf{X}_{j}}\left[\frac{A_{j}}{g\left(\widehat{\boldsymbol{\alpha}}_{k}^{\top}\mathbf{X}_{j}\right)}\left\\{Y_{j}-g\left(\widehat{\boldsymbol{\beta}}_{1,k}^{\top}\mathbf{X}_{j}\right)\right\\}\right.\right.$
$\displaystyle\hskip
30.00005pt\left.\left.-\frac{1-A_{j}}{g\left(-\widehat{\boldsymbol{\alpha}}_{k}^{\top}\mathbf{X}_{j}\right)}\left\\{Y_{j}-g\left(\widehat{\boldsymbol{\beta}}_{0,k}^{\top}\mathbf{X}_{j}\right)\right\\}\right]\mathbf{X}_{j}^{\top}\right)\left\\{\frac{1}{n_{k}}\sum_{j\in\mathcal{I}_{k}}e^{\widehat{\boldsymbol{\gamma}}_{k}^{\top}\mathbf{X}_{j}}\mathbf{X}_{j}\mathbf{X}_{j}^{\top}\right\\}^{-1}$

The variance estimator $\widehat{\sigma}^{2}_{k}$ can be calculated as
$\widehat{\sigma}^{2}_{k}=\frac{1}{n_{k}}\sum_{i\in\mathcal{I}_{k}}\widehat{\xi}_{i,k}^{2}$
through the estimated influence function, where the form of
$\widehat{\xi}_{i,k}$ is given in the Supplement.

In Step 3, we use $\widehat{\Sigma}_{k}$, $\widehat{\mathbf{d}}_{k}$,
$\widehat{\sigma}^{2}_{k}$, $\widehat{\delta}_{\mathcal{T},k}$ and
$\widehat{\delta}_{\mathcal{T},\mathcal{T}}$ to solve the adaptive selection
and aggregation (10), which leads to the FACE estimator
$\widehat{\Delta}_{\mathcal{T},\sf FACE}$ and the confidence interval
$\widehat{\mathcal{C}}_{\alpha}$.

##  4 Theoretical Guarantees

In this section, we provide the theoretical results for our target only
estimator, source site estimators, and adaptive combination estimator. We
start with a high-level theory for a generic choice of models in Section 4.1
and give in Section 4.2 a detailed set of conditions corresponding to the
realization of Section 3.4. In our asymptotic theory, the total sample size
$N$ may grow but the distribution for a single copy of data
$(Y,\mathbf{X}^{\top},A,R)^{\top}$, including the number of sites, stays the
same.

###  4.1 Theory for General FACE

For the flexible FACE framework, we introduce a set of high-level regularity
conditions.

######  Assumption 2.

For absolute constants $M,\varepsilon>0$,

  1. (a)

(Regularity of TATE estimators) The estimators
$\widehat{M}_{\mathcal{T},\mathcal{T}}$,
$\widehat{\delta}_{\mathcal{T},\mathcal{T}}$ and
$\widehat{\Delta}_{\mathcal{T},\sf FACE}$ admit the following approximations
by influence functions

$\displaystyle\sqrt{N_{\mathcal{T}}}(\widehat{M}_{\mathcal{T},\mathcal{T}}-\overline{M}_{\mathcal{T},\mathcal{T}})=\frac{1}{\sqrt{N_{\mathcal{T}}}}\sum_{k\in\mathcal{T}}\sum_{i\in\mathcal{I}_{k}}\zeta_{i}+o_{p}(1),$
$\displaystyle\sqrt{N_{\mathcal{T}}}(\widehat{\delta}_{\mathcal{T},\mathcal{T}}-\overline{\delta}_{\mathcal{T},\mathcal{T}})=\frac{1}{\sqrt{N_{\mathcal{T}}}}\sum_{k\in\mathcal{T}}\sum_{i\in\mathcal{I}_{k}}\xi_{i,\mathcal{T}}+o_{p}(1),$
$\displaystyle\sqrt{n_{k}}(\widehat{\delta}_{\mathcal{T},k}-\overline{\delta}_{\mathcal{T},k})=\frac{1}{\sqrt{n_{k}}}\sum_{i\in\mathcal{I}_{k}}\xi_{i,k}+\sqrt{n_{k}}\overline{\mathbf{d}}^{\top}\left[\overline{\boldsymbol{\psi}}_{\mathcal{T}}-\mathbb{E}\\{\boldsymbol{\psi}(\mathbf{X})\mid
R\in\mathcal{T}\\}\right]+o_{p}(1),$

with asymptotic limits $\overline{M}_{\mathcal{T},\mathcal{T}}$,
$\overline{\delta}_{\mathcal{T},\mathcal{T}}$,
$\overline{\delta}_{\mathcal{T},k}$, $\overline{\mathbf{d}}_{k}$ and mean zero
random variables $\zeta_{i}$, $\xi_{i,\mathcal{T}}$, $\xi_{i,k}$.

  2. (b)

(Stable variance) The variances of $\zeta_{i}$, $\xi_{i,\mathcal{T}}$,
$\xi_{i,k}$ are in the set $[\varepsilon,M]$. The covariates in the density
ratio are bounded $\|\boldsymbol{\psi}(\mathbf{X})\|_{\infty}\leq M$ almost
surely and have stable variance $\operatorname*{{\rm
Var}}(\boldsymbol{\psi}(\mathbf{X})\mid R\in\mathcal{T})\geq\varepsilon$. The
variance-covariance matrix

$\Sigma_{k}=\operatorname*{{\rm
Var}}\left\\{\left(\zeta_{i},\xi_{i,\mathcal{T}},\boldsymbol{\psi}(\mathbf{X}_{i})^{\top}\right)\right\\}$

has eigenvalues all in $[\varepsilon,M]$.

  3. (c)

(Regularity of auxiliary estimators) The estimators $\widehat{\Sigma}_{k}$,
$\widehat{\sigma}^{2}_{k}$, $\widehat{\mathbf{d}}_{k}$ are
$\sqrt{N}$-consistent

$\displaystyle\sum_{k\in\mathcal{T}}\left\|n_{k}\widehat{\Sigma}_{k}-\operatorname*{{\rm
Var}}\left\\{\left(\zeta_{i},\xi_{i,\mathcal{T}},\boldsymbol{\psi}(\mathbf{X}_{i})^{\top}\right)^{\top}\mid
R_{i}=k\right\\}\right\|$
$\displaystyle+\sum_{k\in\mathcal{T}^{c}}\left\\{\left|n_{k}\widehat{\sigma}^{2}_{k}-\operatorname*{{\rm
Var}}(\xi_{i,k}\mid
R_{i}=k)\right|+\left\|\widehat{\mathbf{d}}_{k}-\overline{\mathbf{d}}_{k}\right\|\right\\}=O_{p}\left(N^{-1/2}\right).$

  4. (d)

(Double robustness) For each target site $k\in\mathcal{T}$, at least one of
the two models is correctly specified:

    1. -i

the propensity score model is correct:
$\mathbb{P}(A=a\mid\mathbf{X},R=k)=\pi_{k}(a,\mathbf{X};\overline{\boldsymbol{\alpha}}_{k})$
for some $\|\overline{\boldsymbol{\alpha}}_{k}\|<\infty$ and all
$a\in\\{0,1\\}$, $k\in\mathcal{T}$.

    2. -ii

the outcome regression model is correct:
$\mathbb{E}(Y\mid\mathbf{X},A=a,R=k)=m_{k}(a,\mathbf{X};\overline{\boldsymbol{\beta}}_{a,k})$
for some $\|\overline{\boldsymbol{\beta}}_{a,k}\|<\infty$ and all
$a\in\\{0,1\\}$, $k\in\mathcal{T}$.

Assumptions 2(a) and 2(d) are the typical regularity conditions under
classical parametric models. They can be verified in two steps: 1) asymptotic
normality of model estimators (Van der Vaart, 2000) and 2) local expansion of
the estimators. Assumption 2(b) regulates the scale of variability of the
data, which leads to a stable variance for $\widehat{\Delta}_{\mathcal{T},\sf
FACE}$. Assumption 2(d) ensures identification of the true TATE by anchoring
on $\widehat{\Delta}_{\mathcal{T},\mathcal{T}}$ (Bang & Robins, 2005). Note
that in the setting of multiple target sites, Assumption 2(d) allows for each
target site to have different correct model specifications for either the OR
model or the PS model.

We now state the theory for the general FACE estimation.

######  Theorem 1.

Under Assumptions 1 and 2, the FACE estimator is consistent and asymptotically
normal with consistent variance estimation $\widehat{\mathcal{V}}$,

$\sqrt{N/\widehat{\mathcal{V}}}\left(\widehat{\Delta}_{\mathcal{T},\sf
FACE}-\Delta_{\mathcal{T}}\right)\leadsto N(0,1).$

We use $\leadsto$ for convergence in distribution.

Base on Theorem 1, we justify that the interval estimation (12) provides
asymptotically honest coverage.

######  Corollary 1.

Under Assumptions 1 and 2, the coverage rate of the confidence interval (12)
approaches the nominal level asymptotically

$\lim_{N\to\infty}\mathbb{P}\left(\Delta_{\mathcal{T}}\in\widehat{\mathcal{C}}_{\alpha}\right)=1-\alpha$

A key step in the proof of Theorem 1 is the analysis of the $L_{1}$ penalized
adaptive selection and aggregation (10). The problem is different from the
typical penalized regression, so we develop a new proof strategy. We first
analyze the optimal combination with oracle selection, in which the biased
augmentations are excluded. For unbiased augmentations,
$\widehat{\Delta}_{\mathcal{T},k}-\widehat{\Delta}_{\mathcal{T},\mathcal{T}}=O_{p}\left(N^{-1/2}\right)$,
so the penalty term is asymptotically negligible
$\lambda(\widehat{\Delta}_{\mathcal{T},k}-\widehat{\Delta}_{\mathcal{T},\mathcal{T}})^{2}=o_{p}\left(N^{-1/2}\right)$.
Thus, the estimated combination converges to the asymptotic limit at the
regular $N^{-1/2}$ rate. Finally, we show that the estimated combination with
oracle selection also solves the original problem with large probability.

######  Remark 8.

For consistency of $\widehat{\Delta}_{\mathcal{T},\sf FACE}$, we require that
the propensity score or outcome regression model is correct for the target
site so that the initial doubly robust TATE estimator
$\widehat{\Delta}_{\mathcal{T},\mathcal{T}}$ is consistent. Even if the models
for source site and density ratio are misspecified, we can safely aggregate
the site-specific augmentations by the adaptive selection in (10). To
meaningfully leverage information from source sites for TATE, we would expect
that many $k$ among source sites

  1. 1.

satisfy the ignorability condition 1(d);

  2. 2.

share the same outcome regression model with the target sites

$\mathbb{E}\\{Y\mid A=a,\mathbf{X},R=k\\}=\mathbb{E}\\{Y\mid
A=a,\mathbf{X},R\in\mathcal{T}\\};$

  3. 3.

either the outcome regression model $m_{a,k}$ is correct, or both the
propensity score $\pi_{k}$ and the density ratio $\omega_{k}$ models are
correct.

For source sites satisfying the conditions above, their site-specific
augmentations are unbiased and thus contribute to the efficiency improvement
of $\widehat{\Delta}_{\mathcal{T},\sf FACE}$.

###  4.2 Application for the Realization

For the realization in Section 3.4, we denote the asymptotic parameters as

$\displaystyle\overline{\boldsymbol{\alpha}}_{k}=$ $\displaystyle\mathop{\rm
arg\min}_{\boldsymbol{\alpha}\in\mathbb{R}^{p}}\mathbb{E}\\{\ell(A,\boldsymbol{\alpha}^{\top}\mathbf{X})\mid
R=k\\},$ $\displaystyle\overline{\boldsymbol{\beta}}_{a,k}=$
$\displaystyle\mathop{\rm
arg\min}_{\boldsymbol{\alpha}\in\mathbb{R}^{p}}\mathbb{E}\\{\ell(Y,\boldsymbol{\alpha}^{\top}\mathbf{X})\mid
A=a,R=k\\},$ $\displaystyle\overline{\boldsymbol{\gamma}}_{k}=$
$\displaystyle\mathop{\rm
arg\min}_{\boldsymbol{\gamma}\in\mathbb{R}^{q}}\mathbb{E}\\{\exp(\boldsymbol{\gamma}^{\top}\mathbf{X})-\boldsymbol{\gamma}^{\top}\mathbb{E}(\mathbf{X}\mid
R\in\mathcal{T})\mid R=k\\}.$

We give a mild set of sufficient conditions for Assumption 2.

######  Assumption 3.

For absolute constants $M,\varepsilon>0$,

  1. (a)

(Design) $\|\mathbf{X}\|_{\infty}\leq M$ almost surely, and all eigenvalues of
$\mathbb{E}(\mathbf{X}\mathbf{X}^{\top})$ are in $[\varepsilon,M]$.

  2. (b)

(Overlap) For all $k=1,\dots,K$, $a=0,1$ and $i\in\mathcal{I}_{k}$,
$g(\overline{\boldsymbol{\alpha}}_{k}^{\top}\mathbf{X}_{i})$,
$g^{\prime}(\overline{\boldsymbol{\beta}}_{a,k}^{\top}\mathbf{X}_{i})$ and
$\exp\\{\overline{\boldsymbol{\gamma}}_{k}^{\top}\mathbf{X}_{i}\\}$ are in
$[\varepsilon,1-\varepsilon]$ almost surely.

  3. (c)

(Double robustness) For each target site $k\in\mathcal{T}$, at least one of
the two models is correctly specified:

    1. -i

the propensity score model is correct:
$\mathbb{P}(A=1\mid\mathbf{X},R=k)=g(\overline{\boldsymbol{\alpha}}_{k}^{\top}\mathbf{X});$

    2. -ii

the outcome regression model is correct:
$\mathbb{E}(Y\mid\mathbf{X},A=a,R=k)=g(\overline{\boldsymbol{\beta}}_{a,k}^{\top}\mathbf{X}).$

After verifying that Assumptions 1 and 3 imply the generic Assumption 2, we
can apply Theorem 1 in that realization.

######  Corollary 2.

Under the setting of Section 3.4 and Assumptions 1 and 3, the FACE estimator
is consistent and asymptotically normal with consistent variance estimation
$\widehat{\mathcal{V}}$,

$\sqrt{N/\widehat{\mathcal{V}}}\left(\widehat{\Delta}_{\mathcal{T},\sf
FACE}-\Delta_{\mathcal{T}}\right)\leadsto N(0,1).$

##  5 Simulation Studies

We study the finite sample performance of our proposed data-adaptive
estimators and make comparisons with an estimator that leverages target data
only and a sample-size adjusted estimator that does not adaptively weight
different sites. In the simulation studies, we take the target to be a single
site. We examine the empirical bias, RMSE, and coverage of the $95\%$ CI over
$1000$ simulations. We vary the number of sites $K\in\\{10,20,50\\}$ and the
true outcome regression, propensity score, and density ratio models.

The covariates in each site $\mathbf{X}_{kp}$ are generated from a skewed
normal distribution,
$\mathbf{X}_{kp}\sim\mathcal{SN}(\mathbf{x};\kappa_{kp},\sigma^{2}_{kp},\nu_{kp})$,
where $k=1,...,K$ indexes the sites and $p=1,...,P$ indexes the covariates,
$\kappa_{kp}$ is the location parameter, $\sigma_{kp}$ is the scale parameter,
and $\nu_{kp}$ is the skewness parameter. For all sites, we let
$\kappa_{k\cdot}\in(0.10,0.15)$ and $\sigma_{k\cdot}=(1,...,1)$. For the
target site, we set $\nu_{k\cdot}=0$. For the source sites, we let
$\nu_{k\cdot}\in\\{-4/p,4/p\\}$. We generate the sample size in source sites
to be $n_{k}\sim\mathrm{Unif}(100,300)$.

The true potential outcomes are generated as
$Y_{k}(1)=(\mathbf{X}_{k}-\mu_{1})\boldsymbol{\beta}_{11}+\mathbf{X}_{k}^{\circ
2}\boldsymbol{\beta}_{21}+3+\varepsilon_{k}$ and
$Y_{k}(0)=(\mathbf{X}_{k}-\mu_{1})\boldsymbol{\beta}_{10}+\mathbf{X}_{k}^{\circ
2}\boldsymbol{\beta}_{20}+\varepsilon_{k},$ where
$\varepsilon_{k}\sim\mathcal{N}(0,2.25P^{2})$,
$\boldsymbol{\beta}_{11}=2\cdot(0.4,..,1.2)/P$, and
$\boldsymbol{\beta}_{10}=(0.4,..,1.2)/P$. The true propensity score model is
generated as $A_{k}\mid\mathbf{X}=\mathbf{x}\sim\text{Bernoulli}(\pi_{k})$,
$\pi_{k}=\text{expit}(\mathbf{X}_{k}\boldsymbol{\alpha}_{1k}+\mathbf{X}_{k}^{\circ
2}\boldsymbol{\alpha}_{2k}),$ where for the target site,
$\boldsymbol{\alpha}_{11}=(0.4,...,-0.4)$, $\alpha_{21}=0$. For the source
sites $k$, $\boldsymbol{\alpha}_{1k}=(0.5,...,-0.5)$ and
$\boldsymbol{\alpha}_{2k}=0$ .

In Setting I, we examine the scenario where both the outcome regression and
propensity score models are correctly specified by setting
$\boldsymbol{\beta}_{21}=0$ and $\boldsymbol{\beta}_{20}=0$ and only adjusting
for linear terms. In Setting II, we instead have
$\boldsymbol{\beta}_{21}=(0.2,..,0.4)$, so that the true $Y_{k}(1)$ and
$Y_{k}(0)$ include quadratic terms. We misspecify the outcome regression by
fitting linear models. In Setting III, compared to Setting I, we have
$\boldsymbol{\alpha}_{21}=(0.12,...,-0.12)$. We misspecify the propensity
score model by fitting a linear model. In Setting IV, the outcome model and
the propensity score model are misspecified for half of the source sites.

In each setting, we fix the number of covariates to be $P=10$, choose the
tuning parameter
$\lambda\in\\{0,10^{-4},10^{-3},10^{-2},0.1,0.25,0.5,1,2,5,10\\}$ by
distributed cross validation, where we split the simulated data set in each
site into two equally sized samples and proceed to solve for the
$\boldsymbol{\eta}$ that minimizes the MSE of the $Q(\boldsymbol{\eta})$
function in the validation data set.

In Table 1, we summarize the bias, root mean squared error (RMSE), and
coverage probability of the $95\%$ CI of a target-only estimator (Target), a
sample-size weighted estimator that assumes homogeneity in the covariate
distributions across sites (SS), and our adaptive estimator (FACE) over $1000$
simulations across Settings I-IV. The results show that FACE performs well in
all settings, with minimal bias, substantially reduced RMSE compared to the
target-only estimator and the sample-size weighted estimator, and nominal
coverage. Further, FACE shows decreasing RMSE as $K$ increases.

Table 1: Bias, RMSE, and Coverage probability of the $95\%$ CI of estimators over $1000$ simulations. | Simulation scenarios
---|---
| $K=10$ | $K=20$ | $K=50$
| Bias | RMSE | Cov. | Bias | RMSE | Cov. | Bias | RMSE | Cov.
Setting I |  |  |  |  |  |  |  |  |
Target | 0.09 | 1.67 | 97.00 | 0.09 | 1.67 | 97.00 | 0.09 | 1.67 | 97.00
SS | 0.22 | 0.46 | 99.50 | 0.19 | 0.33 | 99.00 | 0.18 | 0.25 | 94.90
FACE | 0.06 | 0.53 | 99.00 | 0.08 | 0.37 | 100.00 | 0.07 | 0.23 | 98.00
Setting II |  |  |  |  |  |  |  |  |
Target | 0.02 | 3.38 | 93.90 | 0.02 | 3.38 | 93.90 | 0.11 | 2.48 | 96.00
SS | 0.22 | 0.61 | 99.90 | 0.19 | 0.45 | 99.50 | 0.19 | 0.31 | 98.60
FACE | 0.11 | 0.75 | 100.00 | 0.11 | 0.53 | 99.70 | 0.10 | 0.34 | 99.00
Setting III |  |  |  |  |  |  |  |  |
Target | 0.10 | 1.92 | 94.00 | 0.10 | 1.92 | 94.00 | 0.14 | 1.95 | 93.00
SS | 0.30 | 0.96 | 100.00 | 0.21 | 0.64 | 100.00 | 0.14 | 0.42 | 100.00
FACE | 0.02 | 0.54 | 98.00 | 0.02 | 0.38 | 100.00 | 0.03 | 0.22 | 99.00
Setting IV |  |  |  |  |  |  |  |  |
Target | 0.08 | 2.48 | 96.00 | 0.08 | 2.48 | 96.00 | 0.11 | 2.48 | 96.00
SS | 0.48 | 1.67 | 100.00 | 0.36 | 0.96 | 100.00 | 0.33 | 0.64 | 100.00
FACE | 0.12 | 0.84 | 100.00 | 0.08 | 0.54 | 100.00 | 0.10 | 0.37 | 98.00

##  6 Comparative Effectiveness of COVID-19 Vaccines

To illustrate FACE, we study the comparative effectiveness of BNT162b2
(Pfizer) versus mRNA-1273 (Moderna) for the prevention of COVID-19 outcomes in
the VA healthcare system. It is of interest to understand the real-world
effectiveness of these vaccines, but head-to-head comparisons have been rare.
A recent emulated target trial using the EHRs of US veterans showed that the
24-week risk of COVID-19 outcomes was low for patients who received either
vaccine, but lower for veterans assigned to Moderna compared to Pfizer
(Dickerman et al., n.d.). Utilizing FACE, we examine the TATE in a federated
data setting where the target population of interest is one of five sites
(North Atlantic, Southwest, Midwest, Continental, or Pacific) in the VA
healthcare system. Our problem is more challenging than that of Dickerman et
al. (n.d.) due to the federated data setting and the different target
populations of interest that we are able to study.

Inclusion criteria included veteran status, at least 18 years of age by
January 1, 2021, no previously documented COVID-19 infection, no previous
COVID-19 vaccination, and documented two-dose COVID-19 vaccination with either
Pfizer or Moderna between January 1 and March 24, 2021. For each eligible
veteran, follow-up began on the day that the second dose of vaccine was
received (baseline) and ended on the day of death, 120 or 180 days after
baseline, or the end of the study time period (September 24, 2021). The
outcomes of interest were documented SARS-CoV-2 infection either 120 or 180
days after baseline and death with COVID-19 infection either 120 or 180 days
after baseline.

Among the $30437$ eligible veterans, $14775$ received Pfizer and $15662$
received Moderna. Baseline characteristics among the two groups were similar
within site, with some heterogeneity across sites, e.g., a larger Hispanic
population in the Southwest and Pacific (Table 2).

Table 2: Baseline characteristics of veterans assigned to Pfizer or Moderna in
each of the five VA sites

The raw event rates for documented COVID-19 infection within 180 days of
receiving the second dose for Pfizer (Moderna) in the five sites were 2.81%
(1.93%) in the North Atlantic, 3.58% (3.23%) in the Southwest, 2.25% (2.08%)
in the Midwest, 2.97% (2.36%) in the Continental, and 2.80% (1.43%) in the
Pacific. The raw event rates for death with COVID-19 infection within 180 days
of receiving the second dose for Pfizer (Moderna) were 0.37% (0.06%) in the
North Atlantic, 0.36% (0.23%) in the Southwest, 0.18% (0.21%) in the Midwest,
0.21% (0.26%) in the Continental, and 0.11% (0.09%) in the Pacific.

Figure 2 shows the TATE estimates for (a) 120-day COVID-19 infection, (b)
180-day COVID-19 infection, (c) 120-day death with COVID-19 infection, and (d)
180-day death with COVID-19 infection. For each outcome, the target population
is taken to be one of the five sites, and three estimators are plotted along
with their 95% confidence interval: the doubly robust estimator that only uses
target site data (Target Only), the sample-size weighted estimator that
leverages each site where $\eta_{k}$ is taken to be $n_{k}/N$ (SS),
$k=1,...,5$, and the FACE estimator. Our results indicate that the FACE
estimator tracks the Target Only estimator closely with small bias compared to
the SS estimator and substantially tighter confidence intervals compared to
the Target Only estimator, resulting in qualitatively different conclusions in
certain cases, e.g., 180-day COVID-19 infection in the Continental, 120-day
death with COVID-19 infection in the Southwest, and 180-day death with
COVID-19 infection in the Midwest, North Atlantic, and Southwest. Using FACE,
our results show that veterans who received Moderna had an approximately 1%
lower rate of 180-day COVID-19 infection, and this difference appeared
consistent across sites.

(a)

(b)

(c)

(d)

Figure 2: Comparative effectiveness of Moderna vs. Pfizer for (a) COVID-19
positive test 120 days after vaccine dose 2, (b) COVID-19 positive test 180
days after vaccine dose 2, (c) death with COVID-19 120 days after vaccine dose
2, and (d) death with COVID-19 180 days after vaccine dose 2

##  7 Discussion

In this paper, we have developed FACE to leverage heterogeneous data from
multiple study sites to more precisely estimate treatment effects for a target
population of interest. FACE improves upon the precision of the target-
population only estimator by leveraging source population information without
inducing bias. FACE is privacy-preserving and communication-efficient,
requiring only one round of communication between sites. Notably, FACE is able
to solve for site-level weights without requiring patient-level information by
using aggregated summary statistics from the influence functions. Furthermore,
in the selection of the optimal tuning parameter for the penalized regression,
we have developed a distributed way of tuning and cross-validation, which to
the best of our knowledge, is a novel contribution to the literature. In
addition to providing theoretical double robustness and efficiency guarantees,
FACE does not rely on prior knowledge of model stability or correct model
specification, which is a substantial improvement on current federated methods
for causal inference (Xiong et al., 2021). We also obtained promising results
from a real world analysis of COVID-19 outcomes for veterans assigned to
either Pfizer or Moderna vaccines among five federated VA sites.

FACE can easily be generalized to the setting where some sites have RCT data.
In such a setting, one could define the target population as the set of trial
participants. When the RCT data is treated as the anchoring site, the target
site propensity score model is known, so the target site estimator for the
TATE is consistent, and the global adaptive estimator is likely to be more
reliable. Our FACE framework can thus be viewed as a contribution to recent
work on using observational studies to reduce the variance associated with
treatment effect estimates from experimental studies (Rosenman et al., 2018;
Athey et al., 2020). For greater generalizability, participants for whom there
is only observational data can be taken to be the target population. FACE can
also be adapted to target different causal parameters of interest, such as the
average treatment effect of the treated (ATT).

When the target population includes multiple sites, increased robustness can
be obtained at the cost of estimating additional density ratio models. We
detail this alternative method for the target-only estimator in the
Supplementary Materials. Future work may consider focusing on estimands
defined by regions of substantial overlap or alternative methods for estimands
defined by subpopulations of interest. For example, the conditional average
treatment effect (CATE) is an important estimand of real-world interest,
particularly for understanding benefits and dangers of treatments for
underrepresented groups and fairness research.

SUPPLEMENTARY MATERIAL

Appendix:

Proofs of theoretical guarantees and derivation of influence functions.

R package for FACE:

FACE contains code to replicate simulation results.

## References

  * (1)
  * Athey et al. (2020) Athey, S., Chetty, R. & Imbens, G. (2020), ‘Combining experimental and observational data to estimate treatment effects on long term outcomes’, arXiv preprint arXiv:2006.09676 .
  * Bang & Robins (2005) Bang, H. & Robins, J. M. (2005), ‘Doubly robust estimation in missing data and causal inference models’, Biometrics 61(4), 962–973.
  * Brat et al. (2020) Brat, G. A., Weber, G. M., Gehlenborg, N., Avillach, P., Palmer, N. P., Chiovato, L., Cimino, J., Waitman, L. R., Omenn, G. S., Malovini, A. et al. (2020), ‘International electronic health record-derived covid-19 clinical course profiles: the 4ce consortium’, medRxiv .
  * Chen & Xie (2014) Chen, X. & Xie, M.-g. (2014), ‘A split-and-conquer approach for analysis of extraordinarily large data’, Statistica Sinica pp. 1655–1684.
  * Chen et al. (2006) Chen, Y., Dong, G., Han, J., Pei, J., Wah, B. W. & Wang, J. (2006), ‘Regression cubes with lossless compression and aggregation’, IEEE Transactions on Knowledge and Data Engineering 18(12), 1585–1599.
  * Dahabreh et al. (2020) Dahabreh, I. J., Robertson, S. E., Steingrimsson, J. A., Stuart, E. A. & Hernan, M. A. (2020), ‘Extending inferences from a randomized trial to a new target population’, Statistics in medicine 39(14), 1999–2014.
  * Degtiar & Rose (2021) Degtiar, I. & Rose, S. (2021), ‘A review of generalizability and transportability’, arXiv preprint arXiv:2102.11904 .
  * Dickerman et al. (n.d.) Dickerman, B. A., Gerlovin, H., Madenci, A. L., Kurgansky, K. E., Ferolito, B. R., Figueroa Muñiz, M. J., Gagnon, D. R., Gaziano, J. M., Cho, K., Casas, J. P. & Hernán, M. A. (n.d.), ‘Comparative effectiveness of bnt162b2 and mrna-1273 vaccines in u.s. veterans’, New England Journal of Medicine .
  * Dong et al. (2020) Dong, L., Yang, S., Wang, X., Zeng, D. & Cai, J. (2020), ‘Integrative analysis of randomized clinical trials with real world evidence studies’, arXiv preprint arXiv:2003.01242 .
  * Duan et al. (2019) Duan, R., Boland, M. R., Liu, Z., Liu, Y., Chang, H. H., Xu, H., Chu, H., Schmid, C. H., Forrest, C. B., Holmes, J. H., Schuemie, M. J., Berlin, J. A., Moore, J. H. & Chen, Y. (2019), ‘Learning from electronic health records across multiple sites: A communication-efficient and privacy-preserving distributed algorithm’, Journal of the American Medical Informatics Association 27(3), 376–385.
  * Duan et al. (2020b) Duan, R., Ning, Y., Wang, S., Lindsay, B., Carroll, R. & Chen, Y. (2020b), ‘A fast score test for generalized mixture models’, Biometrics 76, 811–820.
  * Hernán & Robins (2020) Hernán, M. A. & Robins, J. M. (2020), ‘Causal inference: What if?’.
  * Hong et al. (2017) Hong, C., Chen, Y., Ning, Y., Wang, S., Wu, H. & Carroll, R. (2017), ‘Plemt: A novel pseudolikelihood based EM test for homogeneity in generalized exponential tilt mixture models.’, Journal of the American Statistical Association .
  * Hripcsak et al. (2016) Hripcsak, G., Ryan, P. B., Duke, J. D., Shah, N. H., Park, R. W., Huser, V., Suchard, M. A., Schuemie, M. J., DeFalco, F. J., Perotte, A. et al. (2016), ‘Characterizing treatment pathways at scale using the ohdsi network’, Proceedings of the National Academy of Sciences 113(27), 7329–7336.
  * Imbens & Rubin (2015) Imbens, G. W. & Rubin, D. B. (2015), Causal inference in statistics, social, and biomedical sciences, Cambridge University Press.
  * Josey et al. (2020) Josey, K. P., Yang, F., Ghosh, D. & Raghavan, S. (2020), ‘A calibration approach to transportability with observational data’, arXiv preprint arXiv:2008.06615 .
  * Lee et al. (2017) Lee, J. D., Liu, Q., Sun, Y. & Taylor, J. E. (2017), ‘Communication-efficient sparse regression’, The Journal of Machine Learning Research 18(1), 115–144.
  * Li et al. (2013) Li, R., Lin, D. K. & Li, B. (2013), ‘Statistical inference in massive data sets’, Applied Stochastic Models in Business and Industry 29(5), 399–409.
  * Lian & Fan (2017) Lian, H. & Fan, Z. (2017), ‘Divide-and-conquer for debiased l 1-norm support vector machine in ultra-high dimensions’, The Journal of Machine Learning Research 18(1), 6691–6716.
  * Neyman (1923) Neyman, J. (1923), ‘On the application of probability theory to agricultural experiments’, Statistical Science 5(5), 463–480.
  * Pan & Yang (2009) Pan, S. J. & Yang, Q. (2009), ‘A survey on transfer learning’, IEEE Transactions on knowledge and data engineering 22(10), 1345–1359.
  * Qin (1998) Qin, J. (1998), ‘Inferences for case-control and semiparametric two-sample density ratio models’, Biometrika 85(3), 619–630.
  * Qin & Liang (2011) Qin, J. & Liang, K.-Y. (2011), ‘Hypothesis testing in a mixture case-control model’, Biometrics 67, 182–193.
  * Rosenbaum & Rubin (1983) Rosenbaum, P. R. & Rubin, D. B. (1983), ‘The central role of the propensity score in observational studies for causal effects’, Biometrika 70(1), 41–55.
  * Rosenman et al. (2018) Rosenman, E., Owen, A. B., Baiocchi, M. & Banack, H. (2018), ‘Propensity score methods for merging observational and experimental datasets’, arXiv preprint arXiv:1804.07863 .
  * Rubin (1974) Rubin, D. B. (1974), ‘Estimating causal effects of treatments in randomized and nonrandomized studies.’, Journal of educational Psychology 66(5), 688.
  * Stuart et al. (2018) Stuart, E. A., Ackerman, B. & Westreich, D. (2018), ‘Generalizability of randomized trial results to target populations: design and analysis possibilities’, Research on social work practice 28(5), 532–537.
  * Stuart et al. (2015) Stuart, E. A., Bradshaw, C. P. & Leaf, P. J. (2015), ‘Assessing the generalizability of randomized trial results to target populations’, Prevention Science 16(3), 475–485.
  * Stuart et al. (2011) Stuart, E. A., Cole, S. R., Bradshaw, C. P. & Leaf, P. J. (2011), ‘The use of propensity scores to assess the generalizability of results from randomized trials’, Journal of the Royal Statistical Society: Series A (Statistics in Society) 174(2), 369–386.
  * Van der Vaart (2000) Van der Vaart, A. W. (2000), Asymptotic statistics, Vol. 3, Cambridge university press.
  * Vo et al. (2021) Vo, T. V., Hoang, T. N., Lee, Y. & Leong, T.-Y. (2021), ‘Federated estimation of causal effects from observational data’, arXiv preprint arXiv:2106.00456 .
  * Wang et al. (2019) Wang, X., Yang, Z., Chen, X. & Liu, W. (2019), ‘Distributed inference for linear support vector machine’, Journal of Machine Learning Research 20(113), 1–41.
  * Xiong et al. (2021) Xiong, R., Koenecke, A., Powell, M., Shen, Z., Vogelstein, J. T. & Athey, S. (2021), ‘Federated causal inference in heterogeneous observational data’, arXiv preprint arXiv:2107.11732 .