# Many Data: Combine Experimental and Observational Data through a Power
Likelihood

Xi Lin  Department of Statistics, University of Oxford  Robin J. Evans
Department of Statistics, University of Oxford

###### Abstract

Randomized controlled trials are commonly regarded as the gold standard for
causal inference and play a pivotal role in modern evidence-based medicine.
However, the sample sizes they use are often too limited to draw significant
causal conclusions for subgroups that are less prevalent in the population. In
contrast, observational data are becoming increasingly accessible in large
volumes but can be subject to bias as a result of hidden confounding. Given
these complementary features, we propose a power likelihood approach to
augmenting RCTs with observational data for robust estimation of heterogeneous
treatment effects. We provide a data-adaptive procedure for maximizing the
Expected Log Predictive Density (ELPD) to select the influence factor that
best regulates the information from the observational data. We conduct a
simulation study to illustrate the efficacy of our method and its favourable
features compared to existing approaches. Lastly, we apply the proposed method
to data from Tennessee’s Student Teacher Achievement Ratio (STAR) Study to
demonstrate its usefulness and practicality in real-world data analysis.

##  1 Introduction

Experimental data and observational data represent two distinct regimes used
for causal inference. Experimental data are collected through designed
experiments with randomized interventions; a typical example is Randomized
Controlled Trials (RCTs). Randomized intervention isolates the causal effect
of interest from unwanted and potentially unobserved confounding factors as
the experiment protocol ensures, in expectation, balanced control and
treatment groups. As a result, empirical researchers can employ straight-
forward estimation strategies such as regression, inverse propensity weighting
(IPW) (Rosenbaum and Rubin, 1983; Robins et al., 1994; Bang and Robins, 2005;
Cao et al., 2009) and matching (Rubin, 1973; Hirano et al., 2003; Abadie and
Imbens, 2016) to consistently estimate causal effects. Therefore, RCTs are
widely used in clinical settings and regulators, such as the U.S. Food and
Drug Administration (FDA) and the European Medicines Agency (EMA), regard them
as the gold standard to demonstrate the efficacy of a proposed treatment.

Observational data, also sometimes referred to as real-world data, are
passively collected without specially designed intervention. Examples include
electronic health records, bank transactions, and user behaviour data on
online platforms. In this Big Data era, they have become ubiquitous and
readily available in large volumes. Due to the lack of randomization,
observational data are subject to hidden confounding, introducing bias in
causal estimates that does not vanish even with infinite samples (Pearl,
2009).111Specifically, we are concerned about hidden common causes shared by
the treatment and outcome. For example, in the case of flu vaccination,
individuals with more active health-seeking behaviour may be more likely to
receive flu shots. This can introduce bias because they are also less likely
to be hospitalized due to the flu regardless of vaccination status. In this
case, the observed effect of the flu shot on hospitalization rates would be
confounded by the underlying health-seeking behaviour of the individuals.
Hidden confounding is a common problem in observational studies, as there may
be numerous unobserved factors that affect both treatment assignment and the
outcome of interest. For that reason, observational data are traditionally
less favoured compared to RCTs in the context of medicine and healthcare.

Nonetheless, RCTs suffer from limitations too. They can be time-consuming and
expensive to conduct and collect data from, which means that they are limited
by slow turnaround times and inadequate sample sizes. The COVID-19 pandemic
highlighted the value of timely and reliable real-world data, which served as
crucial evidence for determining the safety and effectiveness of drugs and
vaccines (Hansen et al., 2021; Pottegård et al., 2020). For rare diseases, for
example, male breast cancer, it may be infeasible to recruit enough patients
needed to conduct clinical trials with adequate statistical power (Wedam et
al., 2020). In such cases, rigorous analyses of reliable real-world data can
provide valuable insights.

Combining the information from the observational and experimental datasets, a
method of data fusion (Bareinboim and Pearl, 2016), is an obvious avenue to
harvesting the best of both realms. The combination is expected to leverage
the internal validity222First introduced in Campbell (1957), the author
motivated the concept of internal validity with the question “did in fact the
experimental stimulus make some significant difference in this specific
instance?”. Alternatively, we can intuitively understand it as the extent to
which we are confident that a causal relationship established in a study
cannot be explained by other factors.  of experimental data and complement it
with the richness of the observational data. Over recent years, both the FDA
and EMA have shared the vision that real-world data and RCTs should be seen as
complementing each other, and their relative importance depends on the
regulatory issue under consideration (U.S Food and Drug Administration, 2018b;
European Medicines Agency, 2021); these agencies have already granted
regulatory approval to a number of studies and methodologies that rigorously
combine evidence from RCTs and real-world data (Wedam et al., 2020; U.S Food
and Drug Administration, 2018a, 2016; European Medicines Agency, 2022)

Essentially, the combination is inherently a bias-variance trade-off. By
introducing information from the large, possibly-biased real-world data to
limited-scaled, unbiased randomized data, we seek to improve the efficiency of
the causal estimator at the expense of an incremental bias and strike a
balance that results in an overall improvement according to a certain
criterion, for example, mean squared error (MSE). In this paper, we present a
novel power likelihood approach for effectively augmenting RCTs with
observational data to improve the estimation of heterogeneous treatment
effects.

The remainder of this paper is organized as follows: in Section 2, we review
data fusion methods proposed in recent literature. In Section 3 we introduce
the setup of the problem and state the main assumptions, before moving on to
our proposed power likelihood approach in Section 4. In Section 5 we conduct a
simulation study to illustrate the efficacy of our approach and its favourable
features compared to existing proposals. We also apply our method to data from
the STAR study in Section 6 to demonstrate its usefulness and practicality in
real-world data analysis. We conclude with a discussion in the final section.
We defer proofs and additional details of simulations to the appendices. The
code to reproduce the simulations and figures is available on
Github333https://github.com/XiLinStats/ManyData.

##  2 Related Literature

In recent years, combining randomized and observational data for causal
inference has attracted a lot of research attention. Before we review existing
works, we first introduce the definitions of average treatment effect (ATE)
and conditional average treatment effect (CATE).

######  Definition 2.1.

(Average treatment effect) We define the average treatment effect as:

$\tau=\mathbb{E}[Y\,|\,do(T=1)]-\mathbb{E}[Y\,|\,do(T=0)]$

where $\mathbb{E}[Y|do(T=t)]$ represents the expected outcome $Y$ if we set a
dichotomous treatment $T$ to $t\in\\{0,1\\}$ through intervention.

######  Definition 2.2.

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

$\tau(c)=\mathbb{E}[Y|do(T=1),C=c]-\mathbb{E}[Y|do(T=0),C=c]$

where $\mathbb{E}[Y\,|\,do(T=t),C=c]$ represents the expected outcome $Y$ of
an individual with covariate values of $c$, if we set $T$ to $t$ through
intervention. If the treatment effect varies by values of $C$, then we say
that the treatment effect is modified by $C$, or there is treatment effect
heterogeneity across values of $C$.

Owing to the concern about introducing bias into the inference, a number of
initial attempts assumed that the observational data is conditionally
unconfounded in the hope that the combined estimator benefits from the
efficiency gain while remaining unbiased (Hartman et al., 2015; Rosenman et
al., 2022). However, one can never prove that all confounders are measured and
the assumption of unconfoundedness almost never holds in practice (Pearl,
2009). Other works attempted to avoid making this stringent assumption and
replaced it with alternative assumptions on the mechanism that drives how the
unmeasured confounding biases the causal effect (Peysakhovich and Lada, 2016;
Kallus et al., 2018; Athey et al., 2020). Specifically, Kallus et al. (2018)
assumed that the influence of hidden confounding on the outcome had a
parametric structure and can be effectively modelled via “experimental
grounding”, based on the difference between the observational data and the
RCT.

Yang and Ding (2020) formulated the combination problem in a setting where the
experimental and observational samples are drawn from a common distribution,
with all confounders measured in the experimental sample but not all measured
in the observational sample. Their setting essentially assumed that the
unmeasured confounders influence the observational data in the same fashion as
the experimental data. They then proposed an estimator minimizing the
asymptotic variance aiming to improve the efficiency of the estimator from the
randomized sample.

The abovementioned assumptions centring around unconfoundedness are likely to
be reasonable in practice, but they remain unfalsifiable nonetheless. This
motivated several “test-then-pool”-styled approaches (Viele et al., 2014; Yang
et al., 2020). Originally proposed by Viele et al. (2014), such approaches
start with the null hypothesis of equality between the causal estimates from
both datasets and only pool the data together if the null hypothesis is not
rejected. This approach is the most basic form of dynamically borrowing
information from an external, or historical, dataset. Typically, “test-and-
pool” approaches require the researcher to specify the size of the test, say
$\alpha=0.025$, which can be somewhat subjective.Yang et al. (2020) proposed a
data-adaptive approach to dynamically choose the threshold for the test static
based on an estimation of the bias. A major issue with “test-then-pool”-styled
approaches are 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 EMA recently approved a historical borrowing method by Schuler et al.
(2021). They proposed to augment the control arm of RCTs by incorporating a
prognostic score fitted using a historical dataset, which is required to be
drawn from the same population as the trial control arm so that the mean
control outcomes are the same. While they suggested that this assumption can
be weakened, the theoretical efficiency gain still relies on the control
outcome estimated from the historical data being consistent, which may not
hold if the external data comes from a previous trial, or for instance, from
another country.

If an unbiased estimate is not absolutely compulsory and there is some
appetite or tolerance for bias, then there are more options of methods aimed
at optimizing the tradeoff between bias and efficiency gain (Rosenman et al.,
2020; Oberst et al., 2022; Dang et al., 2022). Similar to the setting in
Schuler et al. (2021), Dang et al. (2022) proposed a Targeted Maximum
Likelihood Estimation (TMLE) experiment selector to dynamically incorporate
external data into the trial control arm. They novelly incorporated a bias
estimated using a negative control outcome (NCO) into the Mean Squared Error
to be minimized. As a result, their method relaxed the assumption that the
external data is unconfounded and used the NCO to control the bias in the
control data.

Stein’s Paradox showed that while MLE estimators have the smallest risk among
all unbiased estimators, they are not admissible, meaning that there exist
other estimators with lower MSE regardless of the true value (Stein, 1956;
James and Stein, 1961). One example of such estimator is the classical James-
Stein estimator, which shrinks the MLE towards the zero vector. Green and
Strawderman (1991) and Green et al. (2005) extended this idea and considered
the problem of combining unbiased and possibly biased estimators. Rosenman et
al. (2022) applied this idea specifically to causal inference and developed
several shrinkage-style estimators by minimizing Stein’s Unbiased Risk
Estimate (SURE). Their method requires at least four strata, as well as
independence among stratum estimates, to guarantee a reduction in risk. This
means that its usefulness is limited in the case of ATE or two or three-
category CATE estimation.

Without the dimensionality and independence requirements posed by James-Stein-
styled shrinkage estimators, a straight-forward linear combination strategy
has been popular in recent literature (Tarima et al., 2020; Cheng and Cai,
2021; Oberst et al., 2022). Such estimators take the linear form
$\hat{\theta}_{\lambda}=(1-\lambda)\hat{\theta}_{e}+\lambda\hat{\theta}_{o}$,
where $\hat{\theta}_{e}$ is an unbiased estimator from the RCT and
$\hat{\theta}_{o}$ is a biased estimator from the observational data, with
$\lambda$ being a weight parameter. Specifically, Tarima et al. (2020) and
Oberst et al. (2022) proposed to choose the optimal$\lambda$ by minimizing the
theoretical MSE:

$\lambda^{*}=\frac{\operatorname{Var}(\hat{\theta}_{e})-\operatorname{Cov}(\hat{\theta}_{e},\hat{\theta}_{o})}{\delta^{2}+\operatorname{Var}(\hat{\theta}_{e}-\hat{\theta}_{o})}.$

