# Data fusion for efficiency gain in ATE estimation:
A practical review with simulations

Xi Lin  Department of Statistics, University of Oxford, UK  Jens Magelund Tarp
Novo Nordisk, Denmark  Robin J. Evans  Department of Statistics, University of
Oxford, UK

###### Abstract

The integration of real-world data (RWD) and randomized controlled trials
(RCT) is increasingly important for advancing causal inference in scientific
research. This combination holds great promise for enhancing the efficiency of
causal effect estimation, offering benefits such as reduced trial participant
numbers and expedited drug access for patients. Despite the availability of
numerous data fusion methods, selecting the most appropriate one for a
specific research question remains challenging. This paper systematically
reviews and compares these methods regarding their assumptions, limitations,
and implementation complexities. Through simulations reflecting real-world
scenarios, we identify a prevalent risk-reward trade-off across different
methods. We investigate and interpret this trade-off, providing key insights
into the strengths and weaknesses of various methods; thereby helping
researchers navigate through the application of data fusion for improved
causal inference.

Keywords: Causal inference, clinical trials, data fusion, efficiency gain,
external controls

##  1 Introduction

Randomized Controlled Trials (RCT) are widely considered to be the gold
standard for establishing causality, providing high-quality evidence to inform
high-stakes decisions in scientific, medical, econometric, and other research
areas. However, RCTs face significant limitations; for example, in the context
of rare diseases there may be very few sufferers or events, complicating
patient recruitment, limiting statistical power, and raising data protection
issues. These challenges highlight the need for a complementary approach.

Observational data, or real-world data (RWD) are increasingly available in
abundance, as they are routinely collected without designed intervention.
Common examples include electronic health records (eHRs), government
administrative datasets and online user behaviour data. RWD offers advantages
of accessibility, large volumes, and the availability of long-term outcomes.
However, the lack of randomization in treatment allocation can introduce bias
due to confounding variables that are not measured.

The complementary features of RCTs and RWD have motivated the development of
data fusion techniques. Bareinboim and Pearl (2016) define data fusion as the
integration of multiple datasets collected under heterogeneous conditions.
Recognizing the value of real-world evidence, regulatory bodies such as the
U.S Food and Drug Administrations (FDA) and European Medicines Agency (EMA)
support its use, in particular for examining treatment efficacy for rare
diseases. For example, the FDA approved drugs for Fabry disease based on both
clinical trials and RWD due to the rarity of the disease (U.S. Food and Drug
Administration, 2018). Similarly, the EMA granted accelerated approval for a
drug for metastatic Merkel cell carcinoma (mMCC) using evidence from a single-
arm, open-label study compared with historical controls from eHR (Merck,
2017).

Efficiency gains are a primary goal of data fusion. Augmenting RCTs with RWD
can enhance the statistical power of treatment effect estimates for both the
overall population and specific subgroups. While we utilize the richness of
the RWD, a central problem is confounding bias due to the absence of
randomization. To address this, some methods integrate only the control arm of
RWD, known as external or historical controls, to mitigate bias. Common
methods include test-then-pool (Viele et al., 2014), Bayesian dynamic
borrowing (Ibrahim and Chen, 2000; Hobbs et al., 2012; Schmidli et al., 2014),
and more recently, prognostic score adjustment (Schuler et al., 2021). When
treatment data are available in RWD, including both arms can offer greater
efficiency gains. The challenge lies in balancing variance reduction with the
potential increase in confounding bias. Various methodologies seek to optimize
this trade-off from different angles, including James-Stein style shrinkage
estimators (Rosenman et al., 2023a; Green et al., 2005), estimation under the
general semi-parametric framework (Yang et al., 2020; Li and Luedtke, 2023;
Morucci et al., 2023), weighted average of treatment effect estimates from RCT
and RWD (Oberst et al., 2022; Chen et al., 2021), bias-correction (Kallus et
al., 2018; Yang and Ding, 2020; Yang et al., 2020; Cheng and Cai, 2021; van
der Laan et al., 2024), and Bayesian dynamic borrowing through a power
likelihood (Lin et al., 2024). Section 3 provides a detailed review of these
methods.

Apart from efficiency gains, data fusion can also address generalizability,
transportability, and long-term impact studies. A limitation of RCTs is that
they might not be representative of the target population due to eligibility
criteria, which affects the validity of results. Data fusion addresses this by
reweighting subgroups or extrapolating according to a target population
represented within the RWD. Additionally, RWD can complement RCTs for long-
term impact studies, extending the time horizon and predicting outcomes. For
example, Athey et al. (2020) combined an RCT with administrative education
data and concluded long-term benefits of smaller class sizes, and Robbins et
al. (2022) fused a trial with national mortality registry data and showed that
subsidized health insurance improves long-term mortality. This review,
however, focuses primarily on the efficiency gains achieved through data
fusion.

With numerous new data fusion methods being proposed in recent literature, it
is important to understand which methods are suitable in which settings and
how they compare. Although several surveys, such as Colnet et al. (2020); Shi
et al. (2023); Brantner et al. (2023), have reviewed these methods, their
discussions predominantly remain theoretical. These methods have not been
thoroughly compared to one another in simulation studies reflecting realistic
scenarios. Additionally, these methods have not been tested out extensively in
real data analysis.

This paper aims to bridge this critical knowledge gap by presenting
comprehensive simulations to evaluate various methods in terms of
practicality, sensitivity to assumptions, and performance. Our study focuses
on using individual-level RWD to augment an RCT, aiming for more efficient
estimation of the average treatment effect (ATE) in a target population. This
represents a common use case in medical, economic, and public policy research.
In a subsequent paper, building on learnings from this review, we will provide
a detailed practical example, demonstrating the application of these fusion
methods in a real-world context. Specifically, we will augment the PIONEER 6
medical trial with U.S. health insurance claims records to assess
semaglutide’s impact on cardiovascular outcomes.

The remainder of this paper is organized as follows: Section 2 introduces the
notations, causal estimands and common identification assumptions. Section 3
reviews the fusion methods designed to achieve efficiency gains. Section 4
presents the results of our simulations, offering an in-depth analysis and
interpretation of each method’s performance.. Finally, Section 5 concludes
with a discussion and offers concrete advice for practitioners, guiding them
in selecting the most suitable method for their specific purposes.

##  2 Preliminaries

In this section, we introduce the notations used in this paper, various causal
estimands, and assumptions that are commonly made for data fusion
methodologies.

###  2.1 Notations

We assume two data sources: an RCT and some observational data. Let $A$ be the
treatment, $\bm{X}$ be a vector of pre-treatment covariates, $Y$ be the
outcome of interest, and $U$ represent unmeasured confounding. Let $S$ be a
binary indicator of data source (1: RCT, 0: observational). Random variables
are denoted by capital letters (e.g., $\bm{X}$, $Y$) and their instantiations
by corresponding lowercase letters (e.g., $\bm{X}$, $Y$). Bold letters (e.g.,
$\bm{X}$) indicate vectors, as opposed to scalars.

The RCT data consist of i.i.d. observations
$\left\\{V_{i}=\left(Y_{i},A_{i},\bm{X}_{i}\right),S_{i}=1):i\in\mathcal{R}\right\\}$
with sample size $n_{r}$, and the observational data consist of i.i.d.
observations
$\left\\{V_{i}=\left(Y_{i},A_{i},\bm{X}_{i}\right),S_{i}=0):i\in\mathcal{O}\right\\}$
with sample size $n_{o}$, where $\mathcal{R}$ and $\mathcal{O}$ are sample
index sets for the two data sources. Each RCT observation can be modelled as
sampled from a distribution $P_{R}$, and similarly $P_{O}$ for observational
data. Typically, we would expect $n_{r}<n_{o}$ to represent a common scenario
in practice where researchers want to augment a small-scale RCT with a larger
observational dataset, yet this is not required.

$X$$A$$Y$$U$(a)$X$$A$$Y$$U$(b) Figure 1: Causal models for (a) RCT data and
(b) the observational data. Dashed edges denote possible causal relationships.
For example, In (a), if edge $X\to A$ exists, it means that the treatment is
conditionally randomized, otherwise completely randomized.

There are several overlapping frameworks to represent causal relationships,
including potential outcomes (Rubin, 1974), structural causal models (e.g.
Pearl, 2009) based on causal directed graphs (e.g. Spirtes et al., 2000),
single world intervention graphs (Richardson and Robins, 2013), and the
decision-theoretic framework (Dawid, 2021). In this review, we follow the
potential outcome framework to define causal effects. Let $Y(a)$ be the
potential outcome had the subject been assigned treatment $A=a$.

###  2.2 Causal estimands

There exists various causal quantities that researches may be interested in
depending on the specific research questions. In this subsection, we introduce
several of the most common causal estimands.

######  Definition 2.1.

(Average treatment effect) $\tau=\mathbb{E}[Y(1)-Y(0)].$

This is also known as the population average treatment effect (PATE) because
it is an expectation taken across a certain target population. We can further
specify exactly the population to pick out, for example, the trial population
$P_{R}$ from which the RCT samples are drawn, the observational population
$P_{O}$, or some other target population denoted by $P_{T}$.

If we further calculate the PATE conditional on subjects receiving treatment,
then we derive the estimands of population average effect for the treated
(PATT).

######  Definition 2.2.

(Population average effect for the treated)
$\tau^{\text{T}}=\mathbb{E}[Y(1)-Y(0)\mid A=1].$

We can define the population average effect for the controls (PATC) similarly.

An alternative causal parameter is the sample average treatment effect (SATE),
which is defined as the average difference in outcomes between the treated and
untreated groups within a specific sample drawn from the population, rather
than the population at large (Imbens, 2004).

######  Definition 2.3.

(Sample average treatment effect)

$\tau_{S}=\frac{1}{N}\sum_{i=1}^{N}[Y_{i}(1)-Y_{i}(0)].$

We can similarly define the sample average treatment effect for the treated
(SATT) and for the control (SATC) following the same principles as the
population estimands.

Another set of estimands focuses on the ATE conditional on covariate
distributions and emphasizes treatment effect heterogeneity in pre-specified
subgroups (Abadie and Imbens, 2002).

######  Definition 2.4.

(Conditional average treatment effect) We define the conditional average
treatment effect as:

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

If the treatment effect varies by values of $X$, then we say that the
treatment effect is modified by $X$, or there is treatment effect
heterogeneity (HTE) across values of $X$. As a function of $X$, some papers
refer to $\tau(x)$ as the HTE function(e.g. Yang et al., 2023).

The definitions given above calculate the treatment effect as the difference
in average counterfactual outcomes, identifying it as an additive effect.
However, the concept of treatment effect is more flexible than this
description suggests. Especially for binary outcomes, reasonable alternative
estimands include odds ratios and relative risk.

######  Definition 2.5.

(Average treatment effect as relative risk (RR) and odds ratio (OR))

$\displaystyle\tau^{RR}=\mathbb{E}Y(1)/\mathbb{E}Y(0),$
$\displaystyle\tau^{OR}=\frac{\mathbb{E}Y(1)/(1-\mathbb{E}Y(1))}{\mathbb{E}Y(0)/(1-\mathbb{E}Y(0))}.$

In this paper, we primarily focus on the ATE defined as the difference between
potential outcomes, denoted as $\tau=\mathbb{E}[Y(1)-Y(0)]$ in Definition 2.1.
Many methods discussed in Section 3 are developed and demonstrated using this
definition. In Section 3, we comment on how these methods might be adapted to
accommodate alternative ATE estimands.

###  2.3 Common identification assumptions

In this section, we discuss the common identification assumptions and
highlight two key challenges specific to the data fusion for efficiency gains.

######  Assumption 1.

Stable unit treatment value assumption (SUTVA) (Rubin, 1978)

(i) No interference: Unit $i$’s potential outcomes do not depend on other
units’ treatments.

$Y_{i}(a)\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}(A_{1},\dots A_{i-1},A_{i+1},\dots A_{n}).$

(ii) No multiple versions of treatment or control: the notion of being a
treated (or control) subject is well-defined, so the mechanism by which this
comes about does not alter the distribution of the outcome.

SUTVA is fundamental to the potential outcome framework, because it allows us
to write $Y=Y(a)$ if $A=a$, for $a\in\\{0,1\\}$. This means that an individual
with observed treatment $A$ equal to $a$, has observed outcome $Y$ equal to
their counterfactual outcome $Y(a)$. Under the decision-theoretic framework,
Dawid (2021) proposed SUTDA, which, instead of the value, assumes that the
distribution of $Y_{i}$ only depends on the treatment $A_{i}$ assigned to that
individual. SUTDA is a weaker condition than SUTVA and is sufficient for
causal inference.

######  Assumption 2.

(RCT ignorability) (i) Exchangeability: $Y(a)\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}A\,|\,S=1$ for $a\in\\{0,1\\}$, and (ii)
positivity: $0<P(A=a\,|\,X,S=1)<1$ for all $a\in\\{0,1\\}$ and $x$ in the
support of the RCT.

As the first part, exchangeability assumes that the counterfactual outcomes
are independent of the actual treatments, so that individuals in the treatment
and control group are ‘exchangeable’ with respect to their outcomes. A weaker
version is mean exchangeability,
$\mathbb{E}[Y(a)\,|\,A,S=1]=\mathbb{E}[Y(a)\,|\,S=1]$ or
$\mathbb{E}[Y(a)\,|\,A,X,S=1]=\mathbb{E}[Y(a)\,|\,X,S=1]$ for all $a$. The
second part, positivity, requires each participant to have a non-zero
probability of being assigned to treatment. Randomization is highly valued
because it breaks up the confounding effect on treatment assignment, so the
ignorability assumption holds for RCTs by design.

######  Assumption 3.

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

Assumption 3 states that the RCT is representative of the target population in
terms of ATE. For this review, we make this assumption to concentrate on
evaluating the methods for efficiency gains without considering
transportability issues—addressing one aspect at a time.

Under Assumptions 1–3, it is easy to see that the target causal estimand is
identifiable by the RCT as $\tau=\mathbb{E}[Y\mid A=1,S=1]-\mathbb{E}[Y\mid
A=0,S=1]$ and can be estimated through strategies of various complexity such
as difference in means, inverse propensity weighting (IPW) or augmented
inverse propensity weighting (AIPW) (Robins et al., 1994; Robins, 2000). These
assumptions are generally acceptable due to the controlled nature of RCTs. For
the target parameter to be identifiable in RWD, similar assumptions must be in
place. However, these are unlikely to be met in real-world settings,
presenting us with two key challenges of data fusion in achieving efficiency
gains.

######  Assumption 4.

(RWD comparability) $\mathbb{E}[Y(a)\,|\,X,S=0]=\mathbb{E}[Y(a)\,|\,X,S=1]$
for $a\in\\{0,1\\}$, or a weaker version:
$\mathbb{E}[Y(1)-Y(0)\,|\,X,S=0]=\mathbb{E}[Y(1)-Y(0)\,|\,X,S=1]$.