Obviously, the ‘true’ bias $\delta$ and the variance quantities remain unknown
so both papers considered estimating $\lambda^{*}$ using plug-in estimators.
In particular, they proposed using $(\hat{\theta}_{e}-\hat{\theta}_{o})^{2}$
as an estimate of $\delta^{2}$. While this approach has the great advantage of
simplicity, one drawback is that while $\hat{\theta}_{e}-\hat{\theta}_{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 the optimal $\lambda^{*}$ on average when
$\operatorname{Cov}(\hat{\theta}_{e},\hat{\theta}_{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 theoretically
optimal. Furthermore, their method does not have a natural extension to
estimate heterogeneous causal effects.

##  3 Setup, Notation and Assumptions

###  3.1 Setup and notations

We consider two data sources, experimental data $\mathcal{D}_{e}$ of size
$n_{e}$ and observational data $\mathcal{D}_{o}$ of size $n_{o}$. Typically,
we would expect $n_{e}<n_{o}$ to represent a common scenario in practice where
researchers want to augment a small-scaled RCT trial with much larger
observational data, yet this is not required. The causal models for the two
data sources are represented by the graphs in Figure 1, where
$T\in\mathcal{R}$ represents a treatment, $\bm{W}\in\mathcal{R}^{d}$ is the
set of pre-treatment covariates, $Y\in\mathcal{R}$ is the outcome of interest
and $U$ represents unmeasured confounding. In both data sources, we assume
that the observations $X_{i}=\left(Y_{i},T_{i},\bm{W}_{i}\right)$ are
independently and identically distributed, respectively. We are interested in
the causal question “what would happen to $Y$ if we set $T=t$ by external
intervention?”. In particular, we are interested in the heterogeneous
treatment effect, that is “For an individual of certain characteristics, what
is the average effect of a treatment $T$ on their outcome $Y$?”.

There are several overlapping frameworks to represent causal relationships,
including potential outcomes (Rubin, 1974), Structural Causal Models (e.g.
Pearl, 2009) that are 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). The discussion and results in this
paper generally apply independent of the choice of frameworks. For notational
purposes, we use Pearl’s ‘$do(\cdot)$’ operator to indicate interventions. Let
the observations $X_{i}$ be distributed according to $P$ with some density
$p$. The observational and experimental data come from different distributions
with some aspects in common and we use $p_{o}$ and $p_{e}$ to denote the two.

$W$$T$$Y$$U$(a)$W$$T$$Y$$U$(b) Figure 1: Causal models for (a) the
experimental data $\mathcal{D}_{e}$ and (b) the observational data
$\mathcal{D}_{o}$. Dashed edges denote possible causal relationships. For
example, In (a), if edge $W\to T$ exists, it means that the treatment is
conditionally randomized, otherwise completely randomized.

######  Remark 3.1.

In the experimental data $\mathcal{D}_{e}$ following graph (a) in Figure 1,
the covariate set $\bm{W}$ satisfies the back-door criterion and the causal
effect of $T$ on $Y$ is identifiable by

$P_{e}(y\,|\,do(t))=\sum_{\bm{w}}P(y\,|\,t,\bm{w})P(\bm{w}).$

In contrast, in the observational data $\mathcal{D}_{o}$ following graph (b)
in Figure 1, conditional on $\bm{W}$, there is still a spurious back-door path
$T\leftarrow U\to Y$ open because $U$ remains unmeasured. As a result,
$P_{o}\left(y\,|\,do(t)\right)$ is not identifiable in $\mathcal{D}_{o}$.

Based on the setup depicted in Figure 1, we can further divide the set of
observed covariates $\bm{W}$ into two sets, $\bm{Z}$ and $\bm{C}$. The set
$\bm{C}\subseteq\bm{W}$ denotes the causal effect modifiers we are interested
in and $\bm{Z}=\bm{W}\setminus\bm{C}$ represents the remaining covariates that
we want to marginalize over. Specifically, we are interested in the
interventional distribution

$P\left(Y=y\,|\,do(T=t),\bm{C}=\bm{c}\right)$ (1)

which refers to the conditional distribution given $\bm{C}=\bm{c}$, but
marginal over $\bm{Z}$, under an experiment where $T$ is fixed by intervention
to the value $t$. The distribution (1) is sometimes also denoted as the
distribution of the potential outcome $Y^{(t)}\,|\,\bm{C}=\bm{c}$. This
$\bm{C}$-specific interventional distribution allows us to examine the effect
modification by set $\bm{C}$. This conditional distribution describes the
heterogeneous treatment effects, and in the case where treatment $T$ is
dichotomous, we can subsequently derive the CATE.

######  Example 3.2.

We use an example to further illustrate our setup with the distinction of
covariate sets $\bm{C}$ and $\bm{W}$. In a medical context, clinicians may be
interested in the efficacy of a drug on patients with certain characteristics.
Suppose we have access to an RCT and an observational dataset, which both
record whether individuals took the drug $T$, some health outcome $Y$, body
mass index (BMI) $C$ and blood glucose level $Z$. The causal relationships are
represented in Figure 2, where we assume that both $C$ and $Z$ influence the
outcome $Y$. In addition, there is a link between BMI $C$ and blood glucose
level $Z$. To give guidance on the prescription of this drug $T$, we are
interested in the causal effect moderated by BMI $C$ and want to marginalize
it over the distribution of $Z$. Specifically, we are interested in

$\displaystyle\tau({c})$
$\displaystyle=\mathbb{E}[Y\,|\,do(T=1),{c}]-\mathbb{E}[Y\,|\,do(T=0),{c}]$
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[Y\,|\,do(T=1),{c},Z]-\mathbb{E}[Y\,|\,do(T=0),{c},Z\right]\right].$

We will revisit this example in Section 5, where we conduct a simulation study
following this setup.

$T$$Y$$Z$$U$$C$(a)$T$$Y$$Z$$U$$C$(b) Figure 2: Causal models for (a) the
experimental data $\mathcal{D}_{e}$ and (b) the observational data
$\mathcal{D}_{o}$ in Example 3.2.

###  3.2 Assumptions

We then state our main assumptions below.

######  Assumption 1.

(Positivity) $P(T=t\,|\,\bm{C}=\bm{c})>0$ for all values $\bm{c}$ with
$P(\bm{C}=\bm{c})>0$ in $\mathcal{D}_{e}$.

This positivity assumption ensures that for each individual in the
experimental sample, there is a nonzero probability of being assigned to each
of the treatment levels.

######  Assumption 2.

$P_{e}\left(Y=y\,|\,do(T=t),\bm{C}=\bm{c}\right)=P_{o}\left(Y=y\,|\,do(T=t),\bm{C}=\bm{c}\right)$

Assumption 2 assumes that the conditional causal distribution with
intervention on $T$ is the same across $\mathcal{D}_{e}$ and
$\mathcal{D}_{o}$, even though it is generally not identifiable in the latter.
We want to highlight that we only require the conditional causal distribution
of interest, that is, conditional on the subset of covariates $C$ to be
invariant across the two datasets, instead of the full set including $Z$.

As illustrated in Figure 3, we use $\theta$ to denote the set of parameters of
the common conditional distribution
$P\left(Y=y\,|\,do(T=t),\bm{C}=\bm{c}\right)$, while $\mathcal{D}_{e}$ and
$\mathcal{D}_{o}$ may have different sets of parameters, say, $\gamma$ and
$\psi$, describing the rest of the joint distribution. The inference of the
common causal parameters $\theta$ is of interest and we consider $\gamma$ and
$\psi$ as nuisance parameters.

$\mathcal{D}_{e}$$\mathcal{D}_{o}$$\theta$$\gamma$$\psi$ Figure 3: A graphical
interpretation of the combination problem.

######  Remark 3.3.

As discussed in Remark 3.1, the causal distribution
$P\left(Y=y\,|\,do(T=t),\bm{C}=\bm{c}\right)$ is identifiable by
$P\left(Y=y\,|\,T=t,\bm{C}=\bm{c}\right)$ in $\mathcal{D}_{e}$ but not in
$\mathcal{D}_{o}$ due to unmeasured confounder $U$. We approach this
identifiability issue as a model misspecification problem, that is, we view
that the inference for $\theta$ is misspecified using the data model for
$P\left(Y=y\,|\,T=t,Z=z\right)$ because we fail to include the covariates $U$.
As a consequence of model misspecification, the $\theta$ estimated from the
observational data alone, $\hat{\theta}_{o}$, is prone to bias.

Through combining datasets, we want to optimally estimate $\tau(\bm{c})$ with
regard to Mean Squared Error (MSE) as the loss function. One way to look at
this problem is that we are combining $\hat{\theta}_{e}$, which has no bias
but high variance due to the smaller sample size, and $\hat{\theta}_{o}$,
which is possibly biased due to unobserved confounding but low in variance.
Essentially, we aim to optimize a bias-variance trade-off.

##  4 A Power Likelihood Approach

Traditional Bayesian inference assumes that the data model, $f(x;\theta)$ is
correct up to the known parameter value $\theta$. However, Bayesian inference
loses its predictive optimality (Zellner, 1988) when the data models are
misspecified. Bissiri et al. (2016) proposed a framework for general Bayesian
inference where parameters are connected to observations through a loss
function rather than the traditional likelihood functions, which relaxes the
requirement of a ‘true data-generating mechanism’:

$\pi(\theta\mid
x)\propto\exp\left\\{-l\left(\theta,x\right)\right\\}\pi(\theta),$

where $l\left(\theta,x\right)$ is a loss function. Under this overarching
framework, a popular solution to robustly allow for Bayesian learning under
model misspecification is to raise the likelihood to a fractional power. This
approach is discussed as a “power prior” in Ibrahim and Chen (2000), a “data-
modified prior” in Walker and Hjort (2001) and a “power likelihood” in Holmes
and Walker (2017).

$\pi_{\eta}(\theta\mid
X)\propto\pi(\theta)\prod_{i=1}^{n}f(X_{i};\theta)^{\eta}.$

Adapting to our inference problem, we take the joint likelihood but raise the
likelihood of the observational data to a power $\eta$:

$f_{\eta}(\bm{X_{e}},\bm{X_{o}};\eta,\psi,\gamma)=f_{e}(\bm{X_{e}};\theta,\gamma)\times\left(f_{o}\left(\bm{X_{o}};\theta,\psi\right)\right)^{\eta}.$

where $\bm{X_{e}}$ and $\bm{X_{o}}$ represent observations in
$\mathcal{D}_{e}$ and $\mathcal{D}_{o}$, and $p_{e}(\bm{X_{e}};\theta,\gamma)$
and $p_{o}\left(\bm{X_{o}};\theta,\psi\right)$ are the joint densities
respectively. This is equivalent to, under the general Bayesian framework,
defining a loss function:

$l_{\eta}(\theta,\psi,\gamma;\bm{X_{e}},\bm{X_{o}})=-\left\\{\log
f_{e}\left(\bm{X_{e}};\theta,\gamma\right)+\eta\times\log
f_{o}\left(\bm{X_{o}};\theta,\psi\right)\right\\}.$

Then by Bayes’ rule, the posterior distribution of $\theta$ and nuisance
parameters $\psi$ and $\gamma$ becomes

$\displaystyle\pi\left(\theta,\psi,\gamma\mid\bm{X_{e}},\bm{X_{o}}\right)$
$\displaystyle\propto\exp\left\\{-l_{\eta}\left(\theta,\psi,\gamma;\bm{X_{e}},\bm{X_{o}}\right)\right\\}\pi(\theta,\psi,\gamma)$
$\displaystyle=f_{e}(\bm{X_{e}};\theta,\gamma)\cdot
f_{o}\left(\bm{X_{o}};\theta,\psi\right)^{\eta}\cdot\pi(\theta,\psi,\gamma),$
(2)

where $\pi(\theta,\psi,\gamma)$ is the joint prior distribution of the
parameters. We can see that the parameter $\eta$ works like a dial moderating
the influence from the observational data, relative to that from the
experimental data. When $\eta=0$, the loss function is just the negative log
likelihood of the experimental data $\mathcal{D}_{e}$ and the influence of the
observational data is cut off as if $\mathcal{D}_{o}$ is excluded from the
analysis. When $\eta=1$ then it is the conventional Bayesian posterior and
means that we treat $\mathcal{D}_{e}$ and $\mathcal{D}_{o}$ equally for
inference. Intuitively, we want to choose $\eta$ to be between $0$ and $1$.

Two important components of the power likelihood approach are the likelihood
and the power $\eta$. In the remainder of this section, we discuss the
likelihood we propose to use and the robust selection of $\eta$.

###  4.1 Frugal parameterization

To specify the joint densities $p_{e}$ and $p_{o}$ in (2), we propose to adopt
the frugal parameterization (Evans and Didelez, 2021), which is designed
specially for causal inference applications. Following its framework, we can
break down a joint distribution $p(\bm{z},t,y\,|\,\bm{c})$ into three separate
pieces:

  * •

the distribution of ‘the past’:
$p_{\bm{Z}T\,|\,\bm{C}}(\bm{z},t\,|\,\bm{c}):=P(\bm{Z}=\bm{z},T=t\,|\,\bm{C}=\bm{c})$

  * •

the causal quantity of interest: $p_{Y\mid T\bm{C}}^{*}(y\mid
t,\bm{c}):=P\left(Y=y\mid do(T=t),\bm{C}=\bm{c}\right)$, and

  * •

a dependence measure $\phi^{*}_{Y\bm{Z}\,|\,T\bm{C}}$ between $Y$ and $\bm{Z}$
conditional on $T$ and $\bm{C}$.

Examples of such dependence structures include copulas and conditional odds
ratios. The advantages of using frugal parameterization are three-fold.
Firstly, it gives a likelihood with causal parameterization to the
observational data. Secondly, it isolates the causal quantity of interest,
$p_{Y\mid T}^{*}(y\mid t)$, from the rest of the joint distribution, which
means that we can directly target inference $\theta$, which parameterizes this
distribution, and treat the data-specific parameters $\gamma$ and $\psi$ as
nuisance parameters. Finally, the frugal parameterization allows the non-
causal distributions to differ between the RCT and observational data, for
example, the joint distribution of covariates and treatment assignment,
$p\left(\bm{z},t\right)$.

Under the frugal parameterization, we can factorize joint densities $p_{e}$
and $p_{o}$ as below:

$\displaystyle p_{e}\left(y,\bm{z},t\,|\,\bm{c};\theta,\gamma\right)$
$\displaystyle=p_{e,\bm{Z}T\,|\,\bm{C}}\left(\bm{z},t\,|\,\bm{c};\gamma\right)\cdot
p_{Y\mid T\bm{C}}^{*}(y\mid t,\bm{c};\theta)\cdot\phi^{*}_{e,Y\bm{Z}\mid
T\bm{C}}\left(y,\bm{z}\,|\,t,\bm{c};\gamma\right)$ (3) $\displaystyle
p_{o}\left(y,\bm{z},t\,|\,\bm{c};\theta,\psi\right)$
$\displaystyle=p_{o,\bm{Z}T\,|\,\bm{C}}\left(\bm{z},t\,|\,\bm{c};\psi\right)\cdot
p_{Y\mid T\bm{C}}^{*}(y\mid t,\bm{c};\theta)\cdot\phi^{*}_{o,Y\bm{Z}\mid
T\bm{C}}\left(y,\bm{z}\,|\,t,\bm{c};\psi\right),$ (4)

where, by Assumption 2, we have

$p_{Y\mid T\bm{C}}^{*}(y\mid t,\bm{c};\theta):=p_{e,Y\mid T\bm{C}}^{*}(y\mid
t,\bm{c};\theta)=p_{o,Y\mid T\bm{C}}^{*}(y\mid t,\bm{c};\theta).$

Ergo this parameterization explicitly isolates the causal distribution of
interest, $p_{Y\mid T\bm{C}}^{*}(y\mid t,\bm{c};\theta)$ from the remainder of
the joint distribution. Additionally, to ensure a consistent estimation of
CATE conditional on the covariate set $\bm{C}$ from the randomized data, we
further a correctly specified model for $p_{e,Y\mid T\bm{C}}^{*}(y\mid
t,\bm{c};\theta)$.

######  Assumption 3.

The causal marginal distribution $P\left(Y=y\mid do(T=t),\bm{C}=\bm{c}\right)$
in the randomized data $\mathcal{D}_{e}$ is correctly specified by a
parametric model $p_{e,Y\mid T\bm{C}}^{*}(y\mid t,\bm{c};\theta)$.

A convenient choice of such specification is a linear model, although
linearity is not required by our proposed method. Unlike requiring a correctly
specified data model for the full joint distribution, Assumption 3 is strictly
less stringent. Furthermore, Assumption 3 imposes a requirement only on the
randomized data, but makes no restriction on the observational data model.

######  Remark 4.1.

With Assumption 3, our inference takes place in the so-called
$\mathcal{M}$-closed world. However, as famously stated by Box (1976), “all
models are wrong”. The assumption of a correct specification is unlikely to
hold in practice. Instead, we conduct causal inference in the
$\mathcal{M}$-open world (Bernardo and Smith, 1994), where, asymptotically,
MLE estimation minimizes the Kullback-Leibler (KL) divergence between the true
data generative model to the assumed family of models. Additionally, it is
worth noting that, in the case of misspecification, if treatment is binary,
the ATE estimate constructed through marginalizing over all conditioning
variables $\bm{C}$ will be consistent nonetheless.

We then apply this parameterization to specify the power posterior in (2).
Substituting the joint densities $p_{e}(\bm{X_{e}};\theta,\gamma)$ and
$p_{o}\left(\bm{X_{o}};\theta,\psi\right)$ by (3) and (4) gives

$\displaystyle\pi_{\eta}\left(\theta,\psi,\gamma\mid\bm{X_{e}},\bm{X_{o}}\right)$
$\displaystyle\propto p_{e}(\bm{X_{e}};\theta,\gamma)\cdot
p_{o}\left(\bm{X_{o}};\theta,\psi\right)^{\eta}\cdot\pi(\theta,\psi,\gamma)$
$\displaystyle\propto
p_{e,\bm{Z}T\,|\,\bm{C}}\left(\bm{z},t\,|\,\bm{c};\gamma\right)\cdot p_{Y\mid
T\bm{C}}^{*}(y\mid t,\bm{c};\theta)\cdot\phi^{*}_{e,Y\bm{Z}\mid
T\bm{C}}\left(y,\bm{z}\,|\,t,\bm{c};\gamma\right)$
$\displaystyle\quad\times\left(p_{o,\bm{Z}T\,|\,\bm{C}}\left(\bm{z},t\,|\,\bm{c};\psi\right)\cdot
p_{Y\mid T\bm{C}}^{*}(y\mid t,\bm{c};\theta)\cdot\phi^{*}_{o,Y\bm{Z}\mid
T\bm{C}}\left(y,\bm{z}\,|\,t,\bm{c};\psi\right)\right)^{\eta}$
$\displaystyle\quad\times\pi(\theta,\psi,\gamma).$ (5)

###  4.2 Choosing the optimal influence factor $\eta$

We can see that $\eta$ plays a crucial role in controlling how much
information we want to borrow from the observational data. The selection of
$\eta$ needs to be performed carefully: the observational data can be useful
in reducing the variance, however, we do not want to contaminate our causal
estimate with excessive bias by setting $\eta$ too high. Ideally, we want
$\eta$ to be chosen robustly adaptive to the data, where we incorporate more
information when the intrinsic bias in the estimate from the observational
data, $\theta_{o}$, is small or its variability is low, and vice versa.
Specifically, we use MSE as the loss function, and we aim at lowering the risk
of our resultant estimator $\hat{\theta}_{\eta}$,

$R(\hat{\theta}_{\eta},\theta)=\mathbb{E}\left(||\hat{\theta}_{\eta}-\theta||^{2}\right).$

In recent literature, several methods for selecting the $\eta$ in the power
likelihood have been proposed, mainly in the context of addressing model
misspecification, for example, the SafeBayes algorithm (Grünwald and Van
Ommen, 2017), expected information matching (Holmes and Walker, 2017; Lyddon
et al., 2019), and frequentist coverage probability calibration (Syring and
Martin, 2019). A review and empirical comparison of these methods can be found
in Wu and Martin (2022).

Considering our problem setting and objectives, similar to the choice of
Carmona and Nicholls (2020), we propose to select $\eta$ by maximizing the
expected log point-wise predictive density (ELPD):

$\displaystyle\operatorname{ELPD}(\eta)$
$\displaystyle=\mathbb{E}_{\tilde{\mathfrak{X}}}\log
p_{\eta}(\tilde{X}\,|\,\bm{x}_{e},\bm{x}_{o})$
$\displaystyle=\int_{\tilde{\mathfrak{X}}}p_{t}(\tilde{x})\log
p_{\eta}(\tilde{x}\,|\,\bm{x}_{e},\bm{x}_{o})\,d\tilde{x},$ (6)

where

$p_{\eta}(\tilde{x}\,|\,\bm{x}_{e},\bm{x}_{o})=\int_{\Gamma}\int_{\Theta}p(\tilde{x}\,|\,\theta,\gamma)p(\theta,\gamma\,|\,\bm{x}_{e},\bm{x}_{o})\,d\theta\,d\gamma$

is the posterior predictive distribution indexed by $\eta$, and the
expectation is with respect to the ‘true’ data generating process $p_{t}$. As
the ‘true’ distribution is unknown, there are different ways to estimate ELPD.
A commonly used method to approximate the ‘true’ distribution is leave-one-out
cross-validation (LOO), which involves leaving one observation out at a time,
evaluating the posterior predictive density on this observation and then
averaging over all observations. This method is computationally expensive as
one can imagine. Vehtari et al. (2017) introduced an efficient computation of
LOO using Pareto-smoothed importance sampling (LOO-PSIS) which avoids repeated
partition of the dataset.

Another estimation method is to use the widely applicable information
criterion (WAIC) which is proven to be asymptotically equal to LOO (Watanabe,
2010). The WAIC method simply uses the ordinary posterior to estimate the
density of each observation, and subtracts the effective number of parameters,
defined as

$\displaystyle\widehat{\operatorname{ELPD}}(\eta)$
$\displaystyle=\frac{1}{n_{e}}\sum_{i=1}^{n_{e}}\log\hat{p}_{\eta}(x_{i}\,|\,\bm{x})-\hat{d}_{\text{WAIC}},$

where $\hat{p}_{\eta}(x_{i}\,|\,\bm{x})$ is estimated using the posterior
samples of the parameters. $\hat{d}_{\text{WAIC}}$ is the estimated effective
number of parameters. Gelman et al. (1995) provided a mean-based and a
variance-based definitions for $\hat{d}_{\text{WAIC}}$, and consistent with
Vehtari et al. (2017), we use the variance-based definition

$d_{\text{WAIC }}=\sum_{i=1}^{n_{e}}\operatorname{Var}_{\text{post
}}\left(\log p\left(x_{i}\mid\theta\right)\right),$

where $\operatorname{Var}_{\text{post }}$ is the posterior variance of the log
predictive density for each data point $x_{i}$ in the randomized data. We
propose to evaluate the ELPD on $\mathcal{D}_{e}$ as we assume the causal
relationship in the experimental data is correct. We then apply a grid search
in $[0,1]$ to find the $\eta^{*}$ such that

$\eta^{*}=\operatorname*{arg\,max}_{\eta}\widehat{\operatorname{ELPD}}(\eta).$

######  Remark 4.2.

It is worth noting the connection between ELPD and KL divergence measure. We
can rewrite (6) as

$\displaystyle\operatorname{ELPD}(\eta)$
$\displaystyle=\int_{\mathfrak{X}}p_{t}(\tilde{x})\log
p_{\eta}(\tilde{x}\,|\,\bm{x}_{e},\bm{x}_{o})\,d\tilde{x}$
$\displaystyle=\int_{\mathfrak{X}}p_{t}(\tilde{x})\left(\log\frac{p_{\eta}(\tilde{x}\,|\,\bm{x}_{e},\bm{x}_{o})}{p_{t}(\tilde{x})}+\log
p_{t}(\tilde{x})\right)\,d\tilde{x}$
$\displaystyle=\int_{\mathfrak{X}}p_{t}(\tilde{x})\log\frac{p_{\eta}(\tilde{x}\,|\,\bm{x}_{e},\bm{x}_{o})}{p_{t}(\tilde{x})}\,d\tilde{x}+\int_{\mathfrak{X}}p_{t}(\tilde{x})\log
p_{t}(\tilde{x})\,d\tilde{x}$
$\displaystyle=-d_{\operatorname{KL}}\left(p_{t}\,|\,p_{\eta}\right)+c,$ (7)

where $c=\int_{\mathfrak{X}}p_{t}(\tilde{x})\log
p_{t}(\tilde{x})\,d\tilde{x}$, which is the entropy of $p_{t}$, is a constant
independent of $\eta$. Essentially, choosing the $\eta$ via maximizing ELPD is
equivalent to selecting the posterior predictive distribution, indexed by
$\eta$, that is closest to the true data distribution in KL-divergence.

Our proposed method to select $\eta$ is summarized in Algorithm 1.

Algorithm 1 Power likelihood $\eta$ selection

Initialize $\eta^{*}\leftarrow 0$ , $\operatorname{ELPD}^{*}\leftarrow-\infty$

for $i=0,1,2,\ldots,N$ do

Let $\eta=i/N$

Sample
$(\theta^{(i)},\psi^{(i)},\gamma^{(i)})\sim\pi_{\eta}\left(\theta,\psi,\gamma\mid\bm{X}_{e},\bm{X}_{o}\right)$
using any appropriate sampler

Compute $\widehat{\operatorname{ELPD}}(\eta)$ using posterior samples
$(\theta^{(i)},\psi^{(i)},\gamma^{(i)})$

if $\widehat{\operatorname{ELPD}}(\eta)>\operatorname{ELPD}^{*}$ then

$\operatorname{ELPD}^{*}\leftarrow\widehat{\operatorname{ELPD}}(\eta)$

$\eta^{*}\leftarrow\eta$

end if

end for

return $\eta^{*}$

Once an optimal $\eta$ is selected, we can then perform likelihood-based
inference.

###  4.3 Consistency

Intuitively, the $\hat{\eta}$ that maximizes ELPD in (7) depends on the
magnitude of bias introduced by the observational dataset. The more biased
$\mathcal{D}_{o}$ is, the less we should incorporate its influence. We are
particularly interested in studying the asymptotics of $\hat{\eta}$ and the
resulting estimator $\hat{\theta}_{\hat{\eta}}$. Similar to the discussion in
Yang and Ding (2020) and Dang et al. (2022), we assume that the bias in the
observational data $\mathcal{D}_{o}$ due to hidden confounding is not fixed
but depends on its sample size, even though this usually is not practically
true. Let $\delta=\delta(n_{o})=\delta^{*}/n_{o}^{k}$ denote the bias in
$\mathcal{D}_{o}$. Specifically, we consider the following two scenarios:

Scenario 1

$k\geq 1/2$ or $\delta^{*}=0$. This is where the observational data is
unbiased or the bias reduces at a rate faster than $\sqrt{n}$, which becomes
negligible as we apply the Central Limit Theorem (CLT).

Scenario 2

$0\leq k<1/2$ and $\delta^{*}\neq 0$. This is where the non-zero bias is a
constant or the bias reduces at a rate slower than $\sqrt{n}$, which does not
vanish as we apply the CLT.

######  Theorem 4.3 (Consistency).

Consider experimental data $X=(X_{1},\dots,X_{n_{e}})$ where $X_{i}$’s are
i.i.d. following a normal distribution
$X_{i}\sim\mathcal{N}(\theta^{*},\sigma^{2})$, and observational data
$Y=(Y_{1},\dots,Y_{n_{o}})$ i.i.d. following
$Y_{i}\sim\mathcal{N}(\theta^{*}+\delta,\sigma^{2})$. Here $\delta$ denotes
the unknown bias and is a function of $n_{o}$. Suppose the variance
$\sigma^{2}$ is known. The estimator following the proposed power likelihood
approach is $\hat{\theta}_{\hat{\eta}}$. Let $n=n_{e}+n_{o}$ and
$n_{e}/n_{o}=c>0$.

  * •

If $0\leq k<1/2$ and $\delta^{*}\neq 0$, the $\hat{\eta}=O(n^{-1/2})$,
$\hat{\theta}_{\hat{\eta}}=\mu^{*}+O(n^{-1/2})$ and
$\mathbb{E}[(\hat{\theta}_{\eta}-\theta^{*})^{2}]=\sigma^{2}/n_{e}+O(n^{-3/2})$.

  * •

If $k\geq 1/2$ or $\delta^{*}=0$ then $\hat{\eta}$ does not converge.

Theorem 4.3 suggests that when bias is non-zero, the selected influence factor
$\hat{\eta}$ converges to 0 and the estimator $\hat{\mu}_{\eta}$ is
$\sqrt{n}$-consistent. As a result, the MSE of the estimator
$\hat{\mu}_{\eta}$ converges to that from using the experimental data only. To
provide intuition for this theorem, as $n_{o}$ grows, the fraction of its
influence we want to incorporate reduces and as $n_{e}$ grows, the need to
augment it with external data decreases at the same time. However, when the
observational data is inherently unbiased, the selected $\hat{\eta}$ will not
converge to zero. On the contrary, as bias diminishes close to zero, we will
very likely infer that $\hat{\eta}=1$ is the best value as we restrict $\eta$
to be between $0$ and $1$. The proof is presented in Appendix A.

###  4.4 Normal approximations to the posterior distribution

Our proposed method, as described in Algorithm 1, involves sampling from
posterior distributions multiple times in the search for the optimal $\eta$.
While a standard approach is to obtain the posterior samples by Markov-Chain
Monte Carlo (MCMC) with a prior, it can be computationally challenging when
data are massive, $\mathcal{D}_{o}$ in particular. The computational challenge
is exacerbated when the parameter vector $(\theta,\psi,\gamma)$ has high
dimensionality which may affect the efficiency of the sampling procedure.

In such scenarios, we propose to apply a normal approximation of the posterior
and then sample from this distribution (Gelman et al., 1995), derived as

$\pi_{\eta}\left(\theta,\psi,\gamma\mid
X_{e},X_{o}\right)\approx\mathcal{N}\left(\hat{\theta}_{\eta},[\mathcal{I}(\hat{\theta}_{\eta})]^{-1}\right),$

where $\hat{\theta}_{\eta}$ is the MLE estimate with regard to the
$\eta$-powered likelihood function and $\mathcal{I}(\hat{\theta}_{\eta})$ is
the corresponding Fisher information matrix. A detailed derivation of these
quantities can be found in Appendix B. A nice property of this approximation
is that the bigger the data sets, the better the approximation, which will
clear the computational hurdle and make our method scalable for the
application to big data.

##  5 Simulation

We now run a simulation that illustrates how our method robustly moderates the
influence of the observational data through the $\eta$-powered likelihood and
compare the performance of our method with several existing approaches.

###  5.1 Simulation setup

Following the same setup as described in Example 3.2, and represented in
Figure 2. Suppose we observe a continuous covariate $C$, for example, body
mass index (BMI), a continuous covariate $Z$, for instance, blood glucose
level, which is associated with BMI, a binary treatment variable $T$
indicating whether a drug is taken, and a continuous health outcome $Y$ such
as a change in average blood glucouse level, for example, measured by HbA1c.
There also exists an unobserved binary confounder $U$ that influences both the
treatment and the outcome. The influence of $U$ means that the causal estimate
from the observational data is biased. Specifically, we model the
distributions as follows:

$\displaystyle U$ $\displaystyle\sim\operatorname{Bernoulli}(0.5)$
$\displaystyle C$ $\displaystyle\sim N(0,1)$ $\displaystyle Z\mid C$
$\displaystyle\sim N(\mu_{z},\,1)$ $\displaystyle T\mid C,Z,U$
$\displaystyle\sim\operatorname{Bernoulli}(\mu_{t})$ $\displaystyle Y(t)\mid
C$ $\displaystyle\sim N(\mu_{y},\,1).$

where,

$\displaystyle\mu_{z}$ $\displaystyle=0.2+0.6\,C$
$\displaystyle\operatorname{logit}\mu_{t}$
$\displaystyle=0.5+0.1\,C+0.6\,Z+0.4\,C\,Z+\psi U$ $\displaystyle\mu_{y}$
$\displaystyle=0.6+0.2\,C\phantom{.+0.0T}+1.1\,C\,T+\psi U$

Additionally, let the dependence structure between $Y$ and $Z$ given $T=t$ and
$C=c$, $\phi^{*}_{YZ\mid TC}$, be a conditionally bivariate Gaussian copula,
with correlation parameter $\rho_{t}=2\operatorname{expit}(1+2.5t)-1$.

We use the same parameterization for the experimental data $\mathcal{D}_{e}$
except that we assume that treatment is randomly assigned and hence replace
the mean of $T$ with $\mu_{t}=0.5$. We set the sample sizes to 2,500 and 250
for the observational and experimental data, respectively. The results we show
in this section are averaged across 500 sets of synthetic datasets.

######  Remark 5.1.

In this as well as other simulations that we have tested, we compared the ELPD
estimated by LOO-PSIS and WAIC method, implemented through the loo package in
R (Vehtari et al., 2017), and they almost always give very close results. We
use the WAIC method in this simulation. We have also compared obtaining
posterior parameter samples using the Metropolis-Hasting within Gibbs MCMC
sampler, with sampling directly from the approximated normal distribution as
described in Section 4.4. With sample sizes of 250 and 2,500, in a
12-dimensional parameter space, the two approaches lead to very similar
$\eta$’s being selected. For run time reasons, we sample directly from the
approximated normal distribution in this simulation study.

Figure 4: Simulation results for scenarios where $U$ has no influence
($\psi=0$), small influence ($\psi=0.75$) and large influence ($\psi=1$) on
$T$ and $Y$. Left column: average negative ELPD with $\eta\in[0,1]$. The
yellow dotted lines mark the average $\eta$ selected. Middle column: MSE of
the ATE estimate. Right Column: MSE of the CATE estimate in the subgroup with
$C>0$. The MSE is decomposed into bias (red-shaded area) and variance (blue-
shaded area).

###  5.2 Moderating the bias-variance trade-off

Figure 4 presents three scenarios where unmeasured confounding has an
influence of different magnitude: no influence ($\psi=0$), small influence
($\psi=0.75$) and large influence ($\psi=1)$. The first column shows the
average negative ELPD at each value of $\eta$. Our proposed method searches
for the $\eta$ that maximizes ELPD, which corresponds to the lowest point on
the curve. The second and third columns show the MSE of the ATE and
CATE444While we only show the CATE in the subgroup where $C>0$, we have
analyzed the MSE for CATE estimates in the subgroup $C\leq 0$ and it exhibited
a very similar pattern of MSE loss as $\eta$ varied. estimates evaluated on
$\mathcal{D}_{e}$ at different values of $\eta$. We decompose MSE into
variance (in blue) and squared bias (in red) to give an intuitive illustration
of the moderating effect of $\eta$ on the bias-variance trade-off: as we
increase $\eta$, we include more information from $\mathcal{D}_{o}$ and hence
the variance reduces, yet at the expense of an increase in bias.

When $U$ has no influence on $T$ and $Y$, that is, $\mathcal{D}_{o}$ is not
biased and is just ‘as good as’ $\mathcal{D}_{e}$, the average $\eta$ selected
is around 0.8 which is reasonable as we would want to include as much
information in $\mathcal{D}_{o}$ as possible. Unsurprisingly, the
corresponding MSEs of ATE and CATE estimates are both lower than if we do not
include $\mathcal{D}_{o}$ at all.

When $U$ has a small influence on $T$ and $Y$, our method chooses a smaller
$\eta$, averaging at around 0.7. Again, the corresponding MSEs of ATE and CATE
estimates are both lower than at $\eta=0$. When $U$ has a larger influence,
that is, $\mathcal{D}_{o}$ is more biased, our method responds by selecting an
even smaller $\eta$. Although the corresponding MSE of the ATE estimate has
increased, we achieve a slightly lower MSE of CATE estimates than $\eta=0$.

###  5.3 Comparison with existing methods

Within the same simulation, we have also compared the performance of our
method to four existing approaches namely, the MSE-minimizing estimator by
Oberst et al. (2022), shrinkage estimators by Green et al. (2005) and by
Rosenman et al. (2020), and experimental grounding by Kallus et al. (2018).

Figure 5: A comparison of relative overall MSE of CATE estimates of the 10
strata. The reference line at 1 represents the overall MSE of using IPW
estimates in $\mathcal{D}_{e}$.

We first compare the methods’ performance in estimating heterogeneous
treatment effects. We first stratify $\mathcal{D}_{e}$ and $\mathcal{D}_{o}$
by values of $C$ and normal quantiles of $Z$ into 10 strata. Then we measure
the overall MSE of the CATE estimates of the 10 strata. Figure 5 plots the
overall MSE relative to using inverse propensity weighted (IPW) estimators
from $\mathcal{D}_{e}$ only. As the MSE-minimizing estimator of Oberst et al.
(2022) does not have a natural extension to systematically estimating
heterogeneous treatment effects, it is not included in this comparison.
Additionally, the purple dashed line shows the performance of using the
estimates from $\mathcal{D}_{e}$ only but with a correctly-specified
parametric model. The difference between this line and the reference line at 1
is whether a parametric model is used or not.

We can see that all methods lead to lower overall MSE than using the strata
mean differences in $\mathcal{D}_{e}$. Parametric methods that are correctly
specified, that is, ours and experimental grounding, outperform the shrinkage
estimators. Our method leads to the lowest overall MSE and as bias increases,
our method converges back to using the experimental data only.

Figure 6 compares the relative MSE of ATE estimates using different methods.
Generally, there is a trade-off between the base-case MSE reduction and the
worst-case increase in MSE. Our method dominates both shrinkage estimators
with a lower relative MSE regardless of the bias. The MSE of the experimental
grounding estimator of Kallus et al. (2018) is on par with using the
experimental data only.

The MSE-minimizing estimator in Oberst et al. (2022) has the most comparable
performance to ours. Oberst et al. (2022) observed that, in their simulations,
the worst-case relative MSE of their estimator is bounded within a $27\%$
increase, and this claim seems to hold in our results too. Compared to their
estimator, ours achieves a larger reduction in MSE when the bias is small.
This is possibly due to the weight, $\hat{\lambda}_{ober}$, underestimating
the theoretically optimal weight $\lambda$ on average, which means that their
pooled estimator fails to include as much influence from a ‘reasonably good’
observational data as theoretically optimal. In addition, our estimator
exhibits a faster reversion to using the experimental data only when the bias
is large. This means that our method is more conservative against drawing
awfully ‘wrong’ causal conclusions by allowing excessive influence from a very
biased observational dataset contaminating the experimental data. This is an
attractive feature if a ‘safe policy’ is a priority. Policy-makers are usually
risk-averse when making crucial decisions impacting people’s lives, and prefer
a safe policy that lowers the probability of a worse outcome than the status
quo (Ben-Michael et al., 2021).

Figure 6: A comparison of relative overall MSE of ATE estimates. The reference
line at 1 represents the overall MSE of using a correctly specified parametric
model on $\mathcal{D}_{e}$ only.

###  5.4 Incorporating more variables

To examine the robustness of our method with higher dimensional covariates, we
introduce more variables into the simulation. This reflects that, in practice,
clinicians may be interested in more fine-grained CATE that is conditional on
several patient characteristics allowing for more personalized treatment
decisions. Examples of such variables include age, gender, and medical
history.

Specifically, we increased the dimension of $\bm{C}$, the set of variables
that CATE is conditioned on, to five, consisting of three continuous variables
and two binary variables. We have also introduced one more continuous variable
into $\bm{Z}$ reflecting that we have more information about a patient we want
to marginalize. A detailed model specification can be found in Appendix C.2.

Similar to Figure 5, Figure 7 compares the performance of different methods
regarding CATE estimates. We can see that our method’s performance is
comparable to that when the dimensions of the covariate sets $\bm{Z}$ and
$\bm{C}$ are low. Our method performs better than the non-parametric methods
when the model is correctly specified, and compared to only using the
randomized data $\mathcal{D}_{e}$, our method leads to a reduction in MSE when
the confounding bias is small. There is a range of elevated MSE when the bias
is moderate to high, up to an approximately 60% increase, before reverting to
using $\mathcal{D}_{e}$ only. In line with the findings in Oberst et al.
(2022), some trade-off between reduction and inflation in MSE by the magnitude
of bias seems inevitable, unless we have additional external information. For
example, Dang et al. (2022) assume there is a negative control outcome that
can be used to determine whether or not observational data should be
integrated.

Figure 7: A comparison of relative overall MSE of CATE estimates in 10 strata
when we introduce more covariates into the parametric model. The reference
line at 1 represents the overall MSE of using IPW estimates in
$\mathcal{D}_{e}$.

##  6 Analysis of STAR data

To illustrate the validity and efficacy of our proposed approach, we apply it
to a real-world dataset from Tennessee’s Student Teacher Achievement Ratio
Study (STAR). The STAR study is a longitudinal experiment where students were
randomly assigned into classes of different sizes to study the impact of class
size on students’ educational outcomes (Word et al., 1990). The STAR data
itself is from an RCT, and similar to Kallus et al. (2018), we partition the
STAR data into two and construct an observational data by artificially
introducing confounding through selectively removing a biased subset of
samples.

The STAR data covers over 7,000 students in 79 schools who were randomly
assigned into one of three class types: small class (13 to 17 students per
teacher), regular class (22 to 25 students per teacher) and regular-with-aide
class (22 to 25 students with a full-time aide) and we focus on the first two
types of interventions. Weuse standardized test scores at the end of first
grade as the outcome and use the following covariates: gender, race, date of
birth, and an indicator of whether or not a free lunch was provided. We
applied the same data processing as in Kallus et al. (2018), such as filtering
and removing records with missing values, and arrived at the same final sample
of 4,218 students, consisting of 2,811 students enrolled in rural or inner-
city schools, and 1,407 students in urban or suburban schools.

To construct the experimental datasets, we randomly sample $10\%$ of the 2,811
students enrolled in rural or inner-city schools. To construct the confounded
observational dataset, we take all the controls $(T=0)$ that were not sampled
into the experimental data, and for those receiving treatment $(T=1)$ we take
a sample of 1,000 students where we down-weight those with an outcome below
the 30th percentile. By constructing the observational data this way, we
successfully introduce confounding from two sources: 1) no students from rural
or inner-city schools enter the experimental data, this is as if there is an
unobserved exclusion criterion 2) the treatment effect is artificially biased
upwards by down-weighting units in treatment with poor outcomes. The
confounding translates into a naïve ATE estimate of 57.0 rather than the
original dataset’s 38.4.

We then apply our method to combine the constructed experimental and
observational data to estimate the ATE in the experimental data and repeat
this process for 500 different sets of datasets. Figure 8 shows an assessment
of the performance of our method. We can see a similar pattern as in the
simulation study: the shape of the negative ELPD by $\eta$ mirrors the bias-
variance trade-off. The average $\eta$ selected is 0.17 and the average MSE is
175.23, whereas if we only use the experimental data, the average MSE is
193.84 which means approximately $10\%$ reduction.

Figure 8: Results for application on the re-constructed STAR datasets, where
for the observational dataset, we down-weight those with outcomes in the
bottom 30%. Top graph: average -elpd with $\eta\in[0,1]$. Bottom graph: MSE of
the ATE estimate on the randomized data, which is decomposed into bias (red
shaded area) and variance (blue shaded area).

Figure 9 gives some intuition of our method. The figure compares the
distribution of ATE estimates using different methods. We can see that the
distribution of ATE estimates from only using the small unconfounded data,
represented by the red area, is dispersed and 0 is included in the 99%
confidence interval; this is unfortunate if high statistical power is
required. In contrast, the distribution of ATE estimates from using the large
confounded data, represented by the green area, is very concentrated around
its mean, which is fairly far away from that of the unconfounded distribution.
The blue distribution represents the distribution of the ATE estimate from our
proposed method. We can see that it sits in between the other two
distributions and is moderately spread out. In this case, 0 is no longer with
the 99% confidence interval. Recalling Remark 4.2, the learning factor $\eta$
is chosen so that the posterior predictive distribution is closest, in KL
divergence, to the “true” generative distribution of the randomized data. This
gives us a new perspective: the blue density presents the distribution of ATE
estimates that are inferred from combining information from the two datasets
so that the posterior predictive distribution is closest to the “true”
randomized data.