Assumption 4 posits that after adjusting for baseline covariates, the expected
potential outcomes in the RWD are equal to those in the RCT. This assumption
is reasonable when the RCT and RWD samples originate from the same super
population. For example, this can be reasonably assumed to hold when a
clinical trial and the RWD are collected from the same population around the
same time, with no significant differences in standard of care, treatment
adherence, or measurement of outcomes and covariates. The weaker version above
states that the contrasts of potential outcomes—or causal effects—are
equivalent, accommodating potential overall shifts in outcomes due to factors
like time trends or systematic differences in recording outcomes.

######  Assumption 5.

(RWD ignorability) $Y(a)\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}A\,|\,X,S=0$ and $0<P(A=a\,|\,X,S=0)<1$ for
all $a\in\\{0,1\\}$ and $x$ in the support of the RCT.

Assumption 5 states that the treatment assignment in the RWD is ignorable and
that, for individuals with the same covariates $X$ covering the support of the
RCT, there is a non-zero possibility of treatment in the real world. Along
with Assumption 4, the target causal estimate is identifiable by reweighting
the RWD against the RCT, as detailed in Colnet et al. (2022). However, in
reality, treatment selection is never random, and the circumstances of
individuals who seek specific treatments usually differ from those who do not,
resulting in selection bias. Therefore, while comparability might be
assumable, claiming treatment ignorability in the RWD without additional
external knowledge of the environment is impractical.

So, what makes Assumptions 4 and 5 key challenges of data fusion in achieving
efficiency gains? Firstly, the assumptions of comparability and treatment
exchangeability cannot be falsified using RWD alone. In the presence of an
RCT, potential violations can be tested. For instance, a significant
discrepancy between treatment effects estimated from the RWD and those from
the RCT indicates a likely violation. This principle underlies many data
fusion methods discussed in Section 3. Secondly, even when a violation is
detectable through comparison with the RCT, it remains unclear whether
Assumption 4 or 5 is breached. Morucci et al. (2023) present an interesting
discussion on this in their Theorem 1.

###  2.4 Identification and baseline estimator

To establish a baseline for assessing the benefits of data fusion compared to
using RCT data alone, we use the AIPW estimator (Robins et al., 1994;
Rotnitzky et al., 1998). Under Assumptions 1–3, the target estimand $\tau$ is
identifiable using the RCT alone:

$\tau=\mathbb{E}_{R}\left\\{\frac{A\left(Y-m_{1}(X)\right)}{e(X)}-\frac{(1-A)\left(Y-m_{0}(X)\right)}{1-e(X)}+(m_{1}(X)-m_{0}(X))\right\\}$

where $e(X)$ is the propensity score, i.e., $e(X)=\mathbb{E}[A=1\,|\,X]$, and
$m_{a}(X)$ is the conditional outcome mean, i.e.,
$m_{a}(X)=\mathbb{E}[Y\,|\,a,X]$. When Assumptions 4 and 5 hold, $\tau$ is
also identifiable in the RWD.

The AIPW estimator combines propensity score weighting with outcome regression
to provide a doubly robust method for estimating treatment effects, ensuring
consistent results even if one of the models is misspecified. We define
baseline estimators:

$\hat{\tau}_{r}=\frac{1}{n_{r}}\sum_{i\in\mathcal{R}}\left[\frac{A_{i}\left(Y_{i}-\hat{m}_{1,r}(X_{i})\right)}{\hat{e}_{r}(X_{i})}-\frac{(1-A_{i})\left(Y_{i}-\hat{m}_{0,r}(X_{i})\right)}{1-\hat{e}_{r}(X_{i})}+\left(\hat{m}_{1,r}(X_{r})-\hat{m}_{0,r}(X_{i})\right)\right]$

and

$\hat{\tau}_{o}=\frac{1}{n_{o}}\sum_{i\in\mathcal{O}}\left[\frac{A_{i}\left(Y_{i}-\hat{m}_{1,o}(X_{i})\right)}{\hat{e}_{o}(X_{i})}-\frac{(1-A_{i})\left(Y_{i}-\hat{m}_{0,o}(X_{i})\right)}{1-\hat{e}_{o}(X_{i})}+\left(\hat{m}_{1,o}(X_{i})-\hat{m}_{0,o}(X_{i})\right)\right],$

where $\hat{e}_{r}(X)$, $\hat{e}_{o}(X)$, $\hat{\mu}_{a,r}(X)$, and
$\hat{\mu}_{a,o}(X)$ are estimated propensity score and outcome in the RCT and
RWD, respectively. One attractive property of the AIPW is that it is
asymptotically normal. Under Assumptions 1–3, provided either of the the
propensity score model or the outcome regression model is correctly specified,
$\tau=\mathbb{E}(\hat{\tau}_{r})$. For the RWD, if Assumptions 4 and 5 hold,
with the same condition, $\tau=\mathbb{E}(\hat{\tau}_{o})$. However, if either
Assumption 4 or 5 is violated, then $\hat{\tau}_{o}$ is biased. We define
$\delta=\mathbb{E}(\hat{\tau}_{o})-\tau$ and refer to it as the “RWD bias” for
convenience.

##  3 Existing methods

This section provides a detailed review of recent data fusion techniques aimed
at efficiency gains. Lacking a single taxonomy of methods, we have organized
them based on their overarching philosophy. This means looking at the big-
picture strategy each method uses to bring different data sources together.

Alternative grouping criteria include whether methods make parametric
assumptions about the data generative model, whether they integrate external
data with both treatment and control arms or only the control arm, and whether
they focus on enhancing the CATE or directly on the ATE estimation. In the
latter part of this section, we present a table summarizing the various
characteristics of each method.

###  3.1 Test-then-pool

Test-then-pool style approaches start with the null hypothesis of equality
between the causal estimates from multiple sources and only pool if the null
hypothesis is not rejected. These approaches are commonly applied in meta-
analyses or multi-site studies, where results from several studies are
combined. Test-then-pool is the most basic form of dynamic borrowing (Viele et
al., 2014). Typically,test-then-pool approaches require the researcher to
specify the size of the hypothesis test, say $\alpha=0.025$, which offers a
way to cap the Type I error inflation introduced by the external data.
Essentially, test-then-pool is an “all or nothing” approach and a major issue
with such methods is that when the experimental data is small, the hypothesis
test for any discrepancy is under-powered and rarely gets rejected, resulting
in the observational data being pooled even if it is biased (Li et al., 2020).

The Elastic Integration approach proposed by Yang et al. (2023) aims at a
semi-parametric efficient estimation of the HTE function. They constructed a
test statistic to simultaneously falsify Assumption 4 and 5:

$T=\left\\{n_{o}^{-1/2}\sum_{i\in\mathcal{O}}\widehat{S}_{\mathrm{O},\widehat{\psi}_{\mathrm{E}}}\left(V_{i}\right)\right\\}^{\mathrm{T}}\widehat{\Sigma}^{-1}\left\\{n_{o}^{-1/2}\sum_{i\in\mathcal{O}}\widehat{S}_{\mathrm{O},\widehat{\psi}_{\mathrm{E}}}\left(V_{i}\right)\right\\}$

where $\widehat{S}_{\mathrm{O},\widehat{\psi}_{\mathrm{E}}}\left(\cdot\right)$
is the estimated semi-parametric efficiency score (SES) for an assumed HTE
function $\psi$, see Section 3.2 of Yang et al. (2020) for a detailed
definition. If Assumption 4 and 5 hold, then the test statistic $T$ should be
close to zero and a large $T$ indicates a violation. The authors proposed to
only use the RCT if $T$ exceeds a threshold $c_{\gamma}$, otherwise, pool the
two datasets together. They proposed a data-adaptive procedure to select the
threshold $c_{\gamma}$ to minimize MSE.

Dang et al. (2023c) proposed an Experiment-Selector Cross-Validated TMLE
(ESCV-TMLE) algorithm which splits the pooled data into folds, and employs a
cross-validation approach to data-adaptively select between two experiment
designs: RCT only, and pooled RCT + RWD, based on potential MSE reduction.
This can be broadly described as a test-then-pool approach, though it does not
involve specific hypothesis testing. According to their notation, the
‘selector’ operates as follows:

$s_{n}^{\star}=\underset{s}{\operatorname{argmin}}\frac{\hat{\sigma}_{D_{\Psi_{s}}^{*}}^{2}}{n}+\left(\hat{\Psi}_{s}^{\\#}\left(P_{n}\right)\right)^{2}.$

Intuitively, though maybe imprecisely, $s_{n}^{\star}$ is the ‘selector’ to
decide whether to pool the RWD.
$\frac{\hat{\sigma}_{D_{\Psi_{s}}^{*}}^{2}}{n}$ represents the asymptotic
variance derived from the efficient influence curve, and
$\left(\hat{\Psi}_{s}^{\\#}\left(P_{n}\right)\right)^{2}$ is the estimated
squared bias. Essentially, the decision to pool a certain fold of the RWD is
based on optimizing the MSE.

Moreover, when an negative control outcome (NCO) is available, they
innovatively proposed to include the difference in NCO into the squared-bias
term to improve the estimation of bias:

$s_{n}^{\star\star}=\underset{s}{\operatorname{argmin}}\frac{\hat{\sigma}_{D_{\Psi_{s}}^{*}}^{2}}{n}+\left(\hat{\Psi}_{s}^{\\#}\left(P_{n}\right)+\hat{\Phi}_{s}\left(P_{n}\right)\right)^{2}.$

Here, $\hat{\Phi}_{s}\left(P_{n}\right)$ represents the estimated difference
in the NCO between the treatment and control groups in the pooled RWD + RCT
dataset. In an RCT, or in RWD without confounding, if the NCO is valid,
$\Phi_{s}\left(P_{n}\right)$ should be zero. Thus, any deviation of
$\hat{\Phi}_{s}\left(P_{n}\right)$ serves as a proxy of hidden confounding.

For the inclusion of an NCO to function as intended, several conditions must
be met:

  1. (i)

U-compatibility: The unmeasured factors confounding the treatment-outcome
relationship are the same as those confounding the treatment-NCO relationship.
In other words, the unmeasured confounders are fully represented by the NCO,
with no additional confounders between the NCO and treatment.

  2. (ii)

Same direction: The unmeasured factors that affect the treatment-outcome
relationship should affect the treatment-NCO relationship in the same
direction. Otherwise the two confounding biases may cancel each other out.

Scale is another important consideration. Any significant difference in the
scale of the outcome and the NCO can alter the balance between the two biases.
For example, in a case study where the authors applied the ESCV-TMLE method to
study the effect of an anti-diabetic drug (Dang et al., 2023a), they used rate
of fracture at the NCO, which was around 0.7%, compared to a primary outcome
rate of about 4%. In such cases, the inclusion of an NCO might not influence
the trade-off by much.

In the simulations presented in Section 4, we examine the sensitivity of the
ESCV-TMLE estimator to the use of an NCO, as well as to violation of the
U-compatibility assumption.

The advantage of the ESCV-TMLE is its flexibility to combine either both arms
of the RWD or only the control arm. However, one limitation is the
interpretation of the estimand. For instance, if the selector combines the RWD
in 60% of cross-validation folds, it becomes challenging to interpret the
underlying population upon which the treatment effect is assessed.

###  3.2 Bayesian dynamic borrowing

Bayesian dynamic borrowing methods are widely applied to incorporate
historical studies to construct an informative prior for the outcome $Y$ or
the treatment effect $\tau$. The Bayesian framework offers a natural mechanism
to discount external information that may be in conflict with the trial data.
Ibrahim and Chen (2000) introduced a general power prior approach that raises
the likelihood of historical data to a power $\eta$:

$\pi_{\operatorname{PP}}\left(\theta\mid D_{O},\eta\right)\propto
L\left(\theta\mid D_{O}\right)^{\eta}\pi_{0}\left(\theta\right),$

where $\theta$ represents the parameters of interest, $L(\theta\,|\,\cdot)$
denotes the likelihood function, and $\pi_{0}$ is the initial prior, which is
often chosen to be uninformative. In this setup, $\eta$ is an important hyper-
parameter controlling the discounting of external evidence. It is assumed to
be known and can be set by expert opinion, which can be challenging in
practice. The authors extended this formulation with a hierarchical
specification where the uncertainty about $\eta$ is endogenized through a
prior distribution $\pi(\eta)$:

$\pi_{\operatorname{HPP}}\left(\theta,\eta\mid D_{O}\right)\propto
L\left(\theta\mid
D_{O}\right)^{\eta}\pi_{0}\left(\theta\right)\pi\left(\eta\right).$

However, Duan et al. (2006) and Neuenschwander et al. (2009) both cautioned
against the above formulation as it ignores a normalizing term, which is a
function of the unknown $\eta$ and thus violates the likelihood principle.
They showed that this inappropriate omission leads to counter-intuitive
pooling results where the influence of the external data is close to zero even
when the current and historical data are compatible. They therefore proposed
the normalized power prior:

$\pi_{\operatorname{NPP}}\left(\theta,\eta\mid
D_{O}\right)=\frac{L\left(\theta\mid
D_{O}\right)^{\eta}\pi(\theta)\pi(\eta)}{\int_{\Theta}L\left(\theta\mid
D_{O}\right)^{\eta}\pi(\theta)\mathrm{d}\theta}.$

A weakness of the power prior and normalized power prior approaches is that
the commensurability of the external and current data is not explicitly
parameterized. Ideally, we expect the external data is considered with more
weight if it appears more compatible to the current data. To this end, Hobbs
et al. (2011, 2012) introduced a location commensurate power prior (LCPP),
where the different parameters in the current and historical data are
explicitly modeled:

$\pi_{\operatorname{LCPP}}\left(\theta,\eta,\omega\mid
D_{O}\right)\propto\int\frac{\left[L\left(\theta_{o}\mid
D_{O}\right)\right]^{\eta}}{\int\left[L\left(\theta_{o}\mid
D_{0}\right)\right]^{\eta}\pi(\theta)d\theta_{o}}\times
N\left(\theta\mid\theta_{0},\frac{1}{\omega}\right)d\theta_{o}\times\operatorname{Beta}\left(\eta\mid
g(\omega),1\right)\times\pi(\omega).$

Here $\theta_{o}$ is the parameter in the historical data and $\omega$
represents the precision of how $\theta$ is distributed around $\theta_{o}$
and hence parameterizes the commensurability. $g(\cdot)$ is a positive and
increasing function of $\omega$, so that $\eta$ tends to take a large value
when there is strong evidence of commensurability. Extending the commensurate
prior approach, Schmidli et al. (2014) proposed a meta-analytic-predictive
posterior approach to borrow information from multiple historical trials,
using a mixture of standard priors.

Despite the advances, specifying a hyper-parameter or a commensurability
function remains a challenge and the applications are mainly to either
aggregate data or simple data generating models, for example in Bonapersona et
al. (2021). Lin et al. (2024) proposed a novel data-adaptive power likelihood
approach, enhancing the power prior framework in three aspects: i) introduced
a data adaptive method to select the learning parameter $\eta$; ii) utilized
the frugal parameterization (Evans and Didelez, 2024) to support more
sophisticated causal models for individual participant-level data; and iii)
adopted posterior approximation to overcome the computational challenges posed
by MCMC sampling for large datasets. The posterior takes the form of

$q_{\operatorname{PL}}(\theta;\eta)=q\left(\theta\mid
D_{O},D_{R},\eta\right)\propto L\left(\theta\mid D_{R}\right)L\left(\theta\mid
D_{O}\right)^{\eta}\pi_{0}\left(\theta\right)$

and $\eta$ is adaptively selected by optimizing the expected log predictive
density (ELPD) across the trial population $P_{R}$:

$\eta=\operatorname*{arg\,max}_{\eta\in[0,1]}\mathbb{E}_{P_{R}}\log
p_{\eta}\left(V\mid D_{O},D_{R}\right),$

where $p_{\eta}\left(v\mid
D_{O},D_{R}\right)=\int_{\Theta}p(v\mid\theta)q_{\operatorname{PL}}(\theta;\eta)d\theta$
is the posterior predictive distribution. Essentially, the learning rate
$\eta$ is selected such that the parameter posterior predicts the trial data
with the best expected accuracy. In addition to the three aforementioned
improvements on traditional Bayesian methods, the power likelihood approach
offers two other significant advantages. First, it provides the flexibility to
integrate either only external control or both arms of the observational data.
This flexibility is valuable in real-world applications because treatments,
such as new drugs or social policies, are often not available in RWD. Second,
compared to non-Bayesian methods, it automatically includes uncertainty
quantification through a credible interval. However, a common limitation of
Bayesian methods is the challenge in optimizing the inference of a particular
parameter of interest. Specifically, the power likelihood method selects
$\eta$ based on optimizing the joint posterior predictive distribution instead
focusing on outcome means or the ATE. This can be a drawback compared to
frequentist methods, such as targeted learning using TMLE (van der Laan and
Rubin, 2006), which directly targets the estimation of the ATE.

###  3.3 Weighted combination

Another approach to data fusion is to construct an estimator that is a
weighted average of $\hat{\tau}_{r}$ and $\hat{\tau}_{o}$, and the weight is
determined by optimizing a certain loss function. A common choice is MSE.

Tarima et al. (2020) and Oberst et al. (2022) proposed such a convex
combination estimator
$\hat{\tau}_{\operatorname{MM}}=(1-\lambda)\hat{\tau}_{r}+\lambda\hat{\tau}_{o}$
where $\lambda$ ranges from 0 to 1 and represents the weight given to the
treatment effect estimate from the observational data. $\lambda$ is chosen to
minimizing the MSE. That is,

$\displaystyle\lambda$
$\displaystyle=\operatorname*{arg\,min}_{\lambda\in[0,1]}\mathbb{E}\left[\left(\hat{\tau}_{\operatorname{MM}}-\tau_{0})^{2}\right)\right]$
$\displaystyle=\frac{\operatorname{Var}(\hat{\tau}_{r})-\operatorname{Cov}(\hat{\tau}_{r},\hat{\tau}_{o})}{\delta^{2}+\operatorname{Var}(\hat{\tau}_{r}-\hat{\tau}_{o})}.$

Since $\delta$ is unknown, the authors proposed to plug in the difference
between the RCT and observational data estimates, i.e.
$\hat{\delta}=\hat{\tau_{o}}-\hat{\tau_{r}}$ as an estimate of the bias, to
find an estimate of the MSE-minimizing weight $\hat{\lambda}$. Oberst et al.
(2022) gave a bound on the ratio of MSE over using RCT alone in their Theorem
2, which takes a maximum of $2.25$ when $\hat{\tau}_{r}$ and $\hat{\tau}_{o}$
are independent. It is worth noting that while $\hat{\tau}_{r}-\hat{\tau}_{o}$
is an unbiased estimator of the true bias $\delta$, the plug-in estimator
$\hat{\lambda}^{*}$ is not unbiased with regard to the MSE-minimizing weight
$\lambda^{*}$, due to Jensen’s Inequality. In fact, the plug-in estimator
$\hat{\lambda}$ underestimates $\lambda^{*}$ on average when
$\operatorname{Cov}(\hat{\tau}_{r},\hat{\tau}_{o})=0$ due to the convexity of
the expression. This means that when bias is low, their method fails to
include as much influence from the observational data as is theoretically
optimal.

Chen et al. (2021) introduced a dichotomous weighted average approach with a
preliminary hypothesis testing step to account for bias. The approach hinges
on a comparison of bias against a pre-determined threshold. If the bias is
below this threshold, the weight is selected to optimize the MSE, as though
there is no bias at all; otherwise, the integration is optimized with an
estimate of the bias, similar to the MSE-minimizing estimator in Oberst et al.
(2022). Their dichotomous procedure to estimate bias is as follows:

$\hat{\delta}_{\gamma}=\begin{cases}\operatorname{sgn}(\hat{\delta})\left(|\hat{\delta}|-\gamma\cdot\sqrt{\hat{\sigma}_{r}^{2}+\hat{\sigma}_{o}^{2}}\right)&\text{
if
}|\hat{\delta}|\geq\gamma\cdot\sqrt{\hat{\sigma}_{r}^{2}+\hat{\sigma}_{o}^{2}}\\\
0&\text{ otherwise. }\end{cases}$

The selection of $\gamma$, a critical hyper-parameter, impacts the extent of
bias adjustment, with a smaller $\gamma$ means that the bias term is less
aggressively shrunk to zero and thus is more conservative. The authors proved
that, under a certain data generating process, their estimator achieves a
minimax optimal rate, up to a poly-log factor, of convergence to the oracle,
where the true bias $\delta$ is unknown. However, the determination of
$\gamma$’s value, which significantly influences the combination and lacks a
predefined range, presents a practical challenge due to its arbitrary nature
and the difficulty in justifying its selection to stakeholders.

Two significant advantages of the two described linear combination methods
include: i) ease of implementation—requiring only $\hat{\tau}_{r}$, and
$\hat{\tau}_{o}$ and their variances, which can be obtained from estimators
like IPW. The methods only require straightforward calculations with minimal
coding and does not demand extensive computational resources; ii)
flexibility—the simplicity of these approaches allows for easy adaptation to
alternative estimands, such as the odds ratio, by merely adjusting the
definitions and estimators of $\tau_{r}$, and $\tau_{o}$. This inherent
flexibility makes the method adaptable to scenarios where the outcome is
binary.

Two limitations of the two methods are: i) uncertainty quantification—their
model-free nature implies that confidence intervals cannot rely on assumptions
inherent in parametric models. Neither set of authors made any asymptotic
claims, leading users to resort to bootstrap methods for inference; ii)
applicability–by design, the two methods can only apply to integrating RWD
with both treatment and control groups are available. This restricts their
utility in scenarios where treatment, such as a novel drug, is not yet
available in the market and absent in the RWD.

###  3.4 Shrinkage

Stein’s Paradox demonstrates that, for three or more parameters, there exists
a combined estimator that has a lower expected MSE than the traditional method
of estimating each parameter independently, even when the independent
estimators are unbiased (Stein, 1956; James and Stein, 1961). This result
means that we gain in terms of a lower MSE by simultaneously estimating
multiple parameters. The James-Stein estimator, is an application of this
paradox, which shrinks the independent MLEs towards the zero vector.

Green and Strawderman (1991) and Green et al. (2005) extended this idea to
combining unbiased and possibly biased estimators and proposed an estimator
that shrinks the unbiased estimators towards the biased estimators:

$\hat{\bm{\tau}}_{\operatorname{S1}}=\bm{\delta}_{2}^{+}(\hat{\bm{\tau}}_{\bm{r}},\hat{\bm{\tau}}_{\bm{o}})=\hat{\bm{\tau}}_{\bm{o}}+\left(\mathrm{I}_{K}-\frac{a\bm{\Sigma}_{r}^{-1}}{\left(\hat{\bm{\tau}}_{\bm{o}}-\hat{\bm{\tau}}_{\bm{r}}\right)^{T}\bm{\Sigma}_{r}^{-2}\left(\hat{\bm{\tau}}_{\bm{o}}-\hat{\bm{\tau}}_{\bm{r}}\right)}\right)_{+}\left(\hat{\bm{\tau}}_{\bm{o}}-\hat{\bm{\tau}}_{\bm{r}}\right),$

where $\hat{\bm{\tau}}_{\bm{r}}$ and $\hat{\bm{\tau}}_{\bm{o}}$ are now
vectors of dimension $K$, and each element is the ATE from one of $K$-strata
in the corresponding dataset. $\Sigma_{r}$ is the variance-covariance matrix
of $\bm{\hat{}}\tau_{r}$; note, that the components of
$\hat{\bm{\tau}}_{\bm{r}}$ (and $\hat{\bm{\tau}}_{\bm{o}}$) are assumed
mutually independent so $\bm{\Sigma_{r}}$ is a diagonal matrix. The shrinkage
parameter $a$ is customarily set as $K-2$. The plus sign denotes taking the
element-wise positive part of the matrix.

In the same vein, Rosenman et al. (2023a) proposed a general procedure for
deriving shrinkage estimators in this setting, by minimizing Stein’s Unbiased
Risk Estimate (SURE). They developed a new estimator,

$\hat{\bm{\tau}}_{\operatorname{S2}}=\hat{\bm{\tau}}_{\bm{o}}+\left(\bm{I}_{K}-\frac{\operatorname{Tr}\left(\bm{\Sigma}_{r}^{2}\right)\bm{\Sigma}_{r}}{\left(\hat{\bm{\tau}}_{\bm{o}}-\hat{\bm{\tau}}_{\bm{r}}\right)^{T}\bm{\Sigma_{r}}^{2}\left(\hat{\bm{\tau}}_{\bm{o}}-\hat{\bm{\tau}}_{\bm{r}}\right)}\right)_{+}\left(\hat{\bm{\tau}}_{\bm{o}}-\hat{\bm{\tau}}_{\bm{r}}\right).$

Both methods require at least three strata and independence among stratum
estimates, to ensure a reduction in overall MSE. Simulations by Rosenman et
al. (2023a) revealed that the choice of stratification variables and the
number of strata affect performance. However, their findings remain
inconclusive: it is unclear what the optimal number of strata is, whether it
is better to have more or fewer strata, and no single estimator consistently
outperforms the others.

For inference, both papers proposed empirical Bayes (EB) confidence intervals
for the $k$th stratum-specific estimate, $\hat{\tau}_{\operatorname{S2},k}$,
based on the result of the James-Stein estimator variance in Morris (1983). A
key difference between EB CIs and conventional frequentist CIs is that EB CIs
aim to achieve an average coverage rate for the entire $\hat{\bm{\tau}}$,
which means that they may under-cover the causal effects for some strata and
over-cover for others. Simulations in Rosenman et al. (2023a) confirmed that
for large sample sizes, their 90% EB CIs met the desired average coverage
rate, although the coverage rates varied across individual strata.

Combining the weighted-combination and shrinkage concepts, Rosenman et al.
(2023b) constructed a double shrinkage approach, which intuitively speaking,
involves first constructing a MSE-minimizing weighted average of
$\hat{\bm{\tau}}_{\bm{o}}$ and $\hat{\bm{\tau}}_{\bm{r}}$, and then shrinking
it towards zero. To derive the estimator, the authors introduced a Bayesian
hierarchical model for the ‘true’ value of $\bm{\tau}$, the bias $\bm{\delta}$
and the estimates $\hat{\bm{\tau}}_{\bm{o}}$ conditional on $\bm{\tau}$ and
$\bm{\delta}$. The expression for the $k$th entry of their double-shrinkage
estimator is as follows:

$\displaystyle\tau_{DS,k}\left(\phi^{2},\omega^{2}\right)$
$\displaystyle=\underbrace{\left\\{\frac{\omega^{2}\left(\phi^{2}+\sigma_{ok}^{2}+\sigma_{rk}^{2}\right)}{\sigma_{rk}^{2}\left(\phi^{2}+\sigma_{ok}^{2}\right)+\omega^{2}\left(\phi^{2}+\sigma_{ok}^{2}+\sigma_{rk}^{2}\right)}\right\\}}_{a_{k}}\cdot\Bigg{\\{}\underbrace{\frac{\left(\phi^{2}+\sigma_{ok}^{2}\right)}{\phi^{2}+\sigma_{ok}^{2}+\sigma_{rk}^{2}}}_{\lambda_{k}}\hat{\tau}_{rk}+\underbrace{\frac{\sigma_{rk}^{2}}{\phi^{2}+\sigma_{ok}^{2}+\sigma_{rk}^{2}}}_{1-\lambda_{k}}\hat{\tau}_{ok}\Bigg{\\}},$

where $\phi$ and $\omega$ are the prior variances of $\bm{\tau}$ and
$\bm{\delta}$, respectively. The second term resembles the form of a weighted-
combination estimator, while the first term, $a_{k}$, acts as a shrinkage
factor towards zero. Notably, $\phi$ and $\omega$ are crucial in regulating
the combination, and since they are unknown, the authors suggested three
estimation methods: moment matching, empirical Bayes MLE, and URE
minimization. Their simulations indicated that estimating $\phi$ and $\omega$
using moment matching and EB MLE generally outperform URE minimization. Like
the findings in Rosenman et al. (2023a), simulations indicated that the
performance of the double-shrinkage estimator is influenced by the selection
of stratification variables and the number of strata, though the mechanism
behind this influence remains unclear.

###  3.5 Bias correction

Another idea is to use the RCT data, which offers unbiased causal estimates,
to correct the bias in the RWD estimate. The two methods in Kallus et al.
(2018) and Yang et al. (2020) are in the same vein, and the latter is a
generalization of the former. By contrasting the RCT and RWD, these methods
learn the impact of hidden confounding, that is, violation of Assumption 5.
Both of the methods assume RWD comparability (Assumption 4) to hold.

Kallus et al. (2018) aimed at learning the impact of confounding that results
in the CATE function estimated from the RWD to be biased. That is,

$\displaystyle\xi(X)$ $\displaystyle=\\{\mathbb{E}[Y(1)\mid
X,S=0]-\mathbb{E}[Y(1)\mid X,A=1,S=0]\\}$
$\displaystyle-\\{\mathbb{E}[Y(0)\mid X,S=0]-\mathbb{E}[Y(0)\mid
X,A=0,S=0]\\},$