Figure 9: Distribution of ATE estimates from using the unconfounded data only
(red shaded area), confounded data only (green shaded area) and the combined
data using our proposed methodology (blue shaded area).The distribution is
obtained from 500 different sets of confounded and unconfounded datasets.

##  7 Discussion

In this paper, we present a novel power likelihood approach to combining
experimental and observational data that aims at optimizing the bias-variance
tradeoff for the estimation of heterogeneous treatment effects. The power
likelihood approach fits under the general Bayesian framework and introduces
an influence factor $\eta$ that moderates the impact of observational data. We
propose a principled, data-adaptive procedure to select $\eta$ by maximizing
predictive accuracy, measured by ELPD, on the randomized data. This procedure
robustly selects $\eta$ depending on how ‘contaminated’ and ‘useful’ the
observational data is for our inferential goals and has the practical
advantage that it does not require the specification of hyper-parameters.

Through empirical comparisons in Section 5, we show that our method
outperforms existing methods (Green and Strawderman, 1991; Kallus et al.,
2018; Rosenman et al., 2020; Oberst et al., 2022) for substantially lower
overall MSE for CATE estimates, a larger reduction in MSE for the ATE estimate
when bias is modest, and faster reversion to using the experimental data only
when the observational data is drastically biased. Additionally, we
demonstrate the feasibility and robustness of our proposed method as the
dimension of covariates grows. In Section 6, we apply our method to the STAR
data to illustrate its efficacy in a semi-synthetic real-world data analysis.
In the future, we intend to apply this method to actual clinical data to test
its practicality. Furthermore, although for simplicity, we assume that the
treatment is binary throughout the discussion, our proposed method, in fact,
allows for cases where treatment is multivariate or continuous. This is an
advantage over most of the existing methods.

A common limitation of Bayesian inference methods is that it is difficult to
focus on optimizing a particular parameter of interest, for example, Bayesian
dynamic borrowing methods in Ibrahim and Chen (2000) and Hobbs et al. (2012).
In contrast, in a frequentist setting, methods such as targeted learning using
TMLE (van der Laan and Rubin, 2006) and model selection using Focused
Information Criteria (FIC) (Claeskens and Hjort, 2003; Zhang and Liang, 2011),
are exactly for this purpose. In our proposal, we take a typical Bayesian
route and select $\eta$ based on expected prediction quality instead of
directly targeting the CATE estimates. As the result, our method runs the risk
of inefficiency as we include more covariates, which is manifested in the
inflation of MSE described in Section 5.4. The loss of efficiency is more
pronounced as covariates that are weakly correlated, or not correlated with
the outcome get included in the model. Again, this issue is not unique to our
method. As mitigation, we advise selecting only covariates that are predictive
of the outcome in the model and applying dimension reduction techniques such
as LASSO regression or principal component analysis (PCA) where appropriate.

##  Appendix A Proof of Theorem 4.3

######  Theorem 4.2 (Consistency).

Consider experimental data $X=(X_{1},\dots,X_{n_{e}})$ where $X_{i}$’s are
i.i.d. following a normal distribution
$X_{i}\sim\mathcal{N}(\theta^{*},\sigma^{2})$, and observational data
$Y=(Y_{1},\dots,Y_{n_{o}})$ i.i.d. following
$Y_{i}\sim\mathcal{N}(\theta^{*}+\delta,\sigma^{2})$. Here $\delta$ denotes
the unknown bias and is a function of $n_{o}$. Suppose the variance
$\sigma^{2}$ is known. The estimator following the proposed power likelihood
approach is $\hat{\theta}_{\hat{\eta}}$. Let $n=n_{e}+n_{o}$ and
$n_{e}/n_{o}=c>0$.

  * •

If $0\leq k<1/2$ and $\delta^{*}\neq 0$, the $\hat{\eta}=O(n^{-1/2})$,
$\hat{\theta}_{\hat{\eta}}=\theta^{*}+O(n^{-1/2})$ and
$\mathbb{E}[(\hat{\theta}_{\eta}-\theta^{*})^{2}]=\sigma^{2}/n_{e}+O(n^{-3/2})$.

  * •

If $k\geq 1/2$ or $\delta^{*}=0$ then $\hat{\eta}$ does not converge.

###### Proof.

A roadmap for the proof of Theorem 4.3:

  1. 1.