and when Assumption 5 holds, $\xi(X)=0$. Kallus et al. (2018) proposed to
obtain samples of $\xi(x)$ by calculating the difference between the CATE
estimated in the RCT and RWD. As an illustration, they assumed a linear
structure of $\xi(x)$ and derived the least-squares coefficient estimates. The
resulting $\hat{\xi}(x)$ serves as the bias-correction term in this procedure
termed experimental grounding. They demonstrated that under linearity
assumption, the bias-corrected CATE estimate
$\hat{\tau}_{\operatorname{XG}}(x)$ is root-$n$ consistent, although
efficiency gains are not theoretically proven.

Yang et al. (2020) summarized the impact of unmeasured confounding on the
potential outcomes into the confounding function:

$\varphi(X)=\mathrm{E}[Y(0)\mid A=1,X,S=0]-\mathrm{E}[Y(0)\mid A=0,X,S=0].$

They demonstrated that, assuming the parametric forms of the confounding
function $\varphi(X)$ and the CATE function $\tau(X)$ are known, both
$\varphi(X)$ and $\tau(X)$ are non-parametrically identifiable. They also
derived semi-parametric efficient, doubly robust estimators under several
examples of assumed parametric form of $\tau(X)$. In their Theorem 3, the
authors proved efficiency gains in estimating $\tau(x)$ compared to using the
RCT alone, provided the terms in parametric forms of $\tau(x)$ and
$\varphi(x)$ are not co-linear. However, if the terms of $\tau(x)$ are a
linear combination of those in $\varphi(x)$, their is no efficiency gain from
data fusion. This is empirically observed in our simulations. The authors
provided some intuition for this result: when the terms are linearly
dependent, all the observational data is consumed to estimate the coefficients
of $\varphi(X)$ and hence does not contribute to the estimation of $\tau(x)$.
While their method is demonstrated using the additive treatment effect as
defined in Definition 2.1, adapting it to binary outcomes, or alternative
estimands s such as relative risks or odds ratios may present additional
complexities.

While these two methods do not require RWD ignorability (Assumption 5), these
methods make assumptions about the parametric structure of the functional
forms of CATE and confounding, which are subject to misspecification.

van der Laan et al. (2024) proposed a bias-correction method that targets a
slightly different estimand: the average treatment effect averaged over the
covariate distribution of the combined RCT-RWD data. They decomposed the the
the ATE estimand into the difference between a pooled-ATE estimand and a bias
estimand that accounts for the conditional effect of RCT enrollment on the
outcome. To estimate them, they introduced an adaptive targeted minimum loss-
based estimation (A-TMLE) framework. Unlike the other two bias-correction
methods discussed, the A-TMLE estimator does not rely on parametric models for
the biases but instead employs flexible methods for estimation.

###  3.6 Prognostic adjustment

The EMA recently approved the historical borrowing method, PROCOVA, by Schuler
et al. (2021) (European Medicines Agency, 2022). The method enhances the
control arm of RCTs by incorporating a prognostic score fitted using a
historical dataset, which is required to originate from the same population as
the trial control arm, to ensure that mean control outcomes are the same—this
is the stronger version of Assumption 4. Although the authors suggested that
this assumption can be weakened, the anticipated theoretical efficiency gain
relies on the control outcome estimated from the historical data being
consistent. This may not hold if the historical RWD are sourced from different
trials or countries. In addition, if outcomes are linearly related to
covariates, including linearity after transformation such as in generalized
linear models (GLMs), no efficiency benefit is achievable beyond that from the
routine adjustment by pre-treatment covariates. Coincidentally, this resembles
the intuition of the results in Yang et al. (2020) as discussed in Section
3.5; this is because all the information from the raw covariates is already
exploitable by the linear model.

The key advantage of PROCOVA is protection against Type I error, which is
crucial for regulatory decisions. The authors provide a sample size
calculation method based on the correlation between the prognostic score and
outcome residuals after adjusted for pre-treatment covariates. Despite relying
on assumptions such as equal correlation between the outcome and the
prognostic score, and equal variance of the outcome between the RWD and the
RCT, the ability to perform sample size calculation enables the evaluation of
potential efficiency gains, and thus the reduction in required samples size,
before conducting an RCT. The dual advantages of safeguarding against Type I
error and enabling sample size estimation distinguish PROCOVA as the only data
fusion method endorsed by regulatory agencies to date.

###  3.7 A summary table for application

Table 1 provides a practical summary of the methods discussed above, serving
as a reference for researchers selecting an appropriate method. The “data
assumption” column indicates whether a method is parametric (P), non-
parametric (NP) or semi-parametric (SP). Understanding these assumptions helps
users recognize their limitations: semi-parametric methods require correct
model specification of the parametric part to avoid biased estimates, while
non-parametric methods offer flexibility but may have poor finite sample
performance. The next two columns indicate whether a method is applicable to
combining only the control arm or both arms from the RWD. The “inference”
column describes how each approach quantifies uncertainty. The following
column assesses the difficulty of implementation, specifically coding
difficulty. “Easy” denotes straightforward coding in a couple of lines, while
methods like the power likelihood approach are marked as “medium” due to the
need for users to specify a causal model and assemble basic functions despite
the available package. The R packages available are: causl (Evans, 2021) and
ManyData (Lin, 2023) for the power likelihood approach,
ElasticIntegrative(Yang and Gao, 2023) for the elastic integrative estimator,
EScvtmle (Dang et al., 2023b) for the ESCV-TMLE estimator, and
IntegrativeHTEcf (Yang, 2020) for the confounding function method. See
Appendix A for details. The “alternative estimands” column compares the
complexity of adapting methods to estimands such as relative risk. Semi-
parametric methods typically require significant modification, such as
replacing the efficient influence function, and current packages often lack
this flexibility. The final column lists primary limitations for users to
consider.

Category |  Method |  Data assumption |  RWD control-arm only |  RWD both arms |  Inference |  Implement. difficulty |  Alt. estmnds |  Limitations/requirements
---|---|---|---|---|---|---|---|---
Bayesian dynamic borrowing |  Power likelihood (Lin et al., 2024) |  P |  ✓ |  ✓ |  Bayesian credible interval |  Medium |  Easy |  Requires full parameterization of data distribution and correct HTE specification
Test-then-pool |  Elastic integration (Yang et al., 2023) |  NP |  |  ✓ |  Asymptotics/ bootstrap |  Package ready |  Hard |  Requires correct HTE specification
|  ESCV-TMLE (Dang et al., 2023c) |  SP |  ✓ |  ✓ |  Asymptotics |  Package ready |  Hard |  Interpretability is challenging
Weighted combination |  Anchored thresholding (Chen et al., 2021) |  NP |  |  |  Bootstrap |  Easy |  Easy |  No reversion to RCT-only with increasing bias
|  MSE-minimizing (Oberst et al., 2022) |  NP |  |  ✓ |  Bootstrap |  Easy |  Easy |
Bias correction |  Experiment grounding (Kallus et al., 2018) |  SP |  |  ✓ |  bootstrap |  Easy |  Hard |  Assumes linearity; Efficiency gains in HTE but not in ATE estimation
|  Confounding function (Yang et al., 2020) |  SP |  |  ✓ |  Asymptotics/ bootstrap |  Package ready |  Hard |  Requires correct HTE specification; Efficiency gains in HTE but not in ATE estimation
Covariate adjustment |  PROCOVA (Schuler et al., 2021) |  P |  ✓ |  |  Sandwich estimator |  Easy |  Easy |  Limited efficiency gain when outcome is nearly linear
Shrinkage |  Green et al. (2005) |  NP |  |  ✓ |  Bootstrap |  Easy |  Easy |  Requires multiple independent strata estimates
|  Rosenman et al. (2023a) |  NP |  |  ✓ |  Bootstrap |  Easy |  Easy |  Requires multiple independent strata estimates
Table 1: Comparison of different statistical methods.

##  4 Simulation studies

In our simulation, we aim to be as realistic as possible by introducing both
binary and continuous covariates, explicitly specifying hidden confounding,
and accounting for treatment heterogeneity, among other factors. Details of
method implementation can be found in Appendix A.

###  4.1 Simulation setup

We adopt the frugal parameterization proposed in Evans and Didelez (2024) to
parameterize and simulate from specified marginal structural models (MSM). The
primary advantage is its flexibility: we can simulate exactly the causal
effects we would like. We first generate three pre-treatment covariates
$X_{1}$, $X_{2}$ and $X_{3}$, where $X_{3}$ acts as an effect modifier to
introduce treatment effect heterogeneity. We assume that the marginal
distribution of $X_{1}$ is a standard normal distribution, and the marginal
distribution of $X_{2}$ follows a Bernoulli distribution with a probability of
0.5. We model the dependency between $X_{1}$ and $X_{2}$ later using a copula.
We assume that the effect modifier $X_{3}$ is independent of $X_{1}$ and
$X_{2}$ and is normally distributed.

$\displaystyle X_{1}$ $\displaystyle\sim\mathcal{N}(0,1)$ $\displaystyle
X_{2}$ $\displaystyle\sim\operatorname{Bernoulli}(0.5)$ $\displaystyle X_{3}$
$\displaystyle\sim\mathcal{N}(0,1)$ $\displaystyle U$
$\displaystyle\sim\mathcal{N}(0,1.1)$ $\displaystyle Y(a)\mid X_{3},U$
$\displaystyle\sim\mathcal{N}(0.5+0.2a+0.1X_{3}a+X_{3}+U,1)$

Given the mixture of continuous and discrete variables, we model the
dependence between $X_{1}$, $X_{2}$ and $Y$ using a latent Gaussian copula
model, as proposed by Fan et al. (2017). Specifically, we employ use a
trivariate Gaussian copula with correlation coefficients
$\rho_{X_{1}X_{2}}=2\operatorname{expit}(0.1)-1\approx 0.05$,
$\rho_{X_{1}Y}=2\operatorname{expit}(2)-1\approx 0.76$ and
$\rho_{X_{2}Y}=2\operatorname{expit}(1)-1\approx 0.46$. These coefficients
allow us to simulate a weak correlation between $X_{1}$ and $X_{2}$, a strong
influence of the continuous variable $X_{1}$ on the outcome and a moderate
correlation between the binary variable $X_{2}$ and $Y$.

We simulate treatment assignments in the RCT and in the observational data as
below:

$\displaystyle A_{obs}$
$\displaystyle\sim\operatorname{Bernoulli}(\operatorname{expit}\\{-0.5+X_{1}+X_{2}+X_{3}+\psi
U\\})$ $\displaystyle A_{RCT}$
$\displaystyle\sim\operatorname{Bernoulli}(0.5).$

This parameterization assumes that treatment assignment is completely
randomized in the RCT. In this setup, the value of the target estimand,
$\tau=\mathbb{E}[Y(1)-Y(0)]$ is $0.25$. In contrast, in the observational
data, treatment assignment depends on observed covariates $X_{1}$, $X_{2}$ and
$X_{3}$ as well as confounding from $U$. The degree of hidden confounding is
controlled by the value of $\psi$. When $\psi=0$, $U$ is no longer a
confounder, making the target estimand identifiable in the observational data.
However, as the absolute value of $\psi$ increases, the bias in the ATE
estimated from the RWD correspondingly intensifies. In the simulation, we vary
the value of $\psi$ to introduce different magnitudes of hidden confounding
into the observational data. Note that the average proportion of treated
individuals remains at 50% in all simulations. We set the size of the RCT and
observational datasets to 300 and 1200 respectively.

###  4.2 Augmenting both arms of the RWD

In Figure 2, each line represents one data fusion estimator discussed in
Section 3, and the results are averaged over 300 different simulated datasets.
Comparing the root mean square error (RMSE) of the fusion estimator
$\hat{\tau}$ to that of $\hat{\tau}_{r}$, the estimator from the RCT alone, is
particularly useful to determine the benefits of data fusion. The y-axis of
the figure shows the relative RMSE, which is defined similarly to Oberst et
al. (2022):

$\displaystyle
r\operatorname{RMSE}=\frac{\operatorname{RMSE}(\hat{\tau})}{\operatorname{RMSE}(\hat{\tau}_{r})},$

where as a reference estimator, $\hat{\tau}_{r}$ is the AIPW estimator using
the RCT data only. The fusion estimator $\hat{\tau}$ outperforms $\tau_{r}$ if
rRMSE is less than 1 and worse, if rRMSE exceeds 1, which means that the
efficiency gain is not enough to offset the inflation in bias.

As a consequence of violation of Assumption 4 and 5, the RWD bias
$\delta$—defined in Section 2.4—is plotted on the x-axis. Instead of the
absolute value of $\delta$, we plot bias on a relative scale. Specifically, we
define the relative bias as the multiple of the standard deviation of the RCT-
only estimator $\hat{\tau}_{r}$:

$\displaystyle
r\operatorname{Bias}=\frac{\delta}{\operatorname{sd}(\hat{\tau}_{r})}.$

For instance, a relative bias of one indicates that $\tau_{o}$ is one standard
deviation away from $\tau$. The advantage of using relative bias is that it
provides an intuitive metric of confounding bias in relation to the
variability of $\hat{\tau}_{r}$.

We also include an oracle performance curve for reference, which illustrates
the best possible reduction in RMSE that could be achieved if the bias of the
estimator $\hat{\tau}_{o}$ were known when fusing the two datasets. An example
is when $r\operatorname{Bias}=0$, indicating that the observational data
should be integrated entirely and treated equally as the RCT data. In this
scenario, the reduction in RMSE is due to the increase in total sample size,
which is approximately equal to $\sqrt{300/(1200+300)}\approx 0.44$ in our
simulation setting. We include this curve to demonstrate how close each method
can come to the optimal performance. However, note that the oracle curve does
not apply to the PROCOVA estimate as it only incorporates the control arm.

Figure 2: RMSE relative to RCT-only approach (without data fusion) by varied
magnitude of bias induced by hidden confounding, in the baseline scenario with
a continuous outcome Figure 3: Coverage rate and relative CI length of data
fusion estimators, by varied magnitude of bias induced by hidden confounding.

Figure 3 compares the coverage rates of 95% confidence intervals (credible
interval for the power likelihood estimator), and the ratio of the length of
the 95% CI of each data fusion estimator against that of $\hat{\tau}_{r}$.
This comparison illustrates a trade-off: estimators that achieve greater power
gains when bias is small are also at risk of poorer worst-case coverage rates.
This highlights the trade-off between power gains and coverage rate.

####  4.2.1 Anchored thresholding

As $\delta$ increases, almost all the curves converge to the reference line,
except for the anchored thresholding estimator in Chen et al. (2021). The
reason is that as the RWD bias, $\delta$, approaches infinity, their proposed
estimator becomes asymptotically biased. We illustrate this using the
following example.

######  Example 1.

Let $\hat{\tau}_{r}\sim N(\tau,\sqrt{\sigma^{2}/n_{e}})$ and
$\hat{\tau}_{o}\sim N(\tau+\delta,\sqrt{\sigma^{2}/n_{o}})$. We assume that
the total sample size is $n=n_{r}+n_{o}$ and $n_{o}=\rho n$, where
$\rho\in(0,1)$.

Following the methodology outlined in Section 3.2 of Chen et al. (2021), the
weighting factor becomes
$\omega=\frac{\sigma^{2}/n_{r}}{\sigma^{2}/n_{r}+\sigma^{2}/n_{o}}=\rho$, and
the estimator of bias $\hat{\delta}=\hat{\tau}_{o}-\hat{\tau}_{r}\sim
N\left(\delta,\sqrt{\frac{\sigma^{2}}{n_{r}}+\frac{\sigma^{2}}{n_{o}}}\right)$.
We assume $\delta$ is a very large positive number, specifically,
$\delta\gg\sqrt{\frac{\sigma^{2}}{n_{r}}+\frac{\sigma^{2}}{n_{o}}}$. After the
proposed soft-thresholding procedure, the adjustment for bias they proposed is
$\hat{\delta}_{\gamma}=\hat{\delta}-\gamma\sqrt{\frac{\sigma^{2}}{n_{r}}+\frac{\sigma^{2}}{n_{o}}}\sim
N\left(\delta-\gamma\sqrt{\frac{\sigma^{2}}{n_{r}}+\frac{\sigma^{2}}{n_{o}}},\sqrt{\frac{\sigma^{2}}{n_{r}}+\frac{\sigma^{2}}{n_{o}}}\right)$.
Recall that $\gamma$ is an arbitrary positive hyper-parameter that controls
the adjustment. Therefore, their final estimator becomes:

$\displaystyle\hat{\tau}_{\operatorname{AT}}$
$\displaystyle=\omega\hat{\tau}_{r}+\rho\left(\hat{\tau}_{o}-\hat{\delta}_{\gamma}\right)$
$\displaystyle=(1-\rho)\hat{\tau}_{r}+\rho\left(\hat{\tau}_{r}+\hat{\delta}-\hat{\delta}+\gamma\sqrt{\frac{\sigma^{2}}{n_{r}}+\frac{\sigma^{2}}{n_{o}}}\right)$
$\displaystyle=\hat{\tau}_{r}+\rho\gamma\sqrt{\frac{\sigma^{2}}{n_{r}}+\frac{\sigma^{2}}{n_{o}}}$
$\displaystyle\sim
N\left(\tau+\rho\gamma\sqrt{\frac{\sigma^{2}}{n_{r}}+\frac{\sigma^{2}}{n_{o}}},\sigma/\sqrt{n_{r}}\right).$

It is obvious that when $\delta$ is very large, their anchored-thresholding
estimator $\hat{\tau}_{AT}$ is biased by
$\rho\gamma\sqrt{\frac{\sigma^{2}}{n_{r}}+\frac{\sigma^{2}}{n_{o}}}$, with a
variance equal to using the RCT data alone. Additionally, a larger $\gamma$
leads to greater bias, which explains the empirical observations made in
Oberst et al. (2022). This extra bias is the reason why the RMSE curve remains
constantly above the reference line, and the coverage rate remains low as
$\delta$ increases. This means that the anchored thresholding estimator does
not offer the assurance that when the RWD is highly biased, we will not be
worse off than using the RCT alone. For this reason, despite showing the
highest power gains in Figure 3, we do not recommend the anchored thresholding
method for data fusion.

####  4.2.2 Bias-correction methods

Figure 4 focuses on the two bias-corrections methods. It shows that neither
the confounding function (Yang et al., 2020) nor the experimental grounding
(Kallus et al., 2018) reduces the RMSE, regardless of the level of confounding
bias. In fact, using the experimental grounding estimator results in a loss of
efficiency, leading to wider confidence intervals.

This is because both methods essentially distribute the difference between
$\hat{\tau}_{r}$ and $\hat{\tau}_{o}$ across a pre-specified function space
defined by the covariates. When marginalized over the RCT samples, the
estimate reverts to $\hat{\tau}_{r}$, yielding no efficiency gains in
estimating the ATE. van der Laan et al. (2024) further formalized this point
by comparing the non-parametric efficiency bounds, showing that if the target
parameter measures the HTE in the RCT averaged over the covariate distribution
of the RCT, their may still be no efficiency gains, even if the estimand is
efficiently estimated. Additionally, for the experiment grounding estimator,
the loss is due to the added noise, which makes the estimation less efficient.

It is worth noting that the confounding function method is designed for
efficiency gains in estimating the HTE function, with theoretical proof
provided in Theorem 3 in Yang et al. (2020). However, this method is not
suitable for estimating the ATE in our context. If we make an extra assumption
that the observational sample is representative of the target population,
meaning both the RCT and RWD comes from the same population, we can
marginalize the estimated HTE function over the pooled RCT + RWD data to
achieve power gains by effectively increasing the sample size. Even then, the
efficiency gain relies on correct specification of the function spaces.

Figure 4: Relative RMSE, coverage rate and relative CI length of bias-
correction estimators (Kallus et al., 2018; Yang et al., 2020).

####  4.2.3 Methods showing a trade-off

Figure 5: Relative RMSE, coverage rate and relative CI length of estimators
that exhibit a trade-off (Green et al., 2005; Rosenman et al., 2023a; Oberst
et al., 2022; Lin et al., 2024; Dang et al., 2023c; Yang et al., 2023).

Figure 5 focuses on estimators that exhibit a trade-off pattern, showing a
reduction in RMSE in some scenarios and an increase in others. Specifically,
the leftmost panel of Figure 5 demonstrates that when the RWD bias is small,
typically less than twice the standard deviation of $\hat{\tau}_{r}$, the
reduction in variance outweighs the introduced bias leading to an reduction in
RMSE. However, as $\delta$ exceeds two standard deviations, the RMSE increases
compared to not incorporating RWD. This is primarily due to the bias
overwhelming the benefit of reduced variance. As the RWD bias continues to
increase, most methods gradually converge to the reference line, indicating
that these methods recognize the extremely large RWD bias and respond by down-
weighting the influence of the RWD to zero, relying solely on the RCT for
estimation.

The middle and the rightmost panels of Figure 5 provide further details about
the trade-off. Since the trade-off in RMSE involves a trade-off between
squared bias and variance, it manifests in the coverage rate of the 95%
confidence interval and the relative length of these intervals. Additionally,
there is a general trade-off between best-case reward and worst-case risk.
Methods with larger potential efficiency gains are also more susceptible to
higher worse-case RMSE. This is evident in Figure 5, where the order of lines
from high to low before crossing the reference line is exactly reversed
afterwards. This implies that aiming for a larger efficiency gain involves
taking on greater risks.

Beyond the presented simulations, we have conducted additional simulations
with varying parameter values. The relative behaviour remains somewhat
consistent across these simulations. The power likelihood (Lin et al., 2024),
elastic integration (Yang et al., 2023), and the shrinkage estimator by
Rosenman et al. (2023a) are among the “ambitious” methods with substantial
potential efficiency gains but lower coverage rate in certain ranges of RWD
bias. In contrast, the MSE-minimizing and ESCV-TMLE (no NCO) estimators tend
to be “‘conservative”; they offer a smaller potential reduction in uncertainty
but maintain a coverage rate closer to the nominal rate.

In this simulation setting, when the RWD is unconfounded, the power likelihood
estimator achieves considerable efficiency gains, nearly reducing the RMSE to
the best possible level and shortening the CI length by almost 45%. However,
this comes with the risk of CI under-coverage, with the worst-case coverage
rate dropping below 60%. As discussed in Section 3.2, unlike frequentist
methods that directly target the inference of ATE, the power likelihood
approach selects the learning factor $\eta$ based on optimizing the expected
predictive density over the target population. This means $\eta$ is chosen so
that the $\eta$-indexed joint posterior of the model parameters best predicts
the ‘true’ data. Since this involves a joint distribution, $\eta$ is
influenced not only by the ATE, which is primarily a function of the mean
outcomes, but also the other features of the outcome distribution. For
example, if the conditional variance of $Y$ differs between the RCT and the
RWD—which is likely given the highly controlled environment of the RCT versus
the numerous unmeasured uncertainties in RWD—the power likelihood method
accounts for this difference. In such scenarios, stronger unmeasured
confounding potentially increases the noise in $Y$. Consequently, the power
likelihood approach typically selects a smaller $\eta$, rendering a more
cautious inclusion of information from the RWD. The simulations in Lin et al.
(2024) illustrated a scenario where the conditional volatility of $Y$
intensifies with the strength of unmeasured confounding. Their results showed
that in such cases, the power likelihood method achieves efficiency gains
while maintaining a coverage rate close to the nominal.

Key to the MSE-minimizing estimator is the weight $\lambda$ that balances the
RCT and RWD estimates. As discussed in Section 3.3, plug-in estimator
$\hat{\lambda}$ underestimates the MSE-minimizing $\lambda$ due to Jensen’s
inequality. Consequently, less RWD is included than is theoretical optimal.
This likely contributes to the estimator being more conservative compared to
other methods.

The learning factor $\eta$ in the power likelihood estimator and the weight
$\lambda$ in the MSE-minimizing estimators play crucial roles in moderating
the influence from the RWD. Both methods proposed ways to robustly select
these parameters. Indeed, Figure 6 demonstrates that both methods effectively
downweight the RWD in response to increased bias.

Figure 6: The average learning factor $\eta$ selected by power likelihood
method (Lin et al., 2024) and the average weight $\lambda$ assigned to
$\hat{\tau}_{o}$ selected by the MSE-minimizing method (Oberst et al., 2022).
Mote that due to differences in methodologies, the absolute levels of these
two parameters are not directly comparable.

Regarding the shrinkage estimators in Green et al. (2005) and Rosenman et al.
(2023a), as discussed in Section 3.4, these methods are designed to improve
the efficiency of the vector of strata estimates rather than directly
enhancing the ATE. To get an estimate of the target estimand $\tau$, we
introduce an extra step to multiply by a weight vector:

$\displaystyle\hat{\tau}_{\operatorname{ATE},S1}$
$\displaystyle=\bm{w}^{T}\hat{\bm{\tau}}_{S1}$
$\displaystyle\hat{\tau}_{\operatorname{ATE},S2}$
$\displaystyle=\bm{w}^{T}\hat{\bm{\tau}}_{S2},$

where $\bm{w}^{T}$ is a $1\times K$ vector with weights for each stratum in
the target population. We use $\hat{w}_{k}=n_{r,k}/n_{r},k=1,\dots K$ as an
estimate.

Theoretically, these shrinkage estimators should reduce the overall RMSE,
where “overall” refers to the sum, or a weighted sum, of strata estimator
RMSEs. However, this theoretical guarantee of RMSE reduction does not
necessarily extend to $\hat{\tau}_{\operatorname{ATE},S1}$ and
$\hat{\tau}_{\operatorname{ATE},S2}$. This explains why, in certain ranges of
$\delta$, neither the RMSE of the ATE is not consistently reduced, nor the
length of CIs. Furthermore, as discussed in Section 3.4, Green et al. (2005)
and Rosenman et al. (2023a) both suggested using the EB confidence interval
for each entry of the vectors $\hat{\bm{\tau}}_{S1}$ and
$\hat{\bm{\tau}}_{S2}$. However, since our focus is on $\hat{\tau}_{ATE,S1}$
and $\hat{\tau}_{ATE,S2}$, which are linear combinations of the strata
shrinkage estimates, the point-wise CIs are not very useful. It is more
appropriate to use the non-parametric bootstrap approach, which involves
resampling each stratum in the RCT and RWD, respectively, and calculating the
confidence intervals based on, in this implementation, 200 bootstrap samples.

####  4.2.4 ESCV-TMLE (Dang et al., 2023c)

The ESCV-TMLE method offers the option to incorporate an NCO to enhance the
estimation of bias. As discussed in Section 3.1, the NCO-integrated estimator
relies on the U-compatibility assumption. To investigate the sensitivity to
this assumption, we modify the simulation parameterization by splitting the
unmeasured confounder $U$ into two independent variables,
$U_{1}\sim\mathcal{N}(0,1)$ and $U_{2}\sim\mathcal{N}(0,0.1)$, such that
$U_{1}+U_{2}$ is equivalent to $U$. Specifically,

$\displaystyle Y(a)\mid X_{3},U_{1},U_{2}$
$\displaystyle\sim\mathcal{N}(0.5+0.2a+0.1X_{3}a+X_{3}+U_{1}+U_{2},1)$
$\displaystyle A_{obs}$
$\displaystyle\sim\operatorname{Bernoulli}(\operatorname{expit}(-0.5+X_{1}+X_{2}+X_{3}+\psi
U_{1}+\psi U_{2})).$

We further include three NCO variables:

$\displaystyle N_{1}$ $\displaystyle\sim\mathcal{N}(1+X_{2}+X_{3}+U_{1},2)$
$\displaystyle N_{2}$ $\displaystyle\sim\mathcal{N}(1+X_{2}+X_{3}+U_{2},2)$
$\displaystyle N_{3}$ $\displaystyle\sim\mathcal{N}(0,2).$

In broad terms, $N_{1}$ captures $10/11$ of the total variation of unmeasured
confounding, while $N_{2}$ captures the remaining $1/11$. We then refer to
$N_{1}$ as a strong NCO and $N_{2}$ as a weak NCO. $N_{3}$ is simulated as a
noise variable representing an extreme scenario where an NCO is entirely
irrelevant.

Figure 7: Relative RMSE, coverage rate, relative CI length, and average
selected proportion of RWD of the ESVE-TMLE estimator (Dang et al., 2023c)
across scenarios without NCO, with a weak NCO, with a strong NCO and with an
irrelevant NCO. The MSE-minimizing estimator (Oberst et al., 2022) is included
for reference.

Figure 7 compares the performance is different versions of the ESCV-TMLE
estimator: without an NCO (solid line), with a weak NCO (dashed line), with a
strong NCO (double line), and with an irrelevant NCO (dotted line). The MSE-
minimizing estimator (blue solid line) is includes for reference.
Interestingly, the inclusion of NCOs with varying degrees of relevance affects
the estimator’s behaviour differently. With a strong NCO, which captures most
of the unobserved confounding, the performance is optimal. It shows
substantial efficiency gains when RWD is unbiased and maintains a coverage
rate consistently around 95%. In this scenario, the ESCV-TMLE estimator
dominates the MSE-minimizing estimator, as U-comparability is nearly satisfied
the NCO strongly effectively enhances the bias estimation in the trade-off.
However, if the NCO is weak or totally irrelevant, the estimator become overly
ambitious, resulting in higher worst-case RMSE and CIs under-covering. In
these cases, U-compatibility is violated, and the algorithm places undue
confidence in an unreliable NCO.