We start by deriving the posterior predictive distribution and ELPD

  2. 2.

We set the derivative of the ELPD to zero and find the roots

  3. 3.

For the scenario where $k\geq 1/2$ or $\delta^{*}=0$, we study the asymptotic
behaviour of $\hat{\eta}$

  4. 4.

For the scenario where $0<k<1/2$ and $\delta^{*}\neq 0$, we study the
asymptotic behaviour of $\hat{\eta}$ and the resulting estimator
$\hat{\theta}_{\hat{\eta}}$.

###  A.1 Posterior predictive distribution and the ELPD

We are interested in the parameter $\theta\in\mathcal{R}$, for which we have a
prior distribution $\pi(\theta)\sim\mathcal{N}(\theta_{0},\sigma_{0}^{2})$.
The ‘presumed’ data models of $x_{i}$ and $y_{i}$ condition on $\theta$ are:

$\displaystyle p_{e}(x_{i}\,|\,\theta)$
$\displaystyle=\frac{1}{\sqrt{(2\pi\sigma^{2})}}\exp{\left\\{-\frac{(x_{i}-\theta)^{2}}{2\sigma^{2}}\right\\}}$
(8) $\displaystyle p_{o}(y_{i}\,|\,\theta)$
$\displaystyle=\frac{1}{\sqrt{(2\pi\sigma^{2})}}\exp{\left\\{-\frac{(y_{i}-\theta)^{2}}{2\sigma^{2}}\right\\}},$
(9)

and the power posterior is

$\displaystyle\pi_{\eta}(\theta\,|\,x,y)$ $\displaystyle\propto
p_{\eta}(x,y\,|\,\theta)\times\pi(\theta)$ (10)
$\displaystyle\propto\prod_{i=1}^{n_{e}}p_{e}(x_{i}\,|\,\theta)\times\left(\prod_{i=1}^{n_{o}}p_{o}(y_{i}\,|\,\theta)\right)^{\eta}\times\pi(\theta)$
(11)
$\displaystyle\propto\exp({-\frac{n_{e}(\bar{x}-\theta)^{2}}{2\sigma^{2}}})\times\exp({-\frac{\eta
n_{o}(\bar{y}-\theta)^{2}}{2\sigma^{2}}})\times\exp({-\frac{(\theta-\theta_{0})^{2}}{2\sigma_{0}^{2}}})$
(12) $\displaystyle\propto\exp\frac{(\sigma^{2}+n_{e}\sigma_{0}^{2}+\eta
n_{o}\sigma_{0}^{2})\theta^{2}-2(\sigma^{2}\theta_{0}+n_{e}\sigma_{0}^{2}\bar{x}+\eta
n_{o}\sigma_{0}^{2}\bar{y})\theta}{2\sigma^{2}\sigma_{0}^{2}},$ (13)

where $\bar{x}=\frac{1}{n_{e}}\sum_{i=1}^{n_{e}}x_{i}$ and
$\bar{y}=\frac{1}{n_{o}}\sum_{i=1}^{n_{o}}y_{i}$ are the sample means of the
two datasets. It is easy to see that
$\pi_{\eta}(\theta;x,y)\sim\mathcal{N}(\hat{\theta},\hat{\sigma}^{2})$ where,

$\displaystyle\hat{\theta}$
$\displaystyle=\frac{\sigma^{2}\theta_{0}+n_{e}\sigma_{0}^{2}\bar{x}+\eta
n_{o}\sigma_{0}^{2}\bar{y}}{\sigma^{2}+n_{e}\sigma_{0}^{2}+\eta
n_{o}\sigma_{0}^{2}}$ (14) $\displaystyle\hat{\sigma}^{2}$
$\displaystyle=\frac{\sigma^{2}\sigma_{0}^{2}}{\sigma^{2}+n_{e}\sigma_{0}^{2}+\eta
n_{o}\sigma_{0}^{2}}.$ (15)

Asymptotically the influence from the prior distribution assumptions vanishes
and the posterior mean estimator $\hat{\theta}$ becomes

$\hat{\theta}=\hat{\theta}_{\eta}=\frac{n_{e}\bar{x}+\eta
n_{o}\bar{y}}{n_{e}+\eta n_{o}}$ (16)

which is our proposed estimator. We can see that $\hat{\theta}$ depends on the
value of $\eta$ and we write it as $\hat{\theta}_{\eta}$ to emphasise the
dependence. Similarly, the asymptotic posterior variance is

$\hat{\sigma}=\hat{\sigma}_{\eta}=\frac{\sigma^{2}}{n_{e}+\eta n_{o}}.$ (17)

For (17) to make sense, we require $n_{e}+\eta n_{o}>0$. Following the
definition in Vehtari et al. (2017), we can write $\operatorname{ELPD}$ for a
new datum as a function of $\eta$ as

$\displaystyle\operatorname{ELPD}(\eta)=\int_{\mathcal{Z}}p_{e}^{*}(z_{i})\log
p_{\eta}(z_{i}|x,y)dz_{i},$

where $p_{e}^{*}(z_{i})$ is the distribution representing the true data-
generating process of the experimental data, $z_{i}$ is a new data point from
this true distribution, and

$p_{\eta}(z_{i}|x,y)=\int
p_{e}(z_{i}|\theta,\sigma^{2})\pi_{\eta}(\theta|x,y)d\theta$

is a candidate posterior predictive distribution indexed by $\eta$.
Substituting the expressions for $p_{e}(z_{i}|\theta,\sigma^{2})$ and
$\pi_{\eta}(\theta|x,y)$ gives

$\begin{split}p_{\eta}(z_{i}|x,y)&=\int\frac{1}{\sqrt{(2\pi\sigma^{2})}}\exp{\left(-\frac{(z_{i}-\theta)^{2}}{2\sigma^{2}}\right)}\times\frac{1}{\sqrt{(2\pi\hat{\sigma}^{2})}}\exp{\left(-\frac{(\theta-\hat{\theta})^{2}}{2\hat{\sigma}^{2}}\right)}\\\
&=\frac{1}{\sqrt{2\pi(\sigma^{2}+\hat{\sigma}^{2})}}\exp{\left(-\frac{(z_{i}-\hat{\theta})^{2}}{2(\sigma^{2}+\hat{\sigma}^{2})}\right)},\end{split}$

where $\hat{\theta}$ and $\hat{\sigma}$ are the asymptotic posterior mean and
variance as in (16) and (17). We can then write out the log prediction density
and $\operatorname{ELPD}(\eta)$, given $x$ and $y$, as

$\displaystyle\log p_{\eta}(z_{i}|x,y)$
$\displaystyle=-\frac{1}{2}\log(2\pi)-\frac{1}{2}\log(\sigma^{2}+\hat{\sigma}^{2})-\frac{(z_{i}-\hat{\theta})^{2}}{2(\sigma^{2}+\hat{\sigma}^{2})}$
$\displaystyle\operatorname{ELPD}(\eta)$
$\displaystyle=\mathbb{E}_{P_{z}^{*}}\\{\log p_{\eta}(Z_{i}|x,y)\\}$
$\displaystyle=-\frac{1}{2}\log(2\pi)-\frac{1}{2}\log(\sigma^{2}+\hat{\sigma}^{2})-\mathbb{E}_{P_{z}^{*}}\left\\{\frac{(Z-\hat{\theta})^{2}}{2(\sigma^{2}+\hat{\sigma}^{2})}\right\\}.$

###  A.2 Derivative of ELPD

We want to find $\eta=\eta^{*}$ such that
$\eta^{*}=\operatorname*{arg\,max}_{\eta}\operatorname{ELPD}(\eta)$. However,
the true data generative process represented by $p_{e}^{*}$ remains unknown,
we have to empirically estimate ELPD by evaluating it over a calibration
dataset of size $n_{c}$ with i.i.d. samples from the true distribution,

$\displaystyle\widehat{\operatorname{ELPD}(\eta)}$
$\displaystyle=\frac{1}{n_{c}}\sum_{i=1}^{n_{c}}\left\\{-\frac{1}{2}\log(2\pi)-\frac{1}{2}\log(\sigma^{2}+\hat{\sigma}^{2})-\frac{(z_{i}-\hat{\theta})^{2}}{2(\sigma^{2}+\hat{\sigma}^{2})}\right\\}.$

It is important to recognize that $\eta$ is now a parameter of the likelihood
model for the calibrated data. Essentially we want to find the MLE estimator
of $\eta$, $\hat{\eta}$, such that

$\displaystyle\hat{\eta}=\operatorname*{arg\,max}_{\eta}\widehat{\operatorname{ELPD}(\eta)}.$

We consider the score function,

$\displaystyle S_{n_{c}}(\eta)$
$\displaystyle=\frac{1}{n_{c}}\sum_{i=1}^{n_{c}}S(z_{i},\eta)$
$\displaystyle=\frac{1}{n_{c}}\sum_{i=1}^{n_{c}}\frac{\partial}{\partial\eta}\left\\{-\frac{1}{2}\log(2\pi)-\frac{1}{2}\log(\sigma^{2}+\hat{\sigma}^{2})-\frac{(z_{i}-\hat{\theta})^{2}}{2(\sigma^{2}+\hat{\sigma}^{2})}\right\\}$

for which $\mathbb{E}[S_{n_{c}}(\eta_{o})]=0$. Suppose that $\hat{\eta}$ is a
solution to the estimation equation $S_{n_{c}}(\eta)=0$, i.e.
$S_{n_{c}}(\hat{\eta})=0$. By semi-parametric inference theory,

$\sqrt{n_{c}}\left(\hat{\eta}-\eta_{0}\right)\xrightarrow{D}\mathcal{N}\left(0,\left[\mathbb{E}\left\\{\frac{\partial
S\left(Z,\eta_{0}\right)}{\partial\eta^{T}}\right\\}\right]^{-1}\operatorname{Var}\left\\{S\left(Z,\eta_{0}\right)\right\\}\left[\mathbb{E}\left\\{\frac{\partial
S\left(Z,\eta_{0}\right)}{\partial\eta^{T}}\right\\}\right]^{-1^{T}}\right).$

Under suitable regularity conditions, we have the relationship between the
Fisher Information $\mathcal{I}\left(\eta_{0}\right)$ and the terms in the
variance as below:

$\mathcal{I}^{-1}\left(\eta_{0}\right)=\mathbb{E}\left\\{-\frac{\partial
S\left(Z,\eta_{0}\right)}{\partial\eta^{T}}\right\\}^{-1}=\mathbb{E}\left\\{S_{\eta}\left(Z,\eta_{0}\right)S_{\eta}^{T}\left(Z,\eta_{0}\right)\right\\}=\operatorname{Var}\left\\{S\left(Z,\eta_{0}\right)\right\\}.$

Therefore,

$\sqrt{n_{c}}\left(\widehat{\eta}_{n_{c}}-\eta_{0}\right)\xrightarrow{d}\mathcal{N}\left(0,\mathcal{I}\left(\eta_{0}\right)^{-1}\right).$

To find $\eta_{0}$, we want to solve the equation
$\mathbb{E}[S_{n_{c}}(\eta_{o})]=0$. Expanding it, we have

$\displaystyle\mathbb{E}[S_{n_{c}}(\eta_{o})]$
$\displaystyle=\mathbb{E}_{Z}\left\\{\frac{1}{n_{c}}\sum_{i=1}^{n_{c}}\frac{\partial}{\partial\eta}\left[-\frac{1}{2}\log(2\pi)-\frac{1}{2}\log(\sigma^{2}+\hat{\sigma}^{2})-\frac{(Z_{i}-\hat{\theta})^{2}}{2(\sigma^{2}+\hat{\sigma}^{2})}\right]\right\\}$
$\displaystyle=-\frac{1}{2}\frac{\partial}{\partial\eta}\log(\sigma^{2}+\hat{\sigma}^{2})-\frac{1}{2}\mathbb{E}_{Z}\left[\frac{\partial}{\partial\eta}\frac{(Z-\hat{\theta})^{2}}{\sigma^{2}+\hat{\sigma}^{2}}\right]$
$\displaystyle=-\frac{1}{2}\left\\{\frac{\partial}{\partial\eta}\left[\log(\sigma^{2}+\hat{\sigma}^{2})\right]+\frac{\partial}{\partial\eta}\mathbb{E}_{Z}\left[\frac{(Z-\hat{\theta})^{2}}{\sigma^{2}+\hat{\sigma}^{2}}\right]\right\\}$
$\displaystyle=-\frac{1}{2}\left\\{\frac{\partial}{\partial\eta}\log(\sigma^{2}+\hat{\sigma}^{2})+\frac{\partial}{\partial\eta}\frac{(\theta^{*}-\hat{\theta})^{2}+\sigma^{2}}{\sigma^{2}+\hat{\sigma}^{2}}\right\\}$
(18)

The third equality is by Leibniz integral rule and dominated convergence
theorem. We also have the derivatives of $\hat{\theta}$ and $\hat{\sigma}^{2}$
with respect to $\eta$ as follows:

$\displaystyle\frac{\partial\hat{\theta}}{\partial\eta}$
$\displaystyle=\frac{n_{e}n_{o}\left(\bar{Y}-\bar{X}\right)}{\left(n_{e}+\eta
n_{o}\right)^{2}}$ (19)
$\displaystyle\frac{\partial\hat{\sigma}^{2}}{\partial\eta}$
$\displaystyle=-\frac{n_{o}\sigma^{2}}{\left(n_{e}+\eta n_{o}\right)^{2}}.$
(20)

Substituting Equation (19) - (20) into (18), gives

$\begin{split}\mathbb{E}[S_{n_{c}}(\eta_{0})]=-&\frac{n_{o}\sigma^{2}}{2(\sigma^{2}+\hat{\sigma}^{2})^{2}\left(n_{e}+\eta
n_{o}\right)^{4}}\left\\{2\left(n_{e}\left(\bar{X}-\theta^{*}\right)+\eta
n_{o}\left(\bar{Y}-\theta^{*}\right)\right)n_{e}\left(\bar{Y}-\bar{X}\right)\left(n_{e}+\eta
n_{o}+1\right)\right\\}\\\ &+\left(n_{e}\left(\bar{X}-\theta^{*}\right)+\eta
n_{o}\left(\bar{Y}-\theta^{*}\right)\right)^{2}-\left(n_{e}+\eta
n_{o}\right)\sigma^{2}\\}\end{split}$ (21)

Setting $\mathbb{E}[S_{n_{c}}(\eta_{0})]=0$, by the quadratic formula, the two
roots are:

$\begin{split}\eta_{0,1}=&\frac{1}{2n_{o}\left(\bar{Y}-\theta^{*}\right)\left(\bar{Y}{-\theta^{*}}-2n_{e}(\bar{Y}-\bar{X})\right)}\left\\{\sigma^{2}-2n_{e}\left(\bar{Y}^{2}+n_{e}(\bar{Y}^{2}-\bar{X}^{2})-2\left(\bar{Y}+(\bar{Y}-\bar{X})n_{e}\right)\theta^{*}+{\theta^{*}}^{2}\right)\right\\}\\\
&+\sqrt{\sigma^{4}+4(\bar{X}-\bar{Y})^{2}n_{e}^{2}\left(\sigma^{2}+\left(\bar{Y}+n_{e}(\bar{Y}-\bar{X})+\theta^{*}\right)^{2}\right)}\\}\\\
\end{split}$ (22)
$\begin{split}\eta_{0,2}=&\frac{1}{2n_{o}\left(\bar{Y}-\theta^{*}\right)\left(\bar{Y}-\theta^{*}-2n_{e}(\bar{Y}-\bar{X})\right)}\left\\{\sigma^{2}-2n_{e}\left(\bar{Y}^{2}+n_{e}(\bar{Y}^{2}-\bar{X}^{2})-2\left(\bar{Y}+(\bar{Y}-\bar{X})n_{e}\right)\theta^{*}+{\theta^{*}}^{2}\right)\right)\\\
&-\sqrt{\sigma^{4}+4(\bar{X}-\bar{Y})^{2}n_{e}^{2}\left(\sigma^{2}+\left(\bar{Y}+n_{e}(\bar{Y}-\bar{X})+\theta^{*}\right)^{2}\right)}\\}\end{split}$
(23)

###  A.3 Scenario 1: $k\geq 1/2$ or $\delta^{*}=0$

We start with the case of strict inequality, that is, $k>1/2$. For
$X_{i}\sim\mathcal{N}(\theta^{*},\sigma^{2})$ and
$Y_{i}\sim\mathcal{N}(\theta^{*}+\delta^{*}/n^{k},\sigma^{2})$ where $k>1/2$,
that is, $\delta^{*}/n^{k}=o_{p}(\sqrt{n})$, by CLT, we have

$\displaystyle\bar{X}\xrightarrow{d}\mathcal{N}(\theta^{*},\frac{\sigma^{2}}{n_{e}})$
$\displaystyle\bar{Y}\xrightarrow{d}\mathcal{N}(\theta^{*},\frac{\sigma^{2}}{n_{o}})$

so we write

$\displaystyle\bar{X}=\theta^{*}+\frac{\sigma}{\sqrt{n}_{e}}Z_{x}$
$\displaystyle\bar{Y}=\theta^{*}+\frac{\sigma}{\sqrt{n}_{o}}Z_{y},$ (24)

where $Z_{x}$ and $Z_{y}$ are standard normal random variables. This
representation clearly holds if $\delta^{*}=0$. Assuming
$Z_{x}\neq\sqrt{c}Z_{y}$, which is implied by $\bar{X}\neq\bar{Y}$, and
$Z_{y}\neq 0$, substituting (27) in the root expressions (22) and (23) gives

$\begin{split}\eta_{0,1}&=\frac{1}{2Z_{y}\Biggl{(}Z_{y}+2n_{o}\left(-\sqrt{c}Z_{x}+cZ_{y}\right)\Bigr{)}}\Biggl{\\{}1+2cn_{o}Z_{x}^{2}-2c\left(1+cn_{o}\right)Z_{y}^{2}\\\
&+\sqrt{1+4c\left(Z_{x}-\sqrt{c}Z_{y}\right)^{2}\left(n_{o}\left(1+cn_{o}{Z_{x}}^{2}\right)-2\sqrt{c}n_{o}\left(1+cn_{o}\right)Z_{x}Z_{y}+\left(1+cn_{o}\right)^{2}Z_{y}^{2}\right)}\Biggr{\\}}\\\
&=\frac{1}{4n_{o}Z_{y}(cZ_{y}-\sqrt{c}Z_{x})}\left\\{2cn_{o}Z_{x}^{2}-2c^{2}n_{o}Z_{y}^{2}+\sqrt{4c(Z_{x}-\sqrt{c}Z_{y})^{2}\left(cn_{o}^{2}Z_{x}^{2}-2c^{3/2}n_{o}^{2}Z_{x}Z_{y}+c^{2}n_{o}^{2}Z_{y}^{2}+O(n)\right)}\right\\}\\\
&\times\left(1+O\left(\frac{1}{n}\right)\right)\\\
&=\frac{2cn_{o}}{4n_{o}Z_{y}(cZ_{y}-\sqrt{c}Z_{x})}\left\\{Z_{x}^{2}-cZ_{y}^{2}+\sqrt{(Z_{x}-\sqrt{c}Z_{y})^{2}\left(Z_{x}^{2}-2\sqrt{c}Z_{x}Z_{y}+cZ_{y}^{2}+O\left(\frac{1}{n}\right)\right)}\right\\}\left(1+O\left(\frac{1}{n}\right)\right)\\\
&=\frac{\sqrt{c}}{2Z_{y}(\sqrt{c}Z_{y}-Z_{x})}\left\\{Z_{x}^{2}-cZ_{y}^{2}+\sqrt{(Z_{x}-\sqrt{c}Z_{y})^{4}}\right\\}+O\left(\frac{1}{n}\right)\\\
&=\frac{\sqrt{c}\left(Z_{x}^{2}-2\sqrt{c}Z_{x}Z_{y}\right)}{2Z_{y}(\sqrt{c}Z_{y}-Z_{x})}+O\left(\frac{1}{n}\right)\\\
&=\frac{\sqrt{c}}{2Z_{y}(\sqrt{c}Z_{y}-Z_{x})}\left\\{Z_{x}^{2}-cZ_{y}^{2}+(Z_{x}-\sqrt{c}Z_{y})^{2}\right\\}+O\left(\frac{1}{n}\right)\\\
&=-\sqrt{c}\frac{Z_{x}}{Z_{y}}+O\left(\frac{1}{n}\right)\end{split}$ (25)

We observe that $\eta_{0,1}$ asymptotically follows a Cauchy distribution
which suggests that the mean of $\hat{\eta}$ does not converge and hence is
not identifiable.

It is easy to see that, following the same derivation, the other root can be
simplified into

$\begin{split}&\eta_{0,2}=-c+O\left(\frac{1}{n}\right),\end{split}$ (26)

which does not satisfy our requirement that $n_{e}+\eta n_{o}>0$ for (17) to
make sense. It easily follows that the same result holds when $\delta^{*}=0$.

When k = 1/2 and $\delta^{*}\neq 0$, (27) becomes

$\displaystyle\bar{X}=\theta^{*}+\frac{\sigma}{\sqrt{n}_{e}}Z_{x}$
$\displaystyle\bar{Y}=\theta^{*}+\frac{\sigma}{\sqrt{n}_{o}}(Z_{y}+\delta^{*}),.$
(27)

Following the same derivation as (25) gives

$\begin{split}\eta_{0,1}-\sqrt{c}\frac{Z_{x}}{Z_{y}+\delta^{*}}+O\left(\frac{1}{n}\right),\end{split}$
(28)

which suggests that $\eta_{0}$ asymptotically follows the ratio distribution
of two independent normal distributions, with one of them being non-central.
Such distribution is heavy-tailed and has no finite moments (Hinkley, 1969;
Marsaglia, 2006). Therefore, in this case, the mean of $\hat{\eta}$ also does
not converge.

###  A.4 Scenario 1: $0<k<1/2$ or $\delta^{*}\neq 0$

When the bias decreases with sample size at a rate slower than $\sqrt{n}$, the
bias does not vanish as the CLT applies. Then (27) becomes

$\displaystyle\bar{X}=\theta^{*}+\frac{\sigma}{\sqrt{n}_{e}}Z_{x}$
$\displaystyle\bar{Y}=\theta^{*}+\delta(n)+\frac{\sigma}{\sqrt{n}_{o}}Z_{y},$
(29)

Because of the bias term $\delta(n)$, the root expressions (22) and (23) would
become hard to simplify if we substitute (29) directly into them. Instead, we
observe the $\mathbb{E}[S_{n_{c}}(\eta_{0})]$ expression as in (21)

$\displaystyle\mathbb{E}[S_{n_{c}}(\eta_{0})]$
$\displaystyle=\frac{n_{o}\sigma^{2}}{(\sigma^{2}+\hat{\sigma}^{2})^{2}\left(n_{e}+\eta
n_{o}\right)^{4}}\big{\\{}2\left(n_{e}\left(\bar{X}-\theta^{*}\right)+\eta
n_{o}\left(\bar{Y}-\theta^{*}\right)\right)n_{e}\left(\bar{Y}-\bar{X}\right)\left(n_{e}+\eta
n_{o}+1\right)$
$\displaystyle\quad+\left(n_{e}\left(\bar{X}-\theta^{*}\right)+\eta
n_{o}\left(\bar{Y}-\theta^{*}\right)\right)^{2}-\left(n_{e}+\eta
n_{o}\right)\sigma^{2}\big{\\}}$ (30)
$\displaystyle=\frac{2n_{e}n_{o}\sigma^{2}}{(\sigma^{2}+\hat{\sigma}^{2})^{2}\left(n_{e}+\eta
n_{o}\right)^{3}}\left(n_{e}\left(\bar{X}-\theta^{*}\right)+\eta
n_{o}\left(\bar{Y}-\theta^{*}\right)\right)\left(\bar{Y}-\bar{X}\right)+O\left(n^{-1}\right).$
(31)

We can see that asymptotically, $\eta_{0}$ and $\eta_{0,2}$ should be the
solutions to

$\displaystyle n_{e}\left(\bar{X}-\theta^{*}\right)+\eta
n_{o}\left(\bar{Y}-\theta^{*}\right)$ $\displaystyle=0$

Substituting (29) in, we can write out $\eta_{0}$ as

$\displaystyle\eta_{0}=-\frac{c}{\delta}\left(\frac{\sigma}{\sqrt{cn_{o}}}Z_{x}\right)+O(\frac{1}{n}).$

Next, we want to find the Fisher information evaluated at $\eta_{0}$:

$\displaystyle\mathcal{I}\left(\eta_{0}\right)$
$\displaystyle=-\mathbb{E}\left[\frac{\partial}{\partial\eta}S_{n_{c}}(\eta_{0})\right]$
$\displaystyle=-\frac{\partial}{\partial\eta}\mathbb{E}\left[S_{n_{c}}(\eta)\right]\big{\rvert}_{\eta=\eta_{0}}$
$\displaystyle=-\frac{\partial}{\partial\eta}\left[-\frac{1}{2}\frac{2n_{e}n_{o}\sigma^{2}}{(\sigma^{2}+\hat{\sigma}^{2})^{2}\left(n_{e}+\eta
n_{o}\right)^{3}}\left(n_{e}\left(\bar{X}-\theta^{*}\right)+\eta
n_{o}\left(\bar{Y}-\theta^{*}\right)\right)\left(\bar{Y}-\bar{X}\right)+O\left(n^{-1}\right)\middle]\right\rvert_{\eta=\eta_{0}}$
$\displaystyle=\frac{\partial}{\partial\eta}\left[\frac{n_{e}n_{o}\sigma^{2}\left(n_{e}\frac{\sigma}{\sqrt{n}_{e}}Z_{x}+\eta
n_{o}\left(\delta+\frac{\sigma}{\sqrt{n}_{o}}Z_{y}\right)\right)\left(\delta+\frac{\sigma}{\sqrt{n}_{o}}Z_{y}-\frac{\sigma}{\sqrt{n}_{e}}Z_{x}\right)}{(\sigma^{2}+\hat{\sigma}^{2})^{2}\left(n_{e}+\eta
n_{o}\right)^{3}}+O\left(\frac{1}{n}\right)\middle]\right\rvert_{\eta=\eta_{0}}$
$\displaystyle=\frac{\partial}{\partial\eta}\left[\frac{n_{e}n_{o}\sigma^{2}}{(\sigma^{2}+\hat{\sigma}^{2})^{2}\left(n_{e}+\eta
n_{o}\right)^{3}}\eta
n_{o}\delta^{2}+O\left(\frac{1}{\sqrt{n}}\right)\middle]\right\rvert_{\eta=\eta_{0}}$
$\displaystyle=\frac{\partial}{\partial\eta}\left[\frac{n_{e}n_{o}^{2}\sigma^{2}\delta^{2}\eta}{\left(\sigma^{2}+\frac{\sigma^{2}}{n_{e}+\eta
n_{o}}\right)^{2}\left(n_{e}+\eta
n_{o}\right)^{3}}+O\left(\frac{1}{\sqrt{n}}\right)\middle]\right\rvert_{\eta=\eta_{0}}$
$\displaystyle=\frac{\partial}{\partial\eta}\left[\frac{n_{e}n_{o}^{2}\sigma^{2}\delta^{2}\eta}{\sigma^{4}\left(n_{e}+\eta
n_{o}+1\right)^{2}\left(n_{e}+\eta
n_{o}\right)}+O\left(\frac{1}{\sqrt{n}}\right)\middle]\right\rvert_{\eta=\eta_{0}}$
$\displaystyle=\frac{n_{e}n_{o}^{2}\delta^{2}}{\sigma^{2}}\frac{\partial}{\partial\eta}\left[\frac{\eta}{\left(n_{e}+\eta
n_{o}\right)^{3}}\middle]\right\rvert_{\eta=\eta_{0}}+O\left(\frac{1}{\sqrt{n}}\right)$
$\displaystyle=\frac{n_{e}n_{o}^{2}\delta^{2}}{\sigma^{2}}\left[\frac{\left(n_{e}+\eta
n_{o}\right)^{3}-3n_{o}(n_{e}+\eta n_{o})^{2}\eta}{\left(n_{e}+\eta
n_{o}\right)^{6}}\middle]\right\rvert_{\eta=\eta_{0}}+O\left(\frac{1}{\sqrt{n}}\right)$
$\displaystyle=\frac{n_{e}n_{o}^{2}\delta^{2}}{\sigma^{2}}\left[\frac{n_{e}-2n_{o}\eta}{\left(n_{e}+\eta
n_{o}\right)^{4}}\middle]\right\rvert_{\eta=\eta_{0}}+O\left(\frac{1}{\sqrt{n}}\right)$
and using that $\eta_{0}=-\frac{c}{\delta
n_{o}}\left(\frac{\sigma}{\sqrt{n_{e}}}Z_{x}\right)+O(\frac{1}{n})$
$\displaystyle=\frac{n_{e}n_{o}^{2}\delta^{2}}{\sigma^{2}}\frac{n_{e}-2n_{o}\eta_{0}}{\left(n_{e}-n_{e}\left(\frac{\sigma}{\delta\sqrt{n_{e}}}Z_{x}+O(\frac{1}{n})\right)\right)^{4}}+O\left(\frac{1}{\sqrt{n}}\right)$
$\displaystyle=\frac{n_{e}n_{o}^{2}\delta^{2}}{\sigma^{2}}\frac{n_{e}+2n_{o}\left(\frac{n_{e}}{\delta
n_{o}}\left(\frac{\sigma}{\sqrt{n_{e}}}Z_{x}\right)+O(\frac{1}{n})\right)}{n_{e}^{4}}\left(1+\frac{4\sigma}{\delta\sqrt{n_{e}}}Z_{x}+O(\frac{1}{n})\right)+O\left(\frac{1}{\sqrt{n}}\right)$
$\displaystyle=\frac{n_{o}^{2}\delta^{2}}{n_{e}^{2}\sigma^{2}}\left(1+\frac{2\sigma}{\delta\sqrt{n_{e}}}Z_{x}+O(\frac{1}{n})\right)\left(1+\frac{4\sigma}{\delta\sqrt{n_{e}}}Z_{x}+O(\frac{1}{n})\right)+O\left(\frac{1}{\sqrt{n}}\right)$
$\displaystyle=\frac{n_{o}^{2}\delta^{2}}{n_{e}^{2}\sigma^{2}}+O\left(1/{\sqrt{n}}\right).$

Therefore, together with the fact that $n_{e}=cn_{o}$,

$\sqrt{n_{c}}\left(\widehat{\eta}-\eta_{0}\right)\xrightarrow{D}\mathcal{N}\left(0,\frac{\sigma^{2}c^{2}}{\delta^{2}}+O\left(\frac{1}{\sqrt{n}}\right)\right),$

where
$\eta_{0}=-\frac{c}{\delta}\left(\frac{\sigma}{\sqrt{cn_{o}}}Z_{x}\right)+O(\frac{1}{n})$.
Or equivalently, we can write

$w=\widehat{\eta}_{n_{c}}-\eta_{0}=\frac{c}{\delta}\frac{\sigma}{\sqrt{n}_{c}}Z_{w}+O(\frac{1}{n}),$

where $Z_{w}$ is a standard normal variable independent of $Z_{x}$ and
$Z_{y}$. Then it follows that

$\widehat{\eta}=\eta_{0}+w=\frac{c}{\delta}\left(-\frac{\sigma}{\sqrt{cn_{o}}}Z_{x}+\frac{\sigma}{\sqrt{n_{c}}}Z_{w}\right)+O\left(\frac{1}{n}\right).$

We note that $\mathcal{I}\left(\eta_{0}\right)$ is asymptotically positive,
which means the second derivative of $\operatorname{ELPD}(\eta)$ with
respective to $\eta$ is negative in expectation. It confirms that $\eta_{0}$
is indeed a local maxima.

We also observe that as $\delta$ decreases, the probability of $\hat{\eta}$
exceeding 1 increases. This means that if we restrict $\eta$ within $[0,1]$,
it is more likely that we will choose $\eta=1$ as the best value.

Setting $\eta=\widehat{\eta}$ our proposed estimator in (16) becomes

$\displaystyle\hat{\theta}_{\hat{\eta}}$ $\displaystyle=\frac{\frac{\sigma
Z_{x}}{\sqrt{cn_{o}}}+\theta^{*}+\frac{\left(\frac{\sigma
Z_{w}}{\sqrt{n_{c}}}-\frac{\sigma
Z_{x}}{\sqrt{cn_{o}}}+O\left(\frac{1}{n}\right)\right)\left(\delta+\frac{\sigma
Z_{y}}{\sqrt{n_{o}}}+\theta^{*}\right)}{\delta}}{1+\frac{\sigma\left(\frac{Z_{w}}{\sqrt{n_{c}}}-\frac{Z_{x}}{\sqrt{cn_{o}}}+O\left(\frac{1}{n}\right)\right)}{\delta}}$
$\displaystyle=\left(\frac{\sigma
Z_{x}}{\sqrt{cn_{o}}}+\theta^{*}+\frac{\left(\frac{\sigma
Z_{w}}{\sqrt{n_{c}}}-\frac{\sigma
Z_{x}}{\sqrt{cn_{o}}}+O\left(\frac{1}{n}\right)\right)\left(\delta+\frac{\sigma
Z_{y}}{\sqrt{n_{o}}}+\theta^{*}\right)}{\delta}\right)\left(1+O\left(\frac{1}{\sqrt{n}}\right)\right)$
$\displaystyle=\theta^{*}+\frac{\sigma}{\sqrt{cn_{o}}}Z_{x}+\frac{\sigma}{\sqrt{n_{c}}}Z_{w}-\frac{\sigma}{\sqrt{cn_{o}}}Z_{x}+O\left(\frac{1}{n}\right)$
$\displaystyle=\theta^{*}+\frac{\sigma}{\sqrt{n_{c}}}Z_{w}+O\left(\frac{1}{n}\right).$
(32)

From (32) we can see that the proposed estimator
$\widehat{\theta}_{\hat{\eta}}$ asymptotically converges to $\theta^{*}$. As a
result, the MSE is asymptotically,

$\displaystyle\mathcal{R}(\hat{\theta}_{\eta},\theta^{*})$
$\displaystyle=\mathbb{E}_{X,Y}\left(\|\hat{\theta}_{\hat{\eta}}-\theta^{*}\|^{2}\right)=\mathbb{E}\left[\frac{\sigma}{\sqrt{n_{c}}}Z_{w}+O\left(\frac{1}{n}\right)\right]^{2}$
$\displaystyle=\frac{\sigma^{2}}{n_{c}}+O\left(\frac{1}{n^{3/2}}\right).$

Optimistically, if we use a Leave-One-Out Cross Validation approach, it is
close to using a calibration dataset of size $n_{c}=n_{e}$. Therefore, the MSE
of our proposed estimator converges to that using the experimental data alone.
∎

##  Appendix B Normal approximation of posterior

Following the setup outlined in Section 3, we have an experimental dataset of
size $n_{e}$, $X=(X_{1},\dots,X_{n_{e}})$ and an observational dataset of size
$n_{o}$, $Y=(Y_{1},\dots,Y_{n_{o}})$. We assume that they share the same
causal parameters $\theta^{*}$.

We consider the following estimation equations:

$\displaystyle\sum_{i=1}^{n_{e}}f(\theta;X_{i})$ $\displaystyle=0,$ (33)
$\displaystyle\sum_{i=1}^{n_{o}}g(\theta;Y_{i})$ $\displaystyle=0,$ (34)
$\displaystyle\sum_{i=1}^{n_{e}}f(\theta;X_{i})+\eta\sum_{i=1}^{n_{o}}g(\theta;Y_{i})$
$\displaystyle=0,$ (35)

where where $\mathbb{E}f(\theta^{*};X_{i})=0$ and
$\mathbb{E}g(\theta^{*}+\delta;Y_{i})=0$ for some unknown bias $\delta$. From
the standard theory of estimation equations, we know that the solution to
(33), $\hat{\theta}_{e}$, is asymptotically normal following
$\sqrt{n_{e}}(\hat{\theta}_{e}-\theta^{*})\to N(0,V)$, and likewise,
$\sqrt{n_{o}}(\hat{\theta}_{o}-\theta^{*}-\delta)\to N(0,W)$.

Let the solution to (35) be $\hat{\theta}_{\eta}$. We can approximate each sum
in (35) by the Taylor Expansion and use the fact that
$\left[\sum_{i=1}^{N_{e}}f(\hat{\theta}_{e};X_{i})\right]=0$ and
$\left[\sum_{i=1}^{N_{o}}g(\hat{\theta}_{o};X_{i})\right]=0$, then

$\displaystyle n_{e}V^{-1}(\hat{\theta}_{\eta}-\hat{\theta}_{e})+\eta
n_{o}W^{-1}(\hat{\theta}_{\eta}-\hat{\theta}_{o})$ $\displaystyle=0$
$\displaystyle\left(n_{e}V^{-1}+\eta n_{o}W^{-1}\right)\hat{\theta}_{\eta}$
$\displaystyle=n_{e}V^{-1}\hat{\theta}_{e}+\eta n_{o}W^{-1}\hat{\theta}_{o}.$

Then we have an expression for $\hat{\theta}_{\eta}$ as

$\hat{\theta}_{\eta}=\left(n_{e}V^{-1}+\eta
n_{o}W^{-1}\right)^{-1}\left(n_{e}V^{-1}\hat{\theta}_{e}+\eta
n_{o}W^{-1}\hat{\theta}_{o}\right).$

Let $\mathcal{I}(\hat{\theta}_{e})$ and $\mathcal{I}(\hat{\theta}_{o})$ be the
Fisher information matrices for the experimental and observational data
respectively, then the combined Fisher Information matrix is

$\displaystyle\mathcal{I}(\hat{\theta}_{\eta})$
$\displaystyle=-\mathbb{E}\left[\frac{\partial}{\partial\theta}\left(\sum_{i=1}^{n_{e}}f(\theta;X_{i})+\eta\sum_{i=1}^{n_{o}}g(\theta;Y_{i})\right)\right]$
$\displaystyle=-\mathbb{E}\left[\frac{\partial}{\partial\theta}\sum_{i=1}^{n_{e}}f(\theta;X_{i})\right]-\eta\,\mathbb{E}\left[\frac{\partial}{\partial\theta}\sum_{i=1}^{n_{o}}g(\theta;Y_{i})\right]$
$\displaystyle=\mathcal{I}(\hat{\theta}_{e})+\eta\,\mathcal{I}(\hat{\theta}_{o}).$

That is, the combined Fisher Information is a weighted sum of the two
datasets. In the simulations presented in Section 5, we use the sandwich
estimators of $V$ and $W$ and Fisher information calculated using the R
package causl (Evans, 2021). Then, as outlined in Section 4.4, we can
approximate the $\eta$-powered posterior with

$\pi_{\eta}\left(\theta\,|\,X,Y\right)\approx\mathcal{N}\left(\hat{\theta}_{\eta},[\mathcal{I}(\hat{\theta}_{\eta})]^{-1}\right),$

where $\hat{\theta}_{\eta}$ and $\mathcal{I}(\hat{\theta}_{\eta})$ can be
approximated as derived above.

##  Appendix C Details of simulation

###  C.1 Implementation of other methods

In this section, we provide details on how we implement the comparator methods
in the simulation presented in Section 5.

The shrinkage estimators in Green et al. (2005) require more than two strata
and those in Rosenman et al. (2020) require as least four strata to guarantee
a risk reduction, therefore we split both datasets into 10 subgroups. In the
main simulation in Section 5.2, the stratification is based on deciles of $C$
in the observational data; for the simulation presented in Section 5.4, there
are more covariates and so we use stratification based on the two levels of
$C_{5}$ and quintiles of $C_{1}$. We use $\hat{\bm{\theta}}_{\bm{o}}$,
$\hat{\bm{\theta}}_{\bm{e}}\in\mathbb{R}^{K}$ where $K=10$ to denote the
vector of CATE estimates in the 10 strata in $\mathcal{D}_{e}$ and
$\mathcal{D}_{o}$ respectively.

####  C.1.1 Green et al. (2005)

We implemented the estimators proposed in Green et al. (2005), which
generalizes the estimators in Green and Strawderman (1991) to allow for
heteroscedastic variances. Specifically, we used the estimators

$\bm{\delta}_{1}^{+}(\hat{\bm{\theta}}_{\bm{e}},\hat{\bm{\theta}}_{\bm{o}})=\hat{\bm{\theta}}_{\bm{o}}+\left(1-\frac{a}{\left(\hat{\bm{\theta}}_{\bm{o}}-\hat{\bm{\theta}}_{\bm{e}}\right)^{T}\bm{\Sigma}_{e}^{-1}\left(\hat{\bm{\theta}}_{\bm{o}}-\hat{\bm{\theta}}_{\bm{e}}\right)}\right)_{+}\left(\hat{\bm{\theta}}_{\bm{o}}-\hat{\bm{\theta}}_{\bm{e}}\right)$

and

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

where $\Sigma_{e}$ is the variance-covariance matrix of
$\bm{\hat{}}\theta_{e}$; note, of course, that we assume the components of
$\bm{\hat{}}\theta_{e}$ (and $\bm{\hat{}}\theta_{o}$) are mutually
independent. We have also adopted the default choice of $a=K-2$, where $K=10$
in our simulation setup. $\Sigma_{e}$ is a diagonal matrix with the variance
of each component on its diagonal. In the simulation, we implemented both
estimators and they gave very similar results. For better readability, we only
plot $\bm{\delta}_{1}^{+}$ in Figure 5, 6 and 7.

####  C.1.2 Rosenman et al. (2020)

We implemented two estimators proposed in Rosenman et al. (2020). The first
one is the shrinkage estimator which shares a common shrinkage factor across
all components of $\bm{\hat{}}\theta_{e}$ and $\bm{\hat{}}\theta_{o}$.

$\bm{\kappa}_{1+}=\hat{\bm{\theta}}_{\bm{o}}+\left(1-\frac{\operatorname{Tr}\left(\bm{\Sigma}_{e}\bm{D}\right)}{\left(\hat{\bm{\theta}}_{\bm{o}}-\hat{\bm{\theta}}_{\bm{e}}\right)^{T}\bm{D}\left(\hat{\theta}_{\bm{o}}-\hat{\bm{\theta}}_{\bm{e}}\right)}\right)_{+}\left(\hat{\bm{\theta}}_{\bm{o}}-\hat{\bm{\theta}}_{\bm{e}}\right)$

where $\bm{D}$ is a diagonal weight matrix representing the relative
importance of the strata estimates. In our implementation, we set
$\bm{D}=\mathrm{I}_{K}$. Another estimator they proposed uses variance-
weighted shrinkage factors:

$\bm{\kappa}_{2+}=\hat{\bm{\theta}}_{\bm{o}}+\left(\bm{I}_{K}-\frac{\operatorname{Tr}\left(\bm{\Sigma}_{e}^{2}\bm{D}\right)\bm{\Sigma}_{e}}{\left(\hat{\bm{\theta}}_{\bm{o}}-\hat{\bm{\theta}}_{\bm{e}}\right)^{T}\bm{\Sigma_{e}}^{2}\bm{D}\left(\hat{\theta}_{o}-\hat{\theta}_{e}\right)}\right)_{+}\left(\hat{\bm{\theta}}_{\bm{o}}-\hat{\bm{\theta}}_{\bm{e}}\right).$

In the simulation, we implemented both estimators and they gave very similar
results. We only plot $\bm{\kappa}_{1+}$ in Figure 5, 6 and 7 for simplicity.

For all the abovementioned shrinkage estimators, we use the IPW estimates for
$\hat{\bm{\theta}}_{\bm{o}}$ and $\hat{\bm{\theta}}_{\bm{e}}$:

$\displaystyle\hat{\theta}_{e,k}$
$\displaystyle=\sum_{i\in\mathcal{E}_{k}}\frac{T_{i}Y_{i}}{e_{e}(Z_{i},C_{i})}-\frac{(1-T_{i})Y_{i}}{1-e_{e}(Z_{i},C_{i})}$
$\displaystyle\text{ for }k=1,2,\dots,10$ $\displaystyle\hat{\theta}_{o,k}$
$\displaystyle=\sum_{i\in\mathcal{O}_{k}}\frac{T_{i}Y_{i}}{e_{o}(Z_{i},C_{i})}-\frac{(1-T_{i})Y_{i}}{1-e_{o}(Z_{i},C_{i})}$
$\displaystyle\text{ for }k=1,2,\dots,10$

where $e_{e}(Z_{i},C_{i})$ and $e_{o}(Z_{i},C_{i})$ are the propensity scores,
and $\mathcal{E}_{k}$ and $\mathcal{O}_{k}$ denote the observations in the
$k$-th stratum of the randomized and observational data, respectively.

####  C.1.3 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}(\bm{z},\bm{c})=\mathbb{E}[Y\,|\,T=1,\bm{Z}=\bm{z},\bm{C}=\bm{c},D=\mathcal{D}_{o}]$
and
$m^{0}(\bm{z},\bm{c})=\mathbb{E}[Y\,|\,T=0,\bm{Z}=\bm{z},\bm{C}=\bm{c},D=\mathcal{D}_{o}]$,
and calculate their difference as an estimate of CATE estimate
$\hat{\omega}(\bm{z},\bm{c})=\hat{m}^{1}(\bm{z},\bm{c})-\hat{m}^{0}(\bm{z},\bm{c})$.

  2. Step 2:

We learn the bias correction function, assuming linearity,
$\xi(\bm{z},\bm{c})=\beta^{T}\bm{v}$ where $\bm{v}$ is the combined vector of
$(\bm{z},\bm{c})^{T}$ and

$\hat{\beta}=\underset{\beta}{\arg\min}\sum_{i=1}^{n_{e}}\left(q_{i}Y_{i}-\hat{\omega}\left(\bm{V}_{i}\right)-\beta^{T}\bm{V}_{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 strata CATE estimates as

$\hat{\tau_{k}}=\frac{1}{n_{e,k}}\sum_{i\in\mathcal{E}_{k}}\hat{\omega}(\bm{V}_{i})+\hat{\beta}^{T}\bm{V}_{i}\text{
for }k=1,2,\dots,10,$

where $\mathcal{E}_{k}$ denotes the collection of randomized units in stratum
$k$. Specifically, with respect to $m^{1}$, $m^{0}$ and $\xi$, we fit the
models using the formula $Y\sim C+Z$ for the simulation in Section 5.2 and
$Y\sim C_{1}+C_{2}+C_{3}+C_{4}+C_{5}+Z_{1}+Z_{2}$ for the simulation in
Section 5.4. It is worth noting that these formulae correctly specify the
relationships according to the data-generating processes.

####  C.1.4 Oberst et al. (2022)

We construct the estimator proposed in Oberst et al. (2022) as

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

where $\hat{\tau}_{e}$ and $\hat{\tau}_{o}$ are ATE estimates from
$\mathcal{D}_{e}$ and $\mathcal{D}_{o}$ respectively. As we assume that
$\hat{\tau}_{e}$ and $\hat{\tau}_{o}$ are independent, the covariance term in
the original formula of $\hat{\lambda}$ drops out. We use the ATE estimates
and their variance inferred from the fully parametric model as
$\hat{\tau}_{e}$ and $\hat{\tau}_{o}$ so that $\hat{\tau}_{\hat{\lambda}}$ is
as comparable to our approach as possible.

###  C.2 Setup of the simulation in 5.4

Very similar to the setup of the simulation in 5.2, we model the distributions
as follows555For simplicity, we use the shorthand notation $\bm{C}_{1:p}$ to
denote the set $\\{C_{1},\cdots,C_{p}\\}$ :

$\displaystyle\bm{C}_{1:3}$ $\displaystyle\sim N(0,1)$
$\displaystyle\bm{C}_{4:5}$ $\displaystyle\sim\operatorname{Bernoulli}(0.5)$
$\displaystyle Z_{1}\mid C_{1}$ $\displaystyle\sim N(\mu_{z_{1}},\,1)$
$\displaystyle Z_{2}\mid C_{4}$ $\displaystyle\sim N(\mu_{z_{2}},1)$
$\displaystyle T\,|\,C_{1},C_{5},Z_{1}$
$\displaystyle\sim\operatorname{Bernoulli}(\mu_{t})$ $\displaystyle
Y(t)\mid\bm{C}_{1:5}$ $\displaystyle\sim N(\mu_{y},\,1).$

where,

$\displaystyle\mu_{z_{1}}$ $\displaystyle=C_{1}$ $\displaystyle\mu_{z_{2}}$
$\displaystyle=C_{4}$ $\displaystyle\operatorname{logit}\mu_{t}$
$\displaystyle=0.5+0.1\,C_{1}+0.6\,Z_{1}+0.4\,C_{5}\,Z+0.1\,C_{1}Z_{1}+\psi U$
$\displaystyle\mu_{y}$
$\displaystyle=\left[1,C_{1},C_{2},C_{3},C_{4},C_{5}\right]\left[\begin{array}[]{cc}0.1&0.7\\\
0.2&0.8\\\ 0.3&0.9\\\ 0.4&1.0\\\ 0.5&1.1\\\ 0.6&1.2\\\
\end{array}\right]\left[\begin{array}[]{c}1\\\ T\end{array}\right]+\psi U$

Additionally, let the dependence structure between $Y$ and $\bm{Z}_{1:2}$
given $\\{T=t,\bm{C}_{1:5}=\bm{c}_{1:5}\\}$, $\phi^{*}_{Y\bm{Z}\mid T\bm{C}}$,
be a conditionally trivariate Gaussian copula, with correlation matrix

$\displaystyle\bm{\Sigma}=\begin{pmatrix}1&\rho_{Z_{1}Z_{2}}&\rho_{Z_{1}Y}\\\
\rho_{Z_{1}Z_{2}}&1&\rho_{Z_{2}Y}\\\
\rho_{Z_{1}Y}&\rho_{Z_{2}Y}&1\end{pmatrix},$

where
$\rho_{z_{1}z_{2}}=\rho_{z_{1}Y}=\rho_{z_{2}Y}=2\operatorname{expit}(1+t)-1$.
Specifically, the correlations are $0.76$ when $T=1$, and $0.46$ when $T=0$.

We use the same parameterization for the experimental data $\mathcal{D}_{e}$
except that we assume that treatment is randomly assigned and hence replace
the mean of $T$ with $\mu_{t}=0.5$. We set the sample sizes to 2,500 and 250
for the observational and experimental data, respectively. The results in
Section 5.4 are averaged across 500 sets of synthetic datasets.

## References

  * Abadie and Imbens (2016) A. Abadie and G. W. Imbens.  Matching on the estimated propensity score.  _Econometrica_ , 84:781–807, 2016.  ISSN 0012-9682.  doi: 10.3982/ecta11293.
  * 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.
  * Bang and Robins (2005) H. Bang and J. M. Robins.  Doubly robust estimation in missing data and causal inference models.  _Biometrics_ , 61:962–973, 2005.  ISSN 15410420.  doi: 10.1111/j.1541-0420.2005.00377.x.
  * 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.
  * Ben-Michael et al. (2021) E. Ben-Michael, D. J. Greiner, K. Imai, and Z. Jiang.  Safe policy learning through extrapolation: Application to pre-trial risk assessment.  _arXiv preprint arXiv:2109.11679_ , 2021.
  * Bernardo and Smith (1994) J. M. Bernardo and A. F. Smith.  _Bayesian Theory_ , volume 405.  John Wiley & Sons, 1994.
  * Bissiri et al. (2016) P. G. Bissiri, C. C. Holmes, and S. G. Walker.  A general framework for updating belief distributions.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 78(5):1103–1130, 2016.
  * Box (1976) G. E. Box.  Science and statistics.  _Journal of the American Statistical Association_ , 71(356):791–799, 1976.
  * Campbell (1957) D. T. Campbell.  Factors relevant to the validity of experiments in social settings.  _Psychological bulletin_ , 54(4):297, 1957.
  * Cao et al. (2009) W. Cao, A. A. Tsiatis, and M. Davidian.  Improving efficiency and robustness of the doubly robust estimator for a population mean with incomplete data.  _Biometrika_ , 96:723–734, 9 2009.  ISSN 00063444.  doi: 10.1093/biomet/asp033.
  * Carmona and Nicholls (2020) C. Carmona and G. Nicholls.  Semi-modular inference: enhanced learning in multi-modular models by tempering the influence of components.  In _International Conference on Artificial Intelligence and Statistics_ , pages 4226–4235. PMLR, 2020.
  * Cheng and Cai (2021) D. Cheng and T. Cai.  Adaptive combination of randomized and observational data.  _arXiv preprint arXiv:2111.15012_ , 2021.
  * Claeskens and Hjort (2003) G. Claeskens and N. L. Hjort.  The focused information criterion.  _Journal of the American Statistical Association_ , 98(464):900–916, 2003.
  * Dang et al. (2022) 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_ , 2022.
  * Dawid (2021) A. P. Dawid.  Decision-theoretic foundations for statistical causality.  _Journal of Causal Inference_ , 9(1):39–77, 2021\.
  * European Medicines Agency (2021) European Medicines Agency.  A vision for use of real-world evidence in eu medicines regulation.  https://www.ema.europa.eu/en/news/vision-use-real-world-evidence-eu-medicines-regulation, 2021\.  Accessed: 2023-03-13.
  * 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 (2021) R. J. Evans.  _causl: Methods for Specifying, Simulating from and Fitting Causal Models_ , 2021.  R package version 0.4.2.9001.
  * Evans and Didelez (2021) R. J. Evans and V. Didelez.  Parameterizing and simulating from causal models.  _arXiv preprint arXiv:2109.03694_ , 2021.
  * Gelman et al. (1995) A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin.  _Bayesian data analysis_.  Chapman and Hall/CRC, 1995.
  * 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.
  * Grünwald and Van Ommen (2017) P. Grünwald and T. Van Ommen.  Inconsistency of Bayesian inference for misspecified linear models, and a proposal for repairing it.  _Bayesian Analysis_ , 12(4):1069–1103, 2017.
  * Hansen et al. (2021) C. H. Hansen, D. Michlmayr, S. M. Gubbels, K. Mølbak, and S. Ethelberg.  Assessment of protection against reinfection with SARS-CoV-2 among 4 million PCR-tested individuals in Denmark in 2020: a population-level observational study.  _The lancet_ , 397(10280):1204–1212, 2021.
  * Hartman et al. (2015) E. Hartman, R. Grieve, R. Ramsahai, and J. S. Sekhon.  From sample average treatment effect to population average treatment effect on the treated: combining experimental with observational studies to estimate population treatment effects.  _Journal of the Royal Statistical Society: Series A (Statistics in Society)_ , 178(3):757–778, 2015.
  * Hinkley (1969) D. V. Hinkley.  On the ratio of two correlated normal random variables.  _Biometrika_ , 56(3):635–639, 1969.
  * Hirano et al. (2003) K. Hirano, G. W. Imbens, and G. Ridder.  Efficient estimation of average treatment effects using the estimated propensity score.  _Econometrica_ , 71(4):1161–1189, 2003.
  * 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\.
  * Holmes and Walker (2017) C. C. Holmes and S. G. Walker.  Assigning a value to a power likelihood in a general Bayesian model.  _Biometrika_ , 104(2):497–503, 2017.
  * Ibrahim and Chen (2000) J. G. Ibrahim and M.-H. Chen.  Power prior distributions for regression models.  _Statistical Science_ , pages 46–60, 2000.
  * 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 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\.
  * Lyddon et al. (2019) S. P. Lyddon, C. Holmes, and S. Walker.  General Bayesian updating and the loss-likelihood bootstrap.  _Biometrika_ , 106(2):465–478, 2019.
  * Marsaglia (2006) G. Marsaglia.  Ratios of normal variables.  _Journal of Statistical Software_ , 16:1–10, 2006.
  * 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.
  * Peysakhovich and Lada (2016) A. Peysakhovich and A. Lada.  Combining observational and experimental data to find heterogeneous treatment effects.  _arXiv preprint arXiv:1611.02385_ , 2016.
  * Pottegård et al. (2020) A. Pottegård, X. Kurz, N. Moore, C. F. Christiansen, and O. Klungel.  Considerations for pharmacoepidemiological analyses in the SARS-CoV-2 pandemic.  _Pharmacoepidemiology and Drug Safety_ , 29(8):825–831, 2020.
  * 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.
  * 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.
  * Rosenbaum and Rubin (1983) P. R. Rosenbaum and D. B. Rubin.  Assessing sensitivity to an unobserved binary covariate in an observational study with binary outcome.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 45(2):212–218, 1983.
  * Rosenman et al. (2020) E. Rosenman, G. Basse, A. Owen, and M. Baiocchi.  Combining observational and experimental datasets using shrinkage estimators.  _arxiv preprint 2002.06708_ , 2020.
  * Rosenman et al. (2022) E. T. Rosenman, A. B. Owen, M. Baiocchi, and H. R. Banack.  Propensity score methods for merging observational and experimental datasets.  _Statistics in Medicine_ , 41(1):65–86, 2022\.
  * Rubin (1973) D. B. Rubin.  Matching to remove bias in observational studies.  _Biometrics_ , pages 159–183, 1973.
  * Rubin (1974) D. B. Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of Educational Psychology_ , 66(5):688, 1974.
  * 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.
  * 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\.
  * Syring and Martin (2019) N. Syring and R. Martin.  Calibrating general posterior credible regions.  _Biometrika_ , 106(2):479–486, 2019.
  * 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 (2016) U.S Food and Drug Administration.  FDA approves first treatment for rare disease in patients who receive stem cell transplant from blood or bone marrow.  https://www.fda.gov/news-events/press-announcements/fda-approves-first-treatment-rare-disease-patients-who-receive-stem-cell-transplant-blood-or-bone, 2016\.  Accessed: 2023-03-13.
  * U.S Food and Drug Administration (2018a) U.S Food and Drug Administration.  FDA approves new treatment for certain digestive tract cancers.  https://www.fda.gov/news-events/press-announcements/fda-approves-new-treatment-certain-digestive-tract-cancers, 2018a.  Accessed: 2023-03-13.
  * U.S Food and Drug Administration (2018b) U.S Food and Drug Administration.  Framework for FDA’s real-world evidence program.  https://www.fda.gov/media/120060/download, 2018b.  Accessed: 2023-03-13.
  * 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\.
  * Vehtari et al. (2017) A. Vehtari, A. Gelman, and J. Gabry.  Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.  _Statistics and computing_ , 27(5):1413–1432, 2017.
  * 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\.
  * Walker and Hjort (2001) S. Walker and N. L. Hjort.  On Bayesian consistency.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 63(4):811–821, 2001.
  * Watanabe (2010) S. Watanabe.  Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory.  _Journal of machine learning research_ , 11(12), 2010.
  * Wedam et al. (2020) S. Wedam, L. Fashoyin-Aje, E. Bloomquist, S. Tang, R. Sridhara, K. B. Goldberg, M. R. Theoret, L. Amiri-Kordestani, R. Pazdur, and J. A. Beaver.  FDA approval summary: Palbociclib for male patients with metastatic breast cancer.  _Clinical Cancer Research_ , 26(6):1208–1212, 2020.
  * Word et al. (1990) E. Word, J. Johnston, H. P. Bain, B. D. Fulton, J. B. Zaharias, C. M. Achilles, M. N. Lintz, J. Folger, and C. Breda.  The state of Tennessee’s Student/Teacher Achievement Ratio (STAR) project.  _Tennessee Board of Education_ , 1990.
  * Wu and Martin (2022) P. Wu and R. Martin.  A comparison of learning rate selection methods in generalized Bayesian inference.  _Bayesian Analysis_ , 1(1):1–28, 2022.
  * 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 et al. (2020) S. Yang, D. Zeng, and X. Wang.  Elastic integrative analysis of randomized trial and real-world data for treatment heterogeneity estimation.  _arXiv preprint arXiv:2005.10579_ , 2020.
  * Zellner (1988) A. Zellner.  Optimal information processing and Bayes’s theorem.  _The American Statistician_ , 42(4):278–280, 1988\.
  * Zhang and Liang (2011) X. Zhang and H. Liang.  Focused information criterion and model averaging for generalized additive partial linear models.  _The Annals of Statistics_ , pages 174–200, 2011.