The rightmost panel of Figure 7 shows the average proportion of the RWD
included, which helps explain the different behaviours across scenarios.
Compared to not using an NCO, the versions with NCOs tend to select more RWD
when the RWD bias is near zero, which is beneficial. However, the curves
return back to zero more slowly than when not using an NCO. For the versions
with a weak or irrelevant NCO, the proportions selected are consistently
higher than when not using an NCO. This occurs because the squared-bias term
in the selector function is underestimated under the illusion of minimal
confounding bias due to the misleading NCO.

These findings suggest that when a candidate NCO is available, it is crucial
to assess its relevance. A weak or irrelevant NCO can be misleading and may
not function as intended. Therefore, using an NCO is not always better than
not using one.

###  4.3 Augmenting only the control arm of the RWD

####  4.3.1 PROCOVA (Schuler et al., 2021)

As described in Section 3.6, the PROCOVA method only uses the control arm of
the RWD to fit a prognostic score model, $m(x)=\mathbb{E}[Y\,|\,A=0,x,S=0]$,
which is then added as a feature in the RCT dataset. Consequently, the
interpretation of the x-axis in Figure 2 and 3 is more nuanced. Since the
PROCOVA method does not use $\hat{\tau}_{o}$, the impact of unobserved
confounding is more on the expectation of the outcomes of the external
controls . If there is no unobserved confounding, we expect
$\mathbb{E}[Y\,|\,A=0,X,S=0]=\mathbb{E}[Y\,|\,A=0,X,S=1]$. The stronger the
confounding, the larger the difference between these expectations.

In Figure 2 the PROCOVA line overlaps with the reference line, with a coverage
rate constantly close to 95% and no power gain, as shown in Figure 3. In this
case, there is no benefit from data fusion. This is because the simulation
setup, as outlined in Section 4.1, assumes a linear outcome. As discussed in
Section 3.6, the additional prognostic feature does not add any extra
information if the outcome is linear, as noted by Schuler et al. (2021).
Additionally, the relative RMSE curve of the PROCOVA estimator is flat,
indicating that the estimator remains unaffected by the difference between
$\mathbb{E}[Y\,|\,A=0,X,S=0]$ and $\mathbb{E}[Y\,|\,A=0,X,S=1]$.

To understand how the PROCOVA estimator performs with non-linearity in the
outcome, echoing the simulations in Schuler et al. (2021), we introduce
another standard normal variable $X_{4}$ and add a quadratic term $X_{4}^{2}$
to the outcome:

$\displaystyle Y(a)\mid X_{3},U_{1},U_{2}$
$\displaystyle\sim\mathcal{N}\\!\left(0.5+0.2a+0.1a\times X_{3}+X_{3}+\alpha
X_{4}^{2}+U_{1}+U_{2},\,1\right).$

Figure 8: Relative RMSE, coverage rate and relative CI length of PROCOVA
estimator (Schuler et al., 2021) for outcomes with various degrees of
curvature $\alpha=\\{0,0.5,1,2\\}$.

The coefficient $\alpha$ controls the degree of curvature in the outcomes.
When $\alpha=0$, the outcome is linear, as in the main simulation. Figure 8
shows that as $\alpha$ increases, the efficiency gain increases without
compromising the coverage rate. Note that the rightmost panel of Figure 8
plots power on a relative scale, as the ratio against the power of
$\hat{\tau}_{r}$. This is because, as $\alpha$ varies, the variance of $Y$
changes, making it sensible to compare on a relative rather than an absolute
scale. In this specific simulation, when $\alpha=2$, on average, the PROCOVA
estimator shortens the CIs by approximately 12%.

However, the above simulation assumes that the prognostic score is perfect,
i.e. $m(X)=\mathbb{E}[Y(0)|X,S=1]$ using the notations in Schuler et al.
(2021). The authors relaxed this assumption in their Theorem 2, replacing it
with the condition $\left|m(X)-\mu_{0}(X)\right|\xrightarrow{L_{2}}0$, where
$\mu_{0}(X)\equiv\mathbb{E}[Y(0)|X,S=1]$. To test the sensitivity to the
perfect prognostic score assumption, we distinguished the curvature
coefficient $\alpha$ for the RCT and RWD models by $\alpha_{r}$ and
$\alpha_{o}$. Note that in this setting, the condition
$\left|m(X)-\mu_{0}(X)\right|\xrightarrow{L_{2}}0$ still holds. Table 2 gives
the relative RMSE in different ($\alpha_{o}$,$\alpha_{r}$) scenarios. We
observe that the most efficiency gain is generally achieved when
$\alpha_{r}=\alpha_{o}$, and the benefit diminishes when there is a
discrepancy. Interestingly, when the RCT outcome is quadratic w.r.t. $X_{4}$
but the RWD is linear without $X_{4}$, the RMSE of the PROCOVA estimator is
even higher than $\hat{\tau}_{r}$. In this case, the prognostic score
introduces more noise than it explains in $Y$’s residual variation.
Essentially, in a most simplified multivariate regression scenario, the
variance of the regression coefficient of the treatment $A$ can be expressed
as:

$\displaystyle\operatorname{Var}({\hat{\beta}_{A}}\,|\,\textbf{X})$
$\displaystyle=\frac{\sigma_{Y|X_{1}\dots
X_{d}}^{2}}{n\times\sigma^{2}_{A|X_{1},\dots X_{d}}}$
$\displaystyle=\frac{\sigma_{Y|X_{1}\dots
X_{d}}^{2}}{n_{r}\times(1-R^{2}_{A|X_{1},\dots X_{d}})\sigma_{A}^{2}},$

where $\sigma_{Y|X_{1}\dots X_{d}}^{2}$ is the residual variance of $Y$ after
regressing on $X_{1}\dots X_{d}$ and $R^{2}_{A|X_{1},\dots X_{d}}$ denotes the
$R^{2}$ of regressing $A$ on all covariates. Although in an RCT, the treatment
is randomized and $R^{2}_{A|X_{1},\dots X_{d}}$ is expected to be zero in
theory, this may not hold true in finite samples, especially in small RCTs.
Adding an extra variable to the regression, regardless of its relevance,
$\sigma_{Y|X_{1}\dots X_{d}}^{2}$ will either decrease or stay the same, and
$R^{2}_{A|X_{1},\dots X_{d}}$ will either increase or remain the same. If the
additional variable fails to reduce the residual variance of $Y$, but instead
increases $R^{2}_{A|X_{1},\dots X_{d}}$, then variance of $\hat{\beta}_{A}$
will be higher as a result. This explains why, when $\alpha_{r}>0$ but
$\alpha_{o}=0$, the relative RMSE exceeds 100%.

Relative RMSE | $\alpha_{r}=0$ | $\alpha_{r}=0.5$ | $\alpha_{r}=1$ | $\alpha_{r}=2$
---|---|---|---|---
$\alpha_{o}=0$ | 100% | 112% | 130% | 177%
$\alpha_{o}=0.5$ | 100% | 98% | 96% | 89%
$\alpha_{o}=1$ | 100% | 100% | 96% | 88%
$\alpha_{o}=2$ | 100% | 100% | 92% | 84%
Table 2: Relative RMSE of PROCOVA estimator when the prognostic score is
imperfect.

The PROCOVA method proposes a way to calculate the expected reduction in
required RCT sample size prior to conducting the trial. The calculation is
based on the expected correlation between the prognostic score and the
residuals of $Y$ in the RCT. The estimation of the correlation relies entirely
on the RWD and is sensitive to the assumption that the prognostic score is
perfect. For example, when $\alpha_{o}=2$ the estimated correlation is
approximately 0.22, and one would conclude a roughly 4% reduction in the
required RCT sample size. However, if $\alpha_{r}=0.5$, the actual correlation
between the prognostic score and the RCT outcome residuals is close to 0,
indicating that the utility of adding a prognostic score is overestimated. If
we actually have reduced the RCT size by 4% as suggested by the the procedure,
we may end with inadequate power. This example highlights the limitation of
the proposed sample size calculation procedure.

In the above discussion, we explored the sensitivity of the PROCOVA estimator
to the degree of outcome non-linearity and the assumption of perfect
prognostic score. Nonetheless, it maintains the nominal coverage rate
regardless of the magnitude of confounding, safeguarding against Type I error,
which is crucial for informing high-stakes decisions. The ability to calculate
the expected sample size reduction also makes it a practical choice among data
fusion methods.

####  4.3.2 Power likelihood (Lin et al., 2024)

Figure 9 presents the performance of the power likelihood estimator. The
pattern is consistent with its application in integrating both arms: when
confounding bias is small, the estimator effectively reduces variance while
maintaining a coverage rate close to the nominal. However, as RWD bias
intensifies, the estimator’s RMSE exceeds that of $\hat{\tau}_{r}$ due to the
excessive bias introduced by including RWD, resulting in CI undercoverage. As
confounding bias continues to increase, making the RWD increasingly dissimilar
to the RCT, the selected $\hat{\eta}$ to optimize ELPD is reduced to zero.
This avoids extreme bias in the estimate, causing the RMSE to return to the
reference line.

Interestingly, the reduction in variance compared to $\hat{\tau}_{r}$ is more
pronounced as $\alpha$ grows. This might be due to the fact that with stronger
non-linearity or more complex outcome models, the correctly-specified
parametric outcome model outperforms non-parametric machine learning methods
within the AIPW estimator, especially with small sample size. However,
potential model misspecification in practice must be carefully considered.

Figure 9: Relative RMSE, coverage rate and relative CI length of the power
likelihood estimator (Lin et al., 2024) for outcomes with various degrees of
curvature $\alpha=\\{0,0.5,1,2\\}$.

####  4.3.3 ESCV-TMLE (Dang et al., 2023c)

Figure 10 compares different versions of the ESCV-TMLE estimator. The findings
mostly align with those in Section 4.2.4. Including a weak or irrelevant NCO
can expose the estimation to higher risk of bias.

Notably, compared to PROCOVA and power likelihood estimators, the benefits of
the ESCV-TMLE (no NCO) estimator are limited, especially as the outcome
curvature intensifies. In the scenario where $\alpha=2$, there is no reduction
in the RMSE, and in fact, the ESCV-TMLE estimator performs worse than the RCT-
only estimator, exhibiting a higher RMSE. van der Laan et al. (2024) reached a
similar conclusion in their simulations. In their scenario (c) and (d), where
complex non-linear terms were introduced into the generative model, the ESCV-
TMLE failed to achieve any efficiency gain. They suspected that this might be
due to the ESCV-TMLE not accurately estimating the bias structure or its
cross-validation method not fully leveraging the data.

Figure 10: Relative RMSE, coverage rate and relative CI length of different
versions of the ESCV-TMLE estimator (Dang et al., 2023c): without NCO, with a
strong NCO, with a weak NCO and with an irrelevant NCO for outcomes with
various degrees of curvature $\alpha=\\{0,0.5,1,2\\}$.

##  5 Summary

In this paper, we reviewed recent data fusion methods aimed at improving
efficiency.

We first identified key challenges in data fusion, particularly the potential
violations of comparability (Assumption 4) and exchangeability of the RWD
(Assumption 5). These violations result in the ATE estimated from the RWD
biased from the target ATE. When evaluating a data fusion method, it is
crucial to examine these two assumptions and assess their plausibility in the
specific analysis. Additionally, it is important to consider how the method
adapts to such biases.

Secondly, while methodological developments typically provide a theoretical
proof of asymptotic consistency and efficiency gains, the simulations in
Section 4 reveal important common patterns in their finite-sample behaviour.

Thirdly, we discussed the trade-offs between bias and efficiency. There is a
notable trade-off along the scale of bias due to the violation of
comparability and exchangeability assumptions. Potential efficiency gains when
RWD bias is minimal come with the risk of excessive bias in the fusion
estimate within the ‘danger zone’ (approximately 2–5 standard deviations of
$\hat{\tau}_{r}$). This trade-off impacts the coverage rate and relative CI
length, where methods capable of higher efficiency gains also risk producing
estimates that are more biased.

Choosing the appropriate method is critical. Different methods inherently
follow different trade-offs, and no single method dominates all others. What
we can say with certainty is that while bias correction methods (Kallus et
al., 2018; Yang et al., 2020) enhance the efficiency of HTE function
estimation, these efficiency gains vanish when the HTE function is
marginalized over the RCT sample to estimate the ATE. Additionally, the
anchored thresholding method (Chen et al., 2021) is not recommended due to its
lack of robustness against bias from violations of Assumptions 4 and 5. As
bias increases, this method fails to revert to using the RCT alone.

Lastly, we explored the conservatism versus ambition in fusion methods. A more
conservative fusion method is not necessarily better than an ambitious one as
conservative methods often yield less efficiency gain even in optimal
scenarios without confounding bias. If a method does not have the potential to
deliver sufficient power gains in such scenarios, the data fusion analysis
will not fulfil its intended purpose.

Table 1 outlines practical considerations when choosing a data fusion method.
Figure 11 provides a guide for selecting a method based on specific research
needs. A key initial consideration is the availability of treated subjects in
RWD and whether to include them in the fusion. This is crucial when new
treatments or policies are not yet available in the real world. If no unit
received treatment in the RWD, the methods to consider are power likelihood
(Lin et al., 2024), ESCV-TMLE (Dang et al., 2023c) and PROCOVA (Schuler et
al., 2021). Notably, PROCOVA provides Type 1 error protection and sample size
calculation, although its power gains are often limited.

If including both arms from the RWD, the type of outcome dictates the next
step. For continuous outcomes where the typical estimand is the difference in
outcomes, options include power likelihood (Lin et al., 2024), ESCV-TMLE (Dang
et al., 2023c), elastic integration Yang et al. (2023), shrinkage methods
(Green et al., 2005; Rosenman et al., 2023a), and MSE-Minimizing approach
(Oberst et al., 2022). The choice among these should consider the trade-off
between efficiency and bias, as well as implementation complexity.

For binary outcomes, the choice of method depends on the estimand. While the
current implementation of ESCV-TMLE (Dang et al., 2023c) is based on the risk
difference, other methods including the power likelihood (Lin et al., 2024),
the two shrinkage methods (Green et al., 2005; Rosenman et al., 2023a), and
the MSE-Minimizing approach (Oberst et al., 2022) can adapt to alternative
estimands like relative risk.

Figure 11: Decision framework for selecting appropriate data fusion
methodologies based on analysis type, outcome type, and estimand in clinical
research.

Through this paper, our goal is to provide researchers with a reference point
and guidance for navigating various data fusion methods, enabling informed
choices suited to specific research questions. For statisticians interested in
further advancing this topic, we hope this paper presents a comprehensive
overview of the current landscape, highlighting the challenges and
opportunities present. In a subsequent paper, we apply the data fusion methods
compared here to a real-world study, augmenting a medical trial for diabetic
medication with US health claims data. This follow-up will detail practical
considerations such as target trial emulation, handling covariate and outcome
missingness, and variable selection.

## Acknowledgments

This work was conducted during Xi Lin’s internship at Novo Nordisk, Denmark,
in June to September 2023. We are grateful to Kajsa Kvist, Sky Qiu, Henrik
Ravn, Kim Clemmensen, Morten Medici and the Joint Initiative of Causal
Inference (JICI) for their valuable input and helpful discussions that
contributed to this research. We also thank Chenyin Gao for his assistance in
implementing the Elastic Integration method.

##  Appendix A Implementation details of methods in simulation studies

In this appendix, we provide a detailed description of the implementation of
the data fusion methods used in our simulation studies presented in Section 4.

###  A.1 Baseline estimators $\hat{\tau}_{r}$ and $\hat{\tau}_{o}$

As outlined in Section 2.4, we use AIPW estimators $\hat{\tau}_{r}$ and
$\hat{\tau}_{r}$ as reference. The AIPW estimator involves a propensity model
and an outcome model. For both, we employ a discrete SuperLearner (Van der
Laan et al., 2007) ensemble with a library consisting of a mean model and
generalized linear models (GLM) (via R package speedglm (Enea et al., 2009)).
Given the simple generative models for $A$ and $Y$, this library should
suffice. For simulations in Section 4.3 where a quadratic term is added, we
also include generalized additive models (Hastie and Tibshirani, 1987) and
random forest (via R package ranger (Wright and Ziegler, 2017)) into the
library to allow for non-linear flexibility. Additionally, we implemented a
bound of $0.025$ on the propensity score to avoid extreme weights.

We use an empirical sandwich estimation of the variance of $\hat{\tau}_{r}$
following:

$\hat{\sigma}_{r}^{2}=\frac{1}{n_{r}^{2}}\sum_{{i\in\mathcal{R}}}\hat{I}_{i}^{2}$

where

$\displaystyle\hat{I}_{i}=$
$\displaystyle\left[\frac{A_{i}\left(Y_{i}-\hat{m}_{1,r}(X_{i})\right)}{\hat{e}_{r}(X_{i})}-\frac{(1-A_{i})\left(Y_{i}-\hat{m}_{0,r}(X_{i})\right)}{1-\hat{e}_{r}(X_{i})}+\left(\hat{m}_{1,r}(X_{i})-\hat{m}_{0,r}(X_{i})\right)\right]-\hat{\tau}_{r}.$

and similarly for $\hat{\sigma}_{o}^{2}$, estimated variance of
$\hat{\tau}_{o}$.

###  A.2 Shrinkage estimator (Green et al., 2005; Rosenman et al., 2023a)

The shrinkage estimators in Green et al. (2005) and Rosenman et al. (2023a)
require the dimension of the vectors, i.e. $K$, to be greater than two and the
elements to be independent. This necessitates stratifying the data. Rosenman
et al. (2023a) discussed the performance of shrinkage estimators with
different numbers of strata based on empirical evidence from simulations.
However, their findings were inconclusive: the relative performance of the
shrinkage estimators varied with different numbers of strata, and no single
estimator consistently outperformed the others. The authors identified this as
a potential topic for future research.

Despite the lack of definitive guidance, Rosenman et al. (2023a) found
$\hat{\tau}_{S2}$ performed well with relatively few strata (2, 3, or 5).
Based on this, we decided to stratify the RCT and RWD, respectively into four
subgroups according to the values of $X_{1}$ and $X_{3}$. As outlined in
Section 4.1, $X_{1}$ and $X_{3}$ are binary variables, and the ATE estimates
within these subgroups are independent, meeting the required assumptions.

To construct the elements in the shrinkage estimators $\hat{\bm{\tau}}_{S1}$
and $\hat{\bm{\tau}}_{S2}$, we use the AIPW algorithm to get the subgroups ATE
estimates of the RCT, forming $\hat{\bm{\tau}}_{\bm{r}}$ vector. We apply the
same procedure to the RWD. To estimate matrix $\bm{\Sigma}_{r}$, we use the
sandwich estimators of variance of $\tau_{rk}$’s. For $\hat{\bm{\tau}}_{S1}$,
there is a hyper-parameter $a$ and we set $a=K-2$, which is a customary
selection for James-Stein style estimators.

As mentioned in Section 4.2.3, Green et al. (2005) and Rosenman et al. (2023a)
both proposed EB CIs for the elements of $\hat{\bm{\tau}}_{S1}$ and
$\hat{\bm{\tau}}_{S2}$. However, since we are focused on the ATE, these
element-wise CIs are not helpful. Therefore we used non-parametric bootstrap
method, separately resampling the RCT and RWD data respectively, and obtained
the quantiles from 200 iterations.

###  A.3 Bias-correction methods

####  A.3.1 Experiment grounding (Kallus et al., 2018)

We follow the experimental grounding approach as described in Algorithm 1 in
Kallus et al. (2018). Specifically, we implement the following steps:

  1. Step 1:

We fit two linear outcome regression models on the observational data:
$m_{1}(X_{1},X_{2},X_{3})=\mathbb{E}[Y\,|\,X_{1},X_{2},X_{3},A=1,S=0]$ and
$m_{0}(X_{1},X_{2},xX3)=\mathbb{E}[Y\,|\,X_{1},X_{2},X_{3},A=0,S=0]$ , and
calculate their difference as an estimate of CATE estimate
$\hat{\omega}(\bm{z},\bm{c})=\hat{m}_{1}(X_{1},X_{2},X_{3})-\hat{m}_{0}(X_{1},X_{2},X_{3})$.

  2. Step 2:

We learn the bias correction function, assuming linearity,
$\xi(\bm{x})=\beta^{T}\bm{x}$ where $\bm{x}$ is the vector of
$(x_{1},x_{2},X_{3})^{T}$ and

$\hat{\beta}=\underset{\beta}{\arg\min}\sum_{i\in\mathcal{R}}\left(q_{i}Y_{i}-\hat{\omega}\left(\bm{X}_{i}\right)-\beta^{T}\bm{X}_{i}\right)^{2},$

where $q_{i}=2\text{ if }T=1$ and $q_{i}=-2\text{ if }T=0$. That is, we fit
$\eta(\bm{z},\bm{c})$ through least square regression on
$q_{i}Y_{i}-\hat{\omega}(\bm{V}_{i})$ using the randomized data.

  3. Step 3:

We then calculate the ATE estimates as

$\hat{\tau}_{XG}=\frac{1}{n_{r}}\sum_{i\in\mathcal{R}}\hat{\omega}(\bm{X}_{i})+\hat{\beta}^{T}\bm{X}_{i}.$

For inference, we create confidence intervals using non-parametric bootstrap
approach by resampling the RCT and RWD data respectively, obtaining the
quantiles from 200 iterations.

####  A.3.2 Confounding function (Yang et al., 2020)

We implement the confounding function estimator using the R package
IntegrativeHTEcf (Yang, 2020). This package assumes a linear form for both the
confounding and HTE functions, limiting its application to binary outcomes.
Nonetheless, the linearity assumption is appropriate for the simulations
presented. We specify the covariate set as $(X_{1},X_{2},X_{3})$ and the HTE
covariate set to be $X_{3}$, aligning with the correct data generation model.

The main function, IntHTEcf, outputs the estimated linear coefficients of the
HTE and confounding functions, with variance estimated via an internal non-
parametric bootstrap procedure. Since we are interested in the ATE, instead of
the uncertainty around the HTE coefficients, we conduct our own bootstrap
procedure outside the function to estimate the uncertainty of the ATE
estimates.

###  A.4 Adaptive power likelihood (Lin et al., 2024)

The implementation of this method draws upon R packages causl (Evans, 2021)
and ManyData (Lin, 2023). We follow these steps:

First, following the frugal parameterization, we parameterize the causal model
as:

$\displaystyle X_{1}$ $\displaystyle\sim\mathcal{N}(\beta_{1},\phi_{1})$
$\displaystyle X_{2}$
$\displaystyle\sim\operatorname{Bernoulli}(\operatorname{expit(\beta_{2})})$
$\displaystyle Y\mid do(A),X_{3}$
$\displaystyle\sim\mathcal{N}(\beta_{4}+\beta_{5}A+\beta_{6}X_{3}A+\beta_{7}X_{3},\phi_{4})$

and the dependence between $X_{1}$,$X_{2}$,$Y$ follow a multi-variate Gaussian
copula with constant correlation coefficients $\rho_{1}$, $\rho_{2}$ and
$\rho_{3}$. We denote the vector consisting of $\beta$, $\phi$ and $\rho$
parameters as $\bm{\theta}$ . We do not need to assume parametric models for
$A$ and $X_{3}$ as we are interested in the distribution conditional on
treatment and the effect modifier.

Using the fitCausl function in causl, we find $\hat{\bm{\theta}}_{r}$ and
$\hat{\bm{\theta}}_{o}$, which are the MLEs in the RCT and RWD. Next, we
perform a grid search for the optimal learning rate $\eta\in[0,1]$. For each
value of $\eta$ on the grid, we approximate the $\eta$-indexed posterior of
$\bm{\theta}$ as described in Section 4.3 of Lin et al. (2024). We then draw
random samples from this posterior distribution, $\bm{\theta}^{(i)}$, for
$i=0,1,\dots,2000$. For each sample $\bm{\theta}^{(i)}$, we approximate the
ELPD using Leave-One-Out cross validation, implemented via the loo function in
the R package loo (Vehtari et al., 2020). Looping through the grid of $\eta$,
we select the optimal value $\hat{\eta}$ that maximizes the approximated ELPD.

We then obtain the distribution of the ATE estimate based on
$\bm{\theta}^{(i)}$ samples:

$\hat{\tau}_{PL}^{(i)}=\frac{1}{n_{r}}\sum_{j\in\mathcal{R}}\hat{\beta}_{5}^{(i)}+\hat{\beta}_{6}^{(i)}\times
X_{3j},$

for $i=1,2,\dots,2000$. Finally, we derive the estimator
$\hat{\tau}_{PL}=\frac{1}{2000}\sum_{i=1}^{2000}\hat{\tau}_{PL}^{(i)}$, as
well as the credible interval based on the quantiles.

###  A.5 Test-then-pool methods

####  A.5.1 Elastic integration (Yang et al., 2023)

We implement the elastic integration estimator using the R package
ElasticIntegrative (Yang, 2020). This package is use-friendly and
straightforward to use. The main function is elasticHTE and we specify
$(X_{1},X_{2},X_{3})$ as the covariate set for the propensity score model and
$X_{3}$ for the HTE function, aligning with the data generative model.
elasticHTE implemented flexible modeling through SuperLearner (Van der Laan et
al., 2007) and we specify the same libraries used for the AIPW estimation.

For inference, as the estimator is not regular, the variance estimated through
non-parametric bootstrap will be biased; this is formalized with a proof in
the supplementary materials of their paper. Therefore, the authors implemented
a permutation-based bias estimation within the main function. We use the
variance output from the function, although we modified the function slightly
to also output the covariance. As we need to marginalize it over the covariate
distribution of the RCT, we need to account for the correlation between
coefficient estimates.

####  A.5.2 ESCV-TMLE Dang et al. (2023c)

We implement the ESCV-TMLE estimator using the R package EScvtmle (Dang et
al., 2023b). This package is use-friendly and straightforward to use. The main
function ES.cvtmle provides flexibility for combining either both arms or only
the control arm of the RWD, the inclusion of a NCO, as well as different
distribution families for the outcome and the NCO. We specify these settings
accordingly for the simulations in Section 4.2 and 4.3. For the integrated
SuperLearner (Van der Laan et al., 2007) ensemble, we choose the discrete
version with the same libraries used for the AIPW estimation. We use five
folds for cross-validation.

We use the variance output from ES.cvtmle to construct 95% confidence
intervals.

###  A.6 Weighted combination methods

####  A.6.1 MSE-minimizing estimator (Oberst et al., 2022)

We implement the MSE-minimizing estimator following

$\displaystyle\hat{\tau}_{\hat{\lambda}}=\hat{\lambda}\hat{\tau}_{o}+(1-\hat{\lambda})\hat{\tau}_{r}$
where
$\displaystyle\hat{\lambda}=\frac{\hat{\sigma}_{r}^{2}}{\left(\hat{\tau}_{r}-\hat{\tau}_{o}\right)^{2}+\hat{\sigma}_{r}^{2}+\hat{\sigma}_{o}^{2}}.$

As $\hat{\tau}_{e}$ and $\hat{\tau}_{o}$ are independent, the covariance term
in the original formula of $\hat{\lambda}$ drops out. For inference, we
generate confidence intervals by applying a non-parametric bootstrap method,
separately resampling the RCT and RWD data and extracting quantiles from 200
iterations.

####  A.6.2 Anchored-thresholding estimator (Chen et al., 2021)

We implement the anchored-thresholding estimator by first calculating the
estimated bias $\hat{\delta}_{\gamma}$ through soft-shresholding:

$\hat{\delta}_{\gamma}=\begin{cases}\operatorname{sgn}(\hat{\delta})\left(|\hat{\delta}|-\gamma\cdot\sqrt{\hat{\sigma}_{r}^{2}+\hat{\sigma}_{o}^{2}}\right)&\text{
if
}|\hat{\delta}|\geq\gamma\cdot\sqrt{\hat{\sigma}_{r}^{2}+\hat{\sigma}_{o}^{2}}\\\
0&\text{ otherwise. }\end{cases}$

.

As discussed in Section 4.2.3, hyper-parameter $\gamma$ is important in
controlling the combination and the larger $\gamma$ is, the higher the
potential bias is. Chen et al. (2021) state that $\gamma$ should be of
asymptotic order $\sqrt{\log\left(\min(n_{r},n_{o}))\right)}$, so we set
$\gamma=\sqrt{\log 300}$. We then plug $\hat{\delta}_{\gamma}$ into the
estimator:

$\hat{\tau}_{AT}=(1-\hat{\omega})\hat{\tau}_{r}+\hat{\omega}(\hat{\tau}_{o}+\hat{\delta}_{\gamma}),$

where
$\hat{\omega}=\hat{\sigma}_{r}^{2}/(\hat{\sigma}_{r}^{2}+\hat{\sigma}_{o}^{2})$.
We use a non-parametric bootstrap to obtain the 95% confidence interval.

###  A.7 Prognostic adjustment (Schuler et al., 2021)

The implementation of the PROCOVA approach is relatively straight forward. We
first empirically centre the covariates as well as the treatment indicator.
Then employing the SuperLearner ensemble, we fit a prognostic score model and
make predictions on the RCT, creating prognostic scores $\hat{m}(X)$.Then, we
estimate the treatment effect $\hat{\beta}_{A}$ using a linear regression
adjusted for the empirically centered covariates, prognostic score, and their
interactions with the treatment.

We use the sandwich estimator of variance of the $\hat{\beta}_{A}$ to
construct the 95% confidence interval.

## References

  * Abadie and Imbens (2002) A. Abadie and G. W. Imbens.  Simple and bias-corrected matching estimators for average treatment effects.  _NBER Working Paper_ , (no.283), 2002.
  * Athey et al. (2020) S. Athey, R. Chetty, and G. Imbens.  Combining experimental and observational data to estimate treatment effects on long term outcomes.  _arXiv preprint 2006.09676_ , 6 2020.
  * Bareinboim and Pearl (2016) E. Bareinboim and J. Pearl.  Causal inference and the data-fusion problem.  _Proceedings of the National Academy of Sciences_ , 113(27):7345–7352, 2016.
  * Bonapersona et al. (2021) V. Bonapersona, H. Hoijtink, RELACS Consortium, R. Sarabdjitsingh, and M. Joëls.  Increasing the statistical power of animal experiments with historical control data.  _Nature Neuroscience_ , 24(4):470–477, 2021.
  * Brantner et al. (2023) C. L. Brantner, T.-H. Chang, T. Q. Nguyen, H. Hong, L. Di Stefano, and E. A. Stuart.  Methods for integrating trials and non-experimental data to examine treatment effect heterogeneity.  _arXiv preprint arXiv:2302.13428_ , 2023.
  * Chen et al. (2021) S. Chen, B. Zhang, and T. Ye.  Minimax rates and adaptivity in combining experimental and observational data.  9 2021.  URL http://arxiv.org/abs/2109.10522.
  * Cheng and Cai (2021) D. Cheng and T. Cai.  Adaptive combination of randomized and observational data.  _arXiv preprint arXiv:2111.15012_ , 2021.
  * Colnet et al. (2020) B. Colnet, I. Mayer, G. Chen, A. Dieng, R. Li, G. Varoquaux, J.-P. Vert, J. Josse, and S. Yang.  Causal inference methods for combining randomized trials and observational studies: a review.  11 2020.  URL http://arxiv.org/abs/2011.08047.
  * Colnet et al. (2022) B. Colnet, J. Josse, G. Varoquaux, and E. Scornet.  Reweighting the rct for generalization: finite sample analysis and variable selection.  _arXiv preprint arXiv:2208.07614_ , 2022.
  * Dang et al. (2023a) L. E. Dang, E. Fong, J. M. Tarp, K. K. B. Clemmensen, H. Ravn, K. Kvist, J. B. Buse, M. van der Laan, and M. Petersen.  Case study of semaglutide and cardiovascular outcomes: An application of the causal roadmap to a hybrid design for augmenting an rct control arm with real-world data.  _Journal of Clinical and Translational Science_ , 7(1):e231, 2023a.
  * Dang et al. (2023b) L. E. Dang, M. Peterson, and M. van der Laan.  EScvtmle.  https://github.com/Lauren-EylerDang/EScvtmle, 2023b.
  * Dang et al. (2023c) L. E. Dang, J. M. Tarp, T. J. Abrahamsen, K. Kvist, J. B. Buse, M. Petersen, and M. van der Laan.  A cross-validated targeted maximum likelihood estimator for data-adaptive experiment selection applied to the augmentation of RCT control arms with external data.  _arXiv preprint arXiv:2210.05802_ , 2023c.
  * Dawid (2021) A. P. Dawid.  Decision-theoretic foundations for statistical causality.  _Journal of Causal Inference_ , 9(1):39–77, 2021.
  * Duan et al. (2006) Y. Duan, K. Ye, and E. P. Smith.  Evaluating water quality using power priors to incorporate historical information.  _Environmetrics: The Official Journal of the International Environmetrics Society_ , 17(1):95–106, 2006.
  * Enea et al. (2009) M. Enea et al.  Fitting linear models and generalized linear models with large data sets in R.  _Statistical Methods for the Analysis of Large Datasets: book of short papers_ , pages 411–414, 2009.
  * European Medicines Agency (2022) European Medicines Agency.  Qualification opinion for prognostic covariate adjustment (PROCOVA™).  https://www.ema.europa.eu/en/documents/regulatory-procedural-guideline/qualification-opinion-prognostic-covariate-adjustment-procovatm-en.pdf, 2022.  Accessed: 2023-03-13.
  * Evans and Didelez (2024) R. Evans and V. Didelez.  Parameterizing and simulating from causal models.  _Journal of the Royal Statistical Society Series B: Statistical Methodology (in press)_ , 2024.
  * Evans (2021) R. J. Evans.  _causl: Methods for Specifying, Simulating from and Fitting Causal Models_ , 2021.  R package version 0.4.2.9001.
  * Fan et al. (2017) J. Fan, H. Liu, Y. Ning, and H. Zou.  High dimensional semiparametric latent graphical model for mixed data.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 79(2):405–421, 2017.
  * Green and Strawderman (1991) E. J. Green and W. E. Strawderman.  A James-Stein type estimator for combining unbiased and possibly biased estimators.  _Journal of the American Statistical Association_ , 86:1001–1006, 1991.
  * Green et al. (2005) E. J. Green, W. E. Strawderman, R. L. Amateis, and G. A. Reams.  Improved estimation for multiple means with heterogeneous variances.  _Forest Science_ , 51(1):1–6, 2005.
  * Hastie and Tibshirani (1987) T. Hastie and R. Tibshirani.  Generalized additive models: some applications.  _Journal of the American Statistical Association_ , 82(398):371–386, 1987.
  * Hobbs et al. (2011) B. P. Hobbs, B. P. Carlin, S. J. Mandrekar, and D. J. Sargent.  Hierarchical commensurate and power prior models for adaptive incorporation of historical information in clinical trials.  _Biometrics_ , 67(3):1047–1056, 2011.
  * Hobbs et al. (2012) B. P. Hobbs, D. J. Sargent, and B. P. Carlin.  Commensurate priors for incorporating historical information in clinical trials using general and generalized linear models.  _Bayesian Analysis (Online)_ , 7(3):639, 2012.
  * Ibrahim and Chen (2000) J. G. Ibrahim and M.-H. Chen.  Power prior distributions for regression models.  _Statistical Science_ , pages 46–60, 2000.
  * Imbens (2004) G. W. Imbens.  Nonparametric estimation of average treatment effects under exogeneity: A review.  _Review of Economics and statistics_ , 86(1):4–29, 2004.
  * James and Stein (1961) W. James and C. Stein.  Estimation with quadratic loss.  _Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability_ , 1:361–380, 1961.
  * Kallus et al. (2018) N. Kallus, A. M. Puli, and U. Shalit.  Removing hidden confounding by experimental grounding.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Li and Luedtke (2023) S. Li and A. Luedtke.  Efficient estimation under data fusion.  _Biometrika_ , 110(4):1041–1054, 2023.
  * Li et al. (2020) W. Li, F. Liu, and D. Snavely.  Revisit of test-then-pool methods and some practical considerations.  _Pharmaceutical statistics_ , 19(5):498–517, 2020.
  * Lin (2023) X. Lin.  Manydata: Combining experimental and observational data through a power likelihood.  https://github.com/XiLinStats/ManyData, 2023.
  * Lin et al. (2024) X. Lin, J. M. Tarp, and R. J. Evans.  Combine experimental and observational data through a power likelihood.  _arXiv preprint arXiv:2304.02339_ , 2024.
  * Merck (2017) Merck.  European commission approves Bavencio (avelumab) for metastatic merkel cell carcinoma.  https://www.merckgroup.com/en/news/bavencio-ec-approval-2017-09-20.html, 2017.  Accessed: 2024-06-24.
  * Morris (1983) C. N. Morris.  Parametric empirical bayes inference: theory and applications.  _Journal of the American statistical Association_ , 78(381):47–55, 1983.
  * Morucci et al. (2023) M. Morucci, V. Orlandi, H. Parikh, S. Roy, C. Rudin, and A. Volfovsky.  A double machine learning approach to combining experimental and observational data.  _arXiv preprint arXiv:2307.01449_ , 2023.
  * Neuenschwander et al. (2009) B. Neuenschwander, M. Branson, and D. J. Spiegelhalter.  A note on the power prior.  _Statistics in medicine_ , 28(28):3562–3566, 2009.
  * Oberst et al. (2022) M. Oberst, A. D’Amour, M. Chen, Y. Wang, D. Sontag, and S. Yadlowsky.  Bias-robust integration of observational and experimental estimators.  _arXiv preprint arXiv:2205.10467_ , 2022.
  * Pearl (2009) J. Pearl.  _Causality_.  Cambridge University Press, 2009.
  * Richardson and Robins (2013) T. S. Richardson and J. M. Robins.  Single world intervention graphs (SWIGs): A unification of the counterfactual and graphical approaches to causality.  _Center for the Statistics and the Social Sciences, University of Washington Series. Working Paper_ , 128(30):2013, 2013.
  * Robbins et al. (2022) M. W. Robbins, S. Bauhoff, and L. Burgette.  Data fusion for predicting long-term program impacts.  _arXiv preprint arXiv:2205.01804_ , 2022.
  * Robins (2000) J. M. Robins.  Robust estimation in sequentially ignorable missing data and causal inference models.  In _Proceedings of the American Statistical Association_ , volume 1999, pages 6–10. Indianapolis, IN, 2000.
  * Robins et al. (1994) J. M. Robins, A. Rotnitzky, and L. P. Zhao.  Estimation of regression coefficients when some regressors are not always observed.  _Journal of the American statistical Association_ , 89(427):846–866, 1994.
  * Rosenman et al. (2023a) E. T. Rosenman, G. Basse, A. B. Owen, and M. Baiocchi.  Combining observational and experimental datasets using shrinkage estimators.  _Biometrics_ , 79(4):2961–2973, 2023a.
  * Rosenman et al. (2023b) E. T. Rosenman, F. Dominici, and L. Miratrix.  Empirical Bayes double shrinkage for combining biased and unbiased causal estimates.  _arXiv preprint arXiv:2309.06727_ , 2023b.
  * Rotnitzky et al. (1998) A. Rotnitzky, J. M. Robins, and D. O. Scharfstein.  Semiparametric regression for repeated outcomes with nonignorable nonresponse.  _Journal of the american statistical association_ , 93(444):1321–1339, 1998.
  * Rubin (1974) D. B. Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of Educational Psychology_ , 66(5):688, 1974.
  * Rubin (1978) D. B. Rubin.  Bayesian inference for causal effects: The role of randomization.  _The Annals of statistics_ , pages 34–58, 1978.
  * Schmidli et al. (2014) H. Schmidli, S. Gsteiger, S. Roychoudhury, A. O’Hagan, D. Spiegelhalter, and B. Neuenschwander.  Robust meta-analytic-predictive priors in clinical trials with historical control information.  _Biometrics_ , 70(4):1023–1032, 2014.
  * Schuler et al. (2021) A. Schuler, D. Walsh, D. Hall, J. Walsh, C. Fisher, C. P. for Alzheimer’s Disease, A. D. N. Initiative, and A. D. C. Study.  Increasing the efficiency of randomized trial estimates via linear adjustment for a prognostic score.  _The International Journal of Biostatistics_ , 18(2):329–356, 2021.
  * Shi et al. (2023) X. Shi, Z. Pan, and W. Miao.  Data integration in causal inference.  _Wiley Interdisciplinary Reviews: Computational Statistics_ , 15(1):e1581, 2023.
  * Spirtes et al. (2000) P. Spirtes, C. N. Glymour, R. Scheines, and D. Heckerman.  _Causation, prediction, and search_.  MIT press, 2000.
  * Stein (1956) C. Stein.  Inadmissibility of the usual estimator for the mean of a multivariate distribution.  _Proceedings of the Third Berkeley Symposium_ , pages 197–206, 1956.
  * Tarima et al. (2020) S. Tarima, B. Tuyishimire, R. Sparapani, L. Rein, and J. Meurer.  Estimation combining unbiased and possibly biased estimators.  _Journal of Statistical Theory and Practice_ , 14, 6 2020.  ISSN 15598616.  doi: 10.1007/s42519-020-0083-2.
  * U.S. Food and Drug Administration (2018) U.S. Food and Drug Administration.  FDA approves new treatment for a rare genetic disorder, Fabry disease.  https://www.fda.gov/news-events/press-announcements/fda-approves-new-treatment-rare-genetic-disorder-fabry-disease, 2018.  Accessed: 2024-06-24.
  * van der Laan and Rubin (2006) M. J. van der Laan and D. Rubin.  Targeted maximum likelihood learning.  _The international Journal of Biostatistics_ , 2(1), 2006.
  * Van der Laan et al. (2007) M. J. Van der Laan, E. C. Polley, and A. E. Hubbard.  ”super learner”.  _Statistical applications in genetics and molecular biology_ , 6(1), 2007.
  * van der Laan et al. (2024) M. J. van der Laan, S. Qiu, and L. van der Laan.  Adaptive-TMLE for the average treatment effect based on randomized controlled trial augmented with real-world data.  _arXiv preprint arXiv:2405.07186_ , 2024.
  * Vehtari et al. (2020) A. Vehtari, J. Gabry, M. Magnusson, Y. Yao, P. Bürkner, T. Paananen, and A. Gelman.  loo: Efficient leave-one-out cross-validation and waic for bayesian models. r package version 2.4. 1, 2020.
  * Viele et al. (2014) K. Viele, S. Berry, B. Neuenschwander, B. Amzal, F. Chen, N. Enas, B. Hobbs, J. G. Ibrahim, N. Kinnersley, S. Lindborg, et al.  Use of historical control data for assessing treatment effects in clinical trials.  _Pharmaceutical statistics_ , 13(1):41–54, 2014.
  * Wright and Ziegler (2017) M. N. Wright and A. Ziegler.  ranger: A fast implementation of random forests for high dimensional data in C++ and R.  _Journal of Statistical Software_ , 77(i01), 2017.
  * Yang (2020) S. Yang.  IntegrativeHTEcf.  https://github.com/shuyang-stat/IntegrativeHTEcf, 2020.
  * Yang and Ding (2020) S. Yang and P. Ding.  Combining multiple observational data sources to estimate causal effects.  _Journal of the American Statistical Association_ , 115:1540–1554, 7 2020.  ISSN 1537274X.  doi: 10.1080/01621459.2019.1609973.
  * Yang and Gao (2023) S. Yang and C. Gao.  ElasticIntegrative.  https://github.com/Gaochenyin/ElasticIntegrative, 2023.
  * Yang et al. (2020) S. Yang, D. Zeng, and X. Wang.  Improved inference for heterogeneous treatment effects using real-world data subject to hidden confounding.  _arXiv preprint arXiv:2007.12922_ , 2020.
  * Yang et al. (2023) S. Yang, C. Gao, D. Zeng, and X. Wang.  Elastic integrative analysis of randomised trial and real-world data for treatment heterogeneity estimation.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 85(3):575–596, 2023.