# A Bayesian Approach for Selecting Relevant External Data (BASE): Application
to a study of Long-Term Outcomes in a Hemophilia Gene Therapy Trial

Tianyu Pan  Department of Pathology, Stanford University  Xiang Zhang  CSL
Behring  Weining Shen  Department of Statistics, University of California,
Irvine  Ting Ye  Department of Biostatistics, University of Washington

###### Abstract

Gene therapies aim to address the root causes of diseases, particularly those
stemming from rare genetic defects that can be life-threatening or severely
debilitating. While there has been notable progress in the development of gene
therapies in recent years, understanding their long-term effectiveness remains
challenging due to a lack of data on long-term outcomes, especially during the
early stages of their introduction to the market. To address the critical
question of estimating long-term efficacy without waiting for the completion
of lengthy clinical trials, we propose a novel Bayesian framework. This
framework selects pertinent data from external sources, often early-phase
clinical trials with more comprehensive longitudinal efficacy data that could
lead to an improved inference of the long-term efficacy outcome. We apply this
methodology to predict the long-term factor IX (FIX) levels of HEMGENIX
(etranacogene dezaparvovec), the first FDA-approved gene therapy to treat
adults with severe Hemophilia B, in a phase 3 study. Our application showcases
the capability of the framework to estimate the 5-year FIX levels following
HEMGENIX therapy, demonstrating sustained FIX levels induced by HEMGENIX
infusion. Additionally, we provide theoretical insights into the methodology
by establishing its posterior convergence properties.

Keywords: Bayesian model averaging; Data integration; Data fusion; Dynamic
borrowing; Long-term effectiveness

##  1 Introduction

###  1.1 Gene therapy for Hemophilia B

Gene therapy holds great promise as a one-time treatment for life-threatening,
severe-debilitating diseases such as hemophilia, with demonstrated increases
in factor level expression post treatments and substantial reductions in both
bleeds and utilization of factor replacement therapy in treating breakthrough
bleeding (Nathwani et al., 2011, 2014; George et al., 2017; Leebeek and
Miesbach, 2021). However, clinical trial participants have demonstrated
variable expression in factor levels. This uncertainty has generated
significant interest in understanding the long-term efficacy of gene therapy.
Given the high cost of gene therapy, the potential for waning efficacy is
concerning to payers who base the cost on the assumption of long-lasting
effects (Kee and Maio, 2019).

Our paper aims to study the long-term efficacy of HEMGENIX, which is the first
FDA-approved gene therapy designed for treating adults with severe Hemophilia
B. A flourishing body of assessment reports from Health Technology Assessment
(HTA) entities contribute to our understanding, including the Institute for
Clinical and Economic Review (Tice et al., 2022) in the United States and the
National Institute for Health and Care Excellence (Farmer et al., 2023) in the
United Kingdom. The main challenge is that the data on long-term effectiveness
are often lacking during the early marketing period of gene therapy,
necessitating the development of appropriate methods to infer long-term
effectiveness and bridge the evidence gap.

###  1.2 The HOPE-B Study

HOPE-B is a phase 3 study of HEMGENIX, the first FDA-approved gene therapy for
treating Hemophilia B in adults (Pipe et al., 2023, ClinicalTrials.gov number,
NCT03569891). The study involved the enrollment of 54 male participants. At
the point of our data analysis, information spanning up to 36 months post-
treatment was accessible. We opted to utilize only the data from the initial
24 months, reserving the remaining dataset for the purpose of assessing
predictions. The ongoing collection of data will continue for a duration of up
to 15 years as the study progresses. Due to the limited data available, both
from the current study and external sources, predicting the long-term effect
of the gene therapy, such as its effects after 5 years, presents a significant
challenge.

In traditional clinical drug development, multiple early-phase studies (phase
1 and phase 2) are conducted before moving to a larger-scale confirmatory
phase 3 study (Valentinuzzi, 2004). As a result, early-phase studies often
yield data with longer observational periods. In our case, at the time of
analysis, a phase 1 study of 10 patients with 5 years of follow-up (Miesbach
et al., 2021) and a phase 2b study of 3 patients with 2.5 years of follow-up
(Gomez et al., 2021) are crucial additional resources for inferring the long-
term effect. However, it is important to note that the gene vehicle is
different between the phase 1 and phase 2b/3 studies. Specifically, the phase
2b/3 studies employed the Padua gene, while the phase 1 study used a wild-type
factor ix@ gene. These disparities could introduce systematic differences in
the factor level across studies. Indeed, as shown in Figure 1(a), the phase 1
study exhibits noteworthy variations in factor levels and their longitudinal
trajectories in comparison to the other two studies. This highlights the
potential for biased estimates if data across the three studies are naively
pooled. Throughout the rest of the paper, we designate the combined phase 2b
and phase 3 study (which share substantial similarities) as our _internal
study_ , for which we aim to infer the long-term effect, and we consider the
phase 1 study as our _external study_ , from which we seek to borrow
information.

Motivated by this example, a data integration method is necessary for
borrowing information from the external studies (the phase 1 study) to help
infer the long-term efficacy of the internal study (the phase 2b and 3
studies). Moreover, the data integration method should be designed to address
the difference between the internal and external studies and guard against
introducing bias into the inference of the internal studies. To achieve these
two goals, our idea consists of two key steps. First, we perform a
z-transformation on both the external and internal studies, utilizing the
respective mean values and standard deviations computed from the initial 10
weeks of each study. The resulting data are plotted in Figure 1(b). Following
the z-transformation, the external and internal data exhibit greater
similarity. Second, we augment the internal data using a selected subset of
external data that exhibits similar early-stage trajectories to the internal
data in order to gain efficiency of inferring long-term outcomes while
guarding against introducing bias. This approach relies on the key assumption
that similar initial trajectories can predict similar long-term outcomes.
Furthermore, the proportion of the selected external subset offers an
additional interpretation that describes the similarity between the two
datasets. Specifically, a high proportion of the selected external subset
would suggest greater proximity between the two datasets. In summary, the goal
of the data selection procedure is to filter out irrelevant data from the
external studies whilst selecting a subset that aligns with the underlying
data-generating process of the internal trajectories.

Figure 1:  (a) Longitudinal outcome trajectories observed for each study. (b)
The longitudinal trajectories up to the 2nd year after within-study
z-transformation.

###  1.3 Prior work and our contributions

As mentioned in Section 1.2, the goal of combining the external and internal
datasets can be perceived as a data integration problem or, more generally, as
an effort to combine information from multiple data sources for an integrative
and efficient inference.

Over the past decades, data integration has been an active research area.
Various concepts and methods have been developed. To list a few, data fusion
methods combine multiple datasets by assuming a shared latent variable (Liu et
al., 2022) or sufficient conditional overlapping support (Li and Luedtke,
2021) for various purposes such as integrated epigenetic index estimation or
average treatment effect estimation. When datasets are in matrix format,
combining them can sometimes lead to block-wise missingness. To tackle this,
data integration methods consider low-rank matrix recovery (Cai et al., 2016),
spectral clustering (Park et al., 2021), and multiple block-wise imputation
(Xue and Qu, 2021). Multimodal data analysis, where the collected data come in
different types, can also be treated as a form of data selection. For example,
to regress the outcome on the covariates from multiple datasets, people have
considered the conventional linear regression (Li and Li, 2022) and the non-
linear regression (Dai and Li, 2022). From a philosophical standpoint,
Bayesian methods are natural data integration approaches. Among these, the
power prior (Chen et al., 1999, 2000) is a notable example designed for data
integration. It incorporates an uncertain discounting factor $\alpha\in[0,1]$
into the historical likelihood to downgrade its importance while combining it
with the current likelihood (Neuenschwander et al., 2009; Ibrahim et al.,
2015). However, these methods downweight external data equally and do not
effectively address the task of selecting relevant information from external
sources. Numerous other data integration methods have emerged in diverse
directions, including methods for heterogeneous treatment effect estimation
(Yang et al., 2020), long-term treatment effect estimation (Athey et al.,
2019; Imbens et al., 2022), doubly robust estimation for non-probability
samples (Yang et al., 2020; Chen et al., 2020), and adaptively weighted
combination strategy (Chen et al., 2021; Oberst et al., 2022). We refrain from
an exhaustive discussion on this rapidly evolving literature and refer
interested readers to the following reviews (Ritchie et al., 2015; Ibrahim et
al., 2015; Hassler et al., 2023).

Despite the comprehensive developments in data integration, all the
aforementioned methods either rely on incorporating the entire external
information or apply weight scaling to the whole external data to alleviate
the inferential bias when conducting the data integration. However,
considering the non-negligible between-study heterogeneity observed in Figure
1(a), directly combining all external data poses a risk of introducing bias
into early-stage inference. As a remedy, we advocate for a data selection
procedure that selects external subsets that are similar to the internal data
generating mechanism instead of integrating the entire external data. To the
best of our knowledge, the selection of relevant external subsets still
remains an area under investigation. In this article, we propose a novel
Bayesian data selection procedure by assigning a prior to all possible
external data subsets and using the marginal likelihood as the criteria to
favor specific external subsets in the sampling process.

The main contributions of this article are three-fold. First, we propose a
novel data integration procedure, termed data selection, which has been
largely overlooked till now. Second, we provide theoretical justifications for
our method by showing that the relevant external subset can be consistently
selected with high posterior probability. Third, we apply our proposed method
to the hemophilia study, producing promising results that are both
statistically valid and scientifically interpretable. In comparison to the
results given by the direct combination of external and internal data, and the
method without external information incorporation (Shah et al., 2023), our
method offers more reliable and robust conclusions about the long-term
efficacy by selecting relevant external subset.

The remainder of this paper is outlined as follows. We present our methods in
general terms in Section 2, and in the context of the hemophilia study in
Section 3. Results from simulation studies are in Section 4. We return to the
hemophilia application in Section 5. We conclude with a discussion in Section
6. Technical proofs and pseudo-code of MCMC sampling are given in the
Supplementary File.

##  2 A general Bayesian strategy for data selection

###  2.1 Method

In this section, we present a general Bayesian approach for selecting relevant
external data (BASE), i.e., we propose to select subsets from the external
data that are similar to the internal data at the early stage to help infer
the long-term efficacy. We will use the hemophilia study as an example to
concretely describe our idea.

Suppose $Y_{sit}$ is the outcome (i.e., factor level in our application)
measured for patient $i$ at time $t$ in study $s$. The data collected from the
studies are then expressed as
$\\{Y_{sit},i=1,\ldots,N_{s},t\in\mathcal{T}_{s},s=0,1\\}$, where $N_{s}$ is
the number of patients under study $s$, $\mathcal{T}_{s}$ is the available
time points for study $s$, and $s=1$ indicates the internal study and $s=0$
indicates the external study. For a model
$L(\cdot\mid\boldsymbol{\theta}_{s})$ parameterized by
$\boldsymbol{\theta}_{s}$, our data selection procedure proceeds by first
putting a prior $\Pi(\mathcal{C})$ on all the possible subsets of the external
index set denoted as $\mathcal{C}\subseteq\\{1,\ldots,N_{0}\\}$, and then
selecting the external subsets by drawing posterior samples of $\mathcal{C}$
given the outcomes observed before a pre-specified time $T^{*}$, namely,

$\displaystyle\Pi\left(\mathcal{C}\mid\\{Y_{sit},i=1,\ldots,N_{s},t\in\mathcal{T}_{s},t\leq
T^{*},s=0,1\\}\right)$ (1)
$\displaystyle\propto\underbrace{\left(\int\prod_{i=1}^{N_{1}}L(\\{Y_{1it}\\}_{t\in\mathcal{T}_{1},t\leq
T^{*}}\mid\boldsymbol{\theta}_{1})\times\underbrace{\pi(\boldsymbol{\theta}_{1}\mid\\{Y_{0it},i\in\mathcal{C},t\in\mathcal{T}_{0},t\leq
T^{*}\\})}_{\text{the posterior density of $\boldsymbol{\theta}_{1}$ given an
external subset
$\mathcal{C}$}}d\boldsymbol{\theta}_{1}\right)}_{\text{marginal
likelihood}}\times\underbrace{\Pi(\mathcal{C})}_{\text{the prior of
$\mathcal{C}$}},$

where $\pi(\boldsymbol{\theta}_{1})$ and
$\pi(\boldsymbol{\theta}_{1}\mid\cdot)$ denote the prior and posterior density
functions of $\boldsymbol{\theta}_{1}$, and
$L(\\{Y_{1it}\\}_{t\in\mathcal{T}_{1},t\leq
T^{*}}\mid\boldsymbol{\theta}_{1})$ refers to the joint likelihood of the
internal data on the time grid preceding $T^{*}$, parameterized by
$\boldsymbol{\theta}_{1}$. In our application, we designate $T^{*}$ as 2
years, which is the end of the internal study follow-up.

We provide further interpretations of (1) to illustrate the rationale behind
this selection procedure. On the right-hand side of (1), $\Pi(\mathcal{C})$
denotes our prior guess on external subsets $\mathcal{C}$’s. In practice, in
the absence of this information, we can choose, for example, a uniform prior
over all possible $\mathcal{C}$’s. The remaining term on the right-hand side
of (1) is the marginal likelihood (unlike the original definition of the
marginal likelihood where the integration is taken over the prior of
$\boldsymbol{\theta_{1}}$, our definition is more general as it considers the
posterior of $\boldsymbol{\theta_{1}}$ given the external subset
$\mathcal{C}$; and includes the original definition as a special case if
$\mathcal{C}$ is the empty set) that aims to evaluate the similarity between
the generating mechanisms of the external subset $\mathcal{C}$ and the
internal data, considering data prior to time $T^{*}$. In particular, if the
selected external subset $\mathcal{C}$ is generated following the same
mechanism as the internal data, the marginal likelihood value should be large.
On the other hand, if the external subset and internal data follow distinct
generating mechanisms, the marginal likelihood should be smaller. Motivated by
this property, the posterior distribution in (1) uses the product of the prior
and the marginal likelihood to quantify the likeliness that a specific
external subset should be selected. We give a theoretical justification in
Section 2.2.

Since the marginal likelihood plays an important role in our proposed method,
we provide more discussions here. In Bayesian studies, the marginal likelihood
is commonly used to quantify the goodness-of-fit for multiple competing models
(or priors). The models that assign higher probability masses around the true
parameter values (i.e., fit the data better) in general yield higher marginal
likelihood values (Robert et al., 2007). One commonly used Bayesian concept
derived based on this idea is the Bayes factor. It involves calculating the
ratio of the marginal likelihood values between two competitive models $M_{1}$
and $M_{2}$. In most cases, a Bayes factor would asymptotically converge to
either 0 or $+\infty$ (in probability) as the sample size goes to $\infty$,
representing $M_{1}$ or $M_{2}$ is preferred, respectively. In the context of
(1), the posterior density functions induced by different subset indexes
$\mathcal{C}$ can be interpreted as the candidate models. During the data
selection procedure, the external subset that possesses a higher marginal
likelihood can be viewed as having a superior ability to restore the internal
data preceding $T^{*}$. Note that it is possible that the selection process
ends up with selecting all patients in the external dataset. However, instead
of naively utilizing all external data, our approach uses a data-driven way to
determine which of the external data to incorporate, making it more flexible
and robust.

###  2.2 A theoretical perspective

We first consider a toy example to provide theoretical insights of BASE.
Consider the external data $X_{1},\ldots,X_{N_{0}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}p_{\theta_{0}}$ and the internal data
$Y_{1},\ldots,Y_{N_{1}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}p_{\theta_{1}}$, where $p_{\theta}$ is a known probability
density function with parameter $\theta$. Then we show in the following
theorem that by adopting the marginal likelihood in (1) as a criterion for
Bayesian estimation with a uniform prior on $\mathcal{C}$, our method is able
to correctly estimate $\mathcal{C}$ as the entire external subset if
$\theta_{0}$ and $\theta_{1}$ are sufficiently close and the empty set
otherwise.

######  Theorem 2.1.

Suppose $X_{1},\ldots,X_{N_{0}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}p_{\theta_{0}}$ and
$Y_{1},\ldots,Y_{N_{1}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}p_{\theta_{1}}$, and Assumptions (A1)-(A4) detailed in the
Supplementary Materials hold. If
$\|\theta_{0}-\theta_{1}\|_{1}\lesssim\epsilon_{N_{0}}$, the expected Bayes
factor can be controlled as follows,

$\begin{split}&\frac{\int\prod_{k=1}^{N_{1}}p_{\theta}(Y_{k})d\pi(p_{\theta})}{\int_{\|\theta-\theta_{0}\|_{1}\leq
K\epsilon_{N_{0}}}\prod_{k=1}^{N_{1}}p_{\theta}(Y_{k})d\pi(p_{\theta}\mid\left\\{X_{j}\right\\}_{j=1}^{N_{0}})}\stackrel{{\scriptstyle
p}}{{\to}}0,~{}~{}\text{as}~{}N_{0},N_{1}\to\infty,\end{split}$ (2)

On the other hand, if
$\|\theta_{0}-\theta_{1}\|_{1}\gtrsim\epsilon_{N_{0}}\times\psi_{N_{0}}$
instead, for $\psi_{N_{0}}$ diverging to $\infty$ at any rates, it yields that

$\begin{split}&\frac{\int_{\|\theta-\theta_{0}\|_{1}\leq
K\epsilon_{N_{0}}}\prod_{k=1}^{N_{1}}p_{\theta}(Y_{k})d\pi(p_{\theta}\mid\left\\{X_{j}\right\\}_{j=1}^{N_{0}})}{\int\prod_{k=1}^{N_{1}}p_{\theta}(Y_{k})d\pi(p_{\theta})}\stackrel{{\scriptstyle
p}}{{\to}}0,~{}~{}\text{as}~{}N_{0},N_{1}\to\infty,\end{split}$ (3)

where $a\lesssim b$ denotes $a\leq C\times b$ for a universal positive
constant $C$, $K>0$ is defined in the Supplementary File, $\pi(\cdot)$ is a
prior on $p_{\theta}$, $\epsilon_{n}=\left(\frac{M_{n}}{n}\right)^{1/2}$
refers to the posterior contraction rate and $M_{n}=o(\log n)$ is a sequence
of numbers going to $\infty$ as $n\rightarrow\infty$.

Theorem 2.1 can be viewed as a Bayesian testing (model selection) consistency
result. When the external data $\\{X_{i}\\}_{i=1}^{N_{0}}$ is similar to the
internal data $\\{Y_{j}\\}_{j=1}^{N_{1}}$, i.e., $\theta_{0}$ and $\theta_{1}$
are sufficiently close, the marginal likelihood of the internal data, obtained
by integrating over the posterior density given the external data, dominates
the likelihood integrated over the prior. In other words, the external data
provides additional useful information in this case compared to using internal
data only and the magnitude of marginal likelihood is able to reflect that.
Contrarily, when the external and internal data sets are significantly
distinct, the marginal likelihood merely given by the prior will be preferred,
i.e., the marginal likelihood only selects the internal data to use.

Assumptions (A1)-(A4) in Theorem 2.1 can be interpreted as follows: (A1) The
model parameter $\theta$ is finite-dimensional and the prior distribution of
$\theta$ has a large support; (A2) Both $\theta_{0}$ and $\theta_{1}$ can be
consistently estimated based on the external data $\\{X_{j}\\}_{j=1}^{N_{0}}$
and the internal data $\\{Y_{k}\\}_{k=1}^{N_{1}}$; (A3) The external data size
$N_{0}$ is smaller than the internal data size $N_{1}$ at a specific rate; and
(A4) The Kullback-Leibler neighbors include the $\ell_{1}$ neighbors, for
$\theta_{0}$ and $\theta_{1}$. Assumption (A1) is easily satisfied since
$p_{\theta}$ is a parametric model. Assumptions (A2) and (A4) are testable and
widely used in Bayesian asymptotics literature; and they can be easily
verified following the steps of Lemma 8.1 and Theorem 2.1 in Ghosal et al.
(2000) and Lemma B2 in Shen et al. (2013). Assumption (A3) requires the
internal data size to be larger than the external data size, which is
reasonable given that we would like to only borrow information from the
external set while treating the internal data as the main source for
inference.

Next we consider a more realistic scenario where the internal data
$Z_{1},\ldots,Z_{N_{1}}$ is still generated from $p_{\theta_{0}}$, but the
external data consist of $X_{1},\ldots,X_{N_{0}^{*}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}p_{\theta_{0}}$ and
$Y_{1},\ldots,Y_{N_{0}^{\prime}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}p_{\theta_{1}}$. In the next theorem, we show that our method
is capable of identifying the correct external subset
$\\{X_{i}\\}_{i=1}^{N_{0}^{*}}$, with a cost of at most including a small
amount of data from the incorrect external subset
$\\{Y_{j}\\}_{j=1}^{N_{0}^{\prime}}$. The key observation is that if a non-
ignorable proportion of $\\{Y_{j}\\}_{j=1}^{N_{0}^{\prime}}$ is included, the
induced density function, a mixture of $p_{\theta_{0}}$ and $p_{\theta_{1}}$,
diverges from $p_{\theta_{0}}$ at a recognizable distance, hence will be
detected by our method. Theorem 2.2 holds under Assumption (A1), (A3), (B1),
(C1) and (C2). Briefly speaking, (B1) and (C2) are the modifications of (A2)
and (A4) given pooled samples; (C1) assumes that $N_{0}^{*}$ tends to $\infty$
as $N_{0}\equiv N_{0}^{*}+N_{0}^{\prime}$ approaches $\infty$.

######  Theorem 2.2.

Suppose Assumptions (A1), (A3), (B1), (C1) and (C2) in the Supplementary File
hold. Let $m^{\prime}\leq N_{0}^{\prime}$ be the number of subjects selected
from the incorrect subset $\\{Y_{j}\\}_{j=1}^{N_{0}^{\prime}}$. If
$m^{\prime}$ satisfies
$\|p_{\theta_{0}}-p_{\theta_{1}}\|_{1}=o\left(\frac{m^{\prime}+N_{0}^{*}}{m^{\prime}}\times\left(\sqrt{\frac{M_{m^{\prime}+N_{0}^{*}}}{m^{\prime}+N_{0}^{*}}}-\sqrt{\frac{M_{N_{1}}}{N_{1}}}\right)\right)$,
we have

$\begin{split}&\frac{\int_{\|p_{\theta}-p_{u,v}^{*}\|_{1}\leq
K^{*}\epsilon_{u+v}}\prod_{k=1}^{N_{1}}p_{\theta}(Z_{k})d\pi(p_{\theta}\mid\left\\{X_{i}\right\\}_{i=1}^{u},\left\\{Y_{j}\right\\}_{j=1}^{v})}{\int_{\|p_{\theta}-p_{N_{0}^{*},m^{\prime}}^{*}\|_{1}\leq
K^{*}\epsilon_{N_{0}^{*}+m^{\prime}}}\prod_{k=1}^{N_{1}}p_{\theta}(Z_{k})d\pi(p_{\theta}\mid\left\\{X_{i}\right\\}_{i=1}^{N_{0}^{*}},\left\\{Y_{j}\right\\}_{j=1}^{m^{\prime}})}\stackrel{{\scriptstyle
p}}{{\to}}0,~{}~{}\text{as}~{}N_{0}^{*},N_{0},u,v\to\infty\end{split}$ (4)

if either of the following two conditions holds (i)
$\frac{u+v}{v}\times\left(\sqrt{\frac{M_{N_{1}}}{N_{1}}}+\sqrt{\frac{M_{u+v}}{u+v}}\right)=o\left(\|p_{\theta_{0}}-p_{\theta_{1}}\|_{1}\right)$
or (ii) $u+v=o(N_{0}^{*}+m^{\prime})$, where $M_{n}=o(\log n)$ is a sequence
of numbers going to $\infty$ as $n\rightarrow\infty$,
$\epsilon_{n}=\left(\frac{M_{n}}{n}\right)^{1/2}$, $u$ and $v$ denote the
number of subjects selected from the correct and incorrect external subsets,
respectively,
$p^{*}_{n,m}\equiv\frac{n}{n+m}p_{\theta_{0}}+\frac{m}{n+m}p_{\theta_{1}}$
refers to the weighted average of $p_{\theta_{0}}$ and $p_{\theta_{1}}$ given
$n\leq N_{0}^{*}$, $m\leq N_{0}^{\prime}$, and $n,m\in\mathbbm{Z}^{+}$.

Theorem 2.2 states that our model (induced by the posterior
$\pi(\cdot\mid\\{X_{i}\\}_{i=1}^{N_{0}^{*}},\\{Y_{j}\\}_{j=1}^{m^{\prime}})$)
that selects the entire correct external subset
$\\{X_{i}\\}_{i=1}^{N_{0}^{*}}$ plus $m^{\prime}$ subjects from the incorrect
subset $\\{Y_{j}\\}_{j=1}^{N_{0}^{\prime}}$ is preferred over alternative
models $\pi(\cdot\mid\\{X_{i}\\}_{i=1}^{u},\\{Y_{j}\\}_{j=1}^{v})$ in terms of
the marginal likelihood values. The proportion of incorrect selection
$\frac{m^{\prime}}{N_{0}^{*}}$ is decided by the discrepancy between
$p_{\theta_{0}}$ and $p_{\theta_{1}}$. For example, if the discrepancy is
small, e.g., 0, then $m^{\prime}$ can be as large as $N_{0}^{*}$ because there
is no incorrect subset. When the discrepancy is large, then $m^{\prime}$
becomes negligible compared to $N_{0}^{*}$. The number of subjects $u$ and $v$
are assumed to satisfy (i)
$p^{*}_{u,v}\equiv\frac{u}{u+v}p_{\theta_{0}}+\frac{v}{u+v}p_{\theta_{1}}$ is
recognizably distant from $p_{\theta_{0}}$, or (ii)
$\\{X_{i}\\}_{i=1}^{u},\\{Y_{j}\\}_{j=1}^{v}$ is not as informative as
$\\{X_{i}\\}_{i=1}^{N_{0}^{*}},\\{Y_{j}\\}_{j=1}^{m^{\prime}}$. Notably, when
$u+v$ is finite, the result still holds by slightly modifying the proof of
Theorem 2.1.

Theorem 2.1 and 2.2 justify the data selection procedure with the marginal
likelihood criteria. Specifically, the model can select the correct external
subset $\\{X_{i}\\}_{i=1}^{N_{0}^{*}}$ at an order of $O(N_{0}^{*})$, while
discarding the wrong external subset $\\{Y_{j}\\}_{j=1}^{N_{0}^{\prime}}$
except for a negligible proportion. Motivated by these theoretical results
under i.i.d. assumptions, our selection criterion (1) uses the first half
(preceding $T^{*}$) of the follow-up studies to explore the relevant external
subset. This reduces bias in estimating the early-stage parameters, compared
to directly combining the internal and external data. Moreover, introducing
the external trajectories is beneficial for revealing the plateau information.
The validity and performance of this procedure will be further studied in
Section 4 through numerous simulations.

##  3 Model specification for the hemophilia study

###  3.1 Concatenated Cubic Hermite spline

In this section, we introduce a parametric model to characterize the trend of
factor level. Motivated by a shared pattern among similar gene therapy
products, i.e., factor level that increases after the treatment may decrease
over time and eventually reach a plateau value (Nathwani et al., 2018; Shah et
al., 2023), we propose to model the outcome mean trend using a concatenation
of two Cubic Hermite splines. Specifically, for trajectory $i$ in study $s$,
we model the outcome mean trend using $\psi(t;\boldsymbol{\theta}_{si})$, for
$t\in[0,T]$, parametrized by
$\boldsymbol{\theta}_{si}\equiv(\mu_{s0i},m_{s0i},\mu_{s1i},m_{s1i},\mu_{s2i})$,
which is defined using two Cubic Hermite splines, over two intervals
$[0,\alpha]$ and $(\alpha,T]$ and concatenated at a turning point $\alpha$,
i.e.,

$\psi(t;\boldsymbol{\theta}_{si})=\begin{cases}h_{00}(\frac{t}{\alpha})\mu_{s0i}+h_{10}(\frac{t}{\alpha})\alpha
m_{s0i}+h_{01}(\frac{t}{\alpha})\mu_{s1i}+h_{11}(\frac{t}{\alpha})\alpha
m_{s1i},&t\in[0,\alpha],\\\
h_{00}(\frac{t-\alpha}{T-\alpha})\mu_{s1i}+h_{10}(\frac{t-\alpha}{T-\alpha})(T-\alpha)m_{s1i}+h_{01}(\frac{t-\alpha}{T-\alpha})\mu_{s2i},&t\in(\alpha,T],\end{cases}$
(5)

where $T$ is a pre-specified time point at which the outcome has reached a
plateau (in our application, $T$ represents five years). For trajectory $i$ in
study $s$, $\mu_{s0i}$ and $m_{s0i}$ represent the starting value and the
corresponding derivative at $t=0$, respectively, $\mu_{s1i}$ and $m_{s1i}$
respectively denote the factor level value and the corresponding derivative at
the turning point $\alpha$, and $\mu_{s2i}$ represents the final value at the
plateau point $t=T$, where the derivative is 0. In addition,
$h_{00},h_{10},h_{01},h_{11}$ are the Cubic Hermite basis functions, defined
as
$h_{00}(t)=2t^{3}-3t^{2}+1,h_{10}(t)=t^{3}-2t^{2}+t,h_{01}(t)=-2t^{3}+3t^{2},$
and $h_{11}(t)=t^{3}-t^{2}.$ In Figure 2, we illustrate a mean function
$\psi(t;\boldsymbol{\theta}_{s})$ with the parameters annotated.

Figure 2:  The blue dashed lines represent the derivatives at specific time
points, while the black solid line illustrates the mean value. The turning
point is indicated by the red dashed line.

Our choice of using a concatenated Cubic Hermite spline has three-fold
benefits. First, it facilitates incorporating the trend information, e.g., the
factor level that increases after the treatment may decrease over time, as
suggested by Nathwani et al. (2018). For instance, consider a population-level
trend for study $s$, namely,

$\boldsymbol{\theta}_{s}\equiv(\mu_{s0},m_{s0},\mu_{s1},m_{s1},\mu_{s2})=\sum_{i=1}^{n_{s}}\frac{1}{n_{s}}(\mu_{s0i},m_{s0i},\mu_{s1i},m_{s1i},\mu_{s2i}),$

where $n_{s}$ represents the number of trajectories in study $s$. If we impose
a constraint that $m_{s0}$ is positive and $m_{s1}$ is negative, the
concatenated Cubic Hermite spline is able to mimic the initially increasing
and subsequently decreasing trend. Second, the parameterization motivates a
strategy for selecting the relevant external data, i.e., the
$\boldsymbol{\theta}_{s}$ estimated based on the selected external data should
be close to the one given by the internal data. This provides empirical
evidence to explore external subsets perceived as similar to the internal data
preceding $T^{*}$. Here we assume the existence of external subsets with a
$\boldsymbol{\theta}_{s}$ similar to the one from the internal data. In
Section 4, we explore the situations when this assumption is violated. Third,
the concatenated Cubic Hermite spline can be expressed as a linear span of the
basis functions of a Cubic Hermite spline, piece-wisely within each of the two
time intervals segmented at $\alpha$. From a Bayesian computation perspective,
this property allows us to adopt a multivariate normal prior on
$\boldsymbol{\theta}_{s}$ for conjugacy and computational efficiency.

In the following section, we will outline the procedure for selecting the
relevant external data based on this parameterization. It is important to note
that adopting the concatenated Cubic Hermite spline to model the underlying
trend might introduce a potential risk of underfitting in practice due to its
parametric nature. However, the advantage lies in the efficient posterior
sampling. Moreover, considering the limited nature of the external data in the
collected hemophilia dataset, opting for a more complex model could result in
unstable inferences regarding our primary research interests, such as the
plateau value. Therefore, we choose concatenated Cubic Hermite splines to
trade mild underfitting for a substantial increase in modeling efficiency and
robustness.

###  3.2 Full model specification

We present our model before applying the data selection strategy in Section 2.
The proposed model can be outlined in the following hierarchical order.

$\begin{split}&Y_{1it}=\psi(t;\boldsymbol{\theta}_{1i})+\epsilon_{1it},~{}\text{for}~{}i=1,\ldots,N_{1},~{}Y_{0jt}=\psi(t;\boldsymbol{\theta}_{0j})+\epsilon_{0jt},~{}\text{for}~{}j\in\mathcal{C},\\\
&\psi(t;\boldsymbol{\theta}_{si})=\begin{cases}h_{00}(\frac{t}{\alpha})\mu_{s0i}+h_{10}(\frac{t}{\alpha})\alpha
m_{s0i}+h_{01}(\frac{t}{\alpha})\mu_{s1i}+h_{11}(\frac{t}{\alpha})\alpha
m_{s1i},&t\in[0,\alpha]~{}\text{and}~{}t\in\mathcal{T}_{s},\\\
h_{00}(\frac{t-\alpha}{T-\alpha})\mu_{s1i}+h_{10}(\frac{t-\alpha}{T-\alpha})(T-\alpha)m_{s1i}+h_{01}(\frac{t-\alpha}{T-\alpha})\mu_{s2i},&t\in(\alpha,T]~{}\text{and}~{}t\in\mathcal{T}_{s},\end{cases}\\\
&\boldsymbol{\theta}_{1i},\boldsymbol{\theta}_{0j}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}N(\beta^{*},\Sigma^{*}),~{}\text{for}~{}i=1,\ldots,N_{1},~{}j\in\mathcal{C},~{}\text{with}\\\
&\boldsymbol{\theta}_{si}\equiv(\mu_{s0i},m_{s0i},\mu_{s1i},m_{s1i},\mu_{s2i},)^{T},~{}\text{for}~{}s=0,1,\\\
&\beta^{*}\mid\Sigma^{*}\sim N^{(2^{+},4^{-})}(\beta_{0},\Sigma^{*}),\\\
&{\Sigma^{*}}^{-1}\sim W(\nu_{0},\Psi_{0}),\\\
&(\epsilon_{sit})_{t\in\mathcal{T}_{s}}\mid\sigma_{s}^{2},\rho\sim
N(0,\sigma_{s}^{2}\times\Sigma_{\rho}),~{}\sigma_{s}^{-2}\sim\text{Gamma}(a_{0},b_{0}),\\\
&\alpha\sim
tN(\text{mean}=2,\text{sd}=1;\text{lb}=0,\text{ub}=T),~{}\rho\sim\text{Unif}(0,1),\\\
&\mathcal{C}\sim\Pi\left(\mathcal{C}\mid\\{Y_{sit},i=1,\ldots,N_{s},t\in\mathcal{T}_{s},t\leq
T^{*},s=0,1\\}\right),\end{split}$ (6)

where $\theta_{1i}$ and $\theta_{0j}$ refer to the individual-wise parameter
of the concatenated Cubic Hermite spline for both external and internal data,
respectively. Here, $\alpha$ is common for $s=0$ and $s=1$ to increase the
external information borrowing ratio, $N^{(2^{+},4^{-})}$ refers to a
multivariate normal distribution that is truncated above 0 for the second
entry and below 0 for the fourth entry, as discussed in Section 3.1,
$\Sigma_{\rho}$ denotes a covariance matrix with the auto-regressive
structure, with the $(k,l)$ entry being $\rho^{|k-l|}$ for any
$k,l\in\mathcal{T}_{s}$, and
$tN(\text{mean}=2,\text{sd}=1;\text{lb}=0,\text{up}=T)$ refers to a normal
distribution $N(2,1)$ truncated between 0 and $T$. The hyperparameters such as
the mean and the standard deviation are derived from previous studies, i.e.,
Samelson-Jones et al. (2021), available in Section 5 of the Supplementary
File. This suggests that the turning point is likely to occur around the
second year since the treatment, with a one-year standard deviation as
uncertainty. The posterior distribution of $\mathcal{C}$ shown in the last
display of (6) is defined as follows,

$\begin{split}&\Pi\left(\mathcal{C}\mid\\{Y_{sit},i=1,\ldots,N_{s},t\in\mathcal{T}_{s},t\leq
T^{*},s=0,1\\}\right)\\\
&~{}~{}\propto\int\prod_{i=1}^{N_{1}}\prod_{t\in\mathcal{T}_{1},t\leq
T^{*}}L(Y_{1it}\mid\beta^{*},\Sigma^{*},\alpha,\rho,\sigma_{0}^{2},\sigma_{1}^{2})\times\\\
&~{}~{}~{}~{}~{}\pi(\beta^{*},\Sigma^{*},\alpha,\rho,\sigma_{0}^{2},\sigma_{1}^{2}\mid\\{Y_{0it},i\in\mathcal{C},t\in\mathcal{T}_{0},t\leq
T^{*}\\})d\beta^{*}d\Sigma^{*}d\alpha d\rho
d\sigma_{0}^{2}d\sigma_{1}^{2},\end{split}$ (7)

where
$L(\cdot\mid\beta^{*},\Sigma^{*},\alpha,\rho,\sigma_{0}^{2},\sigma_{1}^{2})$
denotes the likelihood function by marginalizing out the individual-wise trend
information, and
$\pi(\beta^{*},\Sigma^{*},\alpha,\rho,\sigma_{0}^{2},\sigma_{1}^{2}\mid\\{Y_{0it},i\in\mathcal{C},t\in\mathcal{T}_{0},t\leq
T^{*}\\})$ refers to the posterior density function given the first half of
the follow-up data.

Model (6) needs more interpretations for clarification. To account for the
high within-study variation observed in Figure 1(a), we consider a mixed-
effects model structure instead of solely relying on the Gaussian error term
$\epsilon_{sit}$ to absorb the variation. This is inspired by the observation
from Figure 1(a) that the internal (phases 2b and 3) trajectories appear to
have an increasing variation level progressively, which is hardly modeled
either by using a large $\sigma_{s}^{2}$ or imposing commonly-used temporal
structures. In addition, it is crucial to highlight that the variables
$Y_{1it}$ and $Y_{0it}$ in (6) have been z-transformed following the steps
specified in Section 1.3. The z-transformation, to some extent, forces the
starting value and the corresponding variance across the studies to be the
same, i.e., the first entry of $\beta^{*}$ to be 0 and the first diagonal term
of $\Sigma^{*}$ to be 1 for both the studies. This normalization potentially
increases the proportion of external data that can be borrowed by the internal
data based on Model (6).

For both simulation and real data analysis, we consider the following hyper-
parameter settings for the priors, $\beta_{0}=(0,0,0,0,0)^{T}$,
$\nu_{0}=0.01$, $\Psi_{0}=100\times\text{diag}(1,1,1,1,1)$, $a_{0}=0.01$,
$b_{0}=0.01$ and $\Pi(\mathcal{C})\propto 1$. Besides, for the real data
analysis, we choose $T^{*}=2$ and $T=6$, which refer to the latest endpoint of
the 2.5-year follow-up study and assume the factor level will approach the
plateau value at the sixth year after the treatment. The simulation settings
will be given in Section 4. We will conduct sensitivity analyses to study
whether it significantly changes the inference result by increasing $T$. For
the sampling procedure, we consider $1,000$ MCMC iterations, with the first
half being the burn-in samples. Empirically, the chosen priors are
sufficiently non-informative, and the iteration number is enough to approach
the stationary posterior distribution with a well-mixing. Additionally, the
most time-consuming task, when deployed on a server that operates at 3.80 GHz,
takes approximately 16 hours. This duration is considered bearable for a
Bayesian method. The efficiency of the task is attributed to the conjugacy
achieved by adopting the multivariate-normal inverse-Wishart distribution and
inverse-Gamma distribution, as defined in (6). In Section 4, we will further
investigate the performance of our proposed model.

##  4 Simulation

We conduct simulations to evaluate the performance of our method. Throughout,
we generate the internal data as follows: 40 trajectories are generated given
$\beta^{*}=(20,1,10,-0.1,8)$, $\Sigma^{*}=\text{diag}(25,0.1,9,0.1,4)$,
$\alpha=2.5$, $T=5$ and $\rho\in\\{0,0.5\\}$. The number of observations for
each trajectory $n_{i}$ follows $\text{Pois}(20)$ with the time grids
$t_{1},\ldots,t_{n_{i}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\text{Unif}(0,2.5)$.

We consider four data-generating processes (DGPs) for generating $K$
trajectories in the external study: (DGP 1) $\beta^{*}=(20,1,10,-0.1,8)$; (DGP
2) $\beta^{*}=(20,1,10,-0.1,12)$; (DGP 3) the first $K/2$ trajectories from
$(20,1,10,-0.1,8)$ and the second $K/2$ trajectories from
$\beta^{*}=(10,4,20,-1,8)$; (DGP 4) the first $K/2$ trajectories from
$(20,1,10,-0.1,8)$ and the second $K/2$ trajectories from $(10,4,20,-1,12)$.

In DGP 1 and DGP 2, we investigate the ability of our model to borrow long-
term information from the external data. We generate the external data using
the same process as the internal data before a specific time point, i.e.,
$(\mu_{s0},m_{s0},\mu_{s1},m_{s1})$ are the same for $s=0,1$. This allows us
to assess whether the internal data can benefit from the external subset
selection in achieving a precise plateau value estimate. Furthermore, we let
the plateau value be the same (DGP 1) or different (DGP 2) across the two
studies to examine the robustness of the plateau value estimated through the
selection procedure. The overall $\theta_{s}$ is different across the external
and internal data for DGP2 due to different $\mu_{s2}$ values. In summary, DGP
1 is to examine our method under the true model, and DGP 2 aims to check the
$\mu_{s2}$ estimation when $\mu_{s2}^{\prime}s$ are distinct.

In DGP 3 and DGP 4, we consider a mix of trajectories from the external data,
half following the same data-generating process as the internal data and the
other half do not, i.e., $m^{\prime}=N_{0}^{\prime}/2$ preceding $T^{*}$. The
main purpose is to examine whether our method has the ability to select the
correct subset, as well as the efficiency of information borrowing when the
external data are significantly distinct from the internal data. We let the
plateau value be the same (DGP 3) or different (DGP 4) across the two studies.
For DGP3, we let $\mu_{s2}$ stay the same regardless of the
$\theta_{s}^{\prime}$ values, and for DGP4 only allow subjects with the
correct $\theta_{s}^{\prime}$ to have the same $\mu_{s2}$ with the one in the
internal data design. The rest of the parameters are set the same for all
DGPs: $\Sigma^{*}=\text{diag}(25,0.1,9,0.1,4)$, $\alpha=2.5$,
$K\in\\{20,40\\}$, $T=5$ and $\rho\in\\{0,0.5\\}$. The number of observations
for each trajectory $n_{i}$ follows $\text{Pois}(20)$ with the time grid
$t_{1},\ldots,t_{n_{i}}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\text{Unif}(0,5)$.

To emulate the observation from the real data that the internal data are
censored after the 2.5-year time point, we terminate the collection of the
internal data by generating the time grid for each observation from
$\text{Unif}(0,2.5)$, while the total follow-up duration is $T=5$. The upper
bound of $2.5$ is intentionally set to be the same with $\alpha$ such that the
estimation of the fifth entry of $\beta^{*}$, the plateau value (beyond
$\alpha$), becomes challenging due to the lack of data without the
incorporation of external data. In Figure 3, we generate one simulated data
for each scenario. DGP 1 and 2 assume the same distribution between external
and internal data before 2.5 years, while DGP 3 and 4 do not. DGP 1 and 3 use
the same plateau values between external and interval data, while DGP 2 and 4
do not.

Figure 3:  Simulated sample data under each data generation process (DGP).

We assess the performance of our method using two metrics: the relative error
and the external proportion selected during the posterior sampling. The
relative error is utilized to measure the accuracy of both the starting value
and the plateau value estimation. We use $\ell_{\text{Start}}$ and
$\ell_{\text{Plateau}}$ to denote them, respectively, and define them as
follows,

$\begin{split}&\ell_{\text{Start}}=\text{median}\left\\{\sum_{k=1}^{2}\left|\frac{\hat{\beta}_{k,l}-\beta_{k}^{*}}{\beta_{k}^{*}}\right|\right\\}_{l=1}^{M},\\\
&\ell_{\text{Plateau}}=\text{median}\left\\{\left|\frac{\hat{\beta}_{5,l}-\beta_{5}^{*}}{\beta_{5}^{*}}\right|\right\\}_{l=1}^{M},\end{split}$
(8)

where $M$ is the number of posterior samples remaining after burn-in,
$\beta^{*}_{j}$ denotes the $j$-th entry of $\beta^{*}$ and
$\hat{\beta}_{j,l}$ is the corresponding posterior sample obtained at the
$l$-th iteration after the burn-in. To evaluate the efficiency of information
borrowing, we examine the external proportion, denoted by,

$\begin{split}&p_{Ext}=\text{median}\left\\{\frac{|\mathcal{C}_{l}|}{K}\right\\}_{l=1}^{M},\end{split}$
(9)

where $|\mathcal{C}_{l}|$ is the size of the external subset (number of
trajectories) selected at the $l$-th iteration after the burn-in, and $K$
refers to the size of the entire external data (total trajectories). Our
method is compared with the direct combination of the internal and external
data under $\ell_{\text{Start}}$ and $\ell_{\text{Plateau}}$, while $p_{Ext}$
is only considered for our model.

We also compare the metrics $\ell_{\text{Start}}$ and $\ell_{\text{Plateau}}$
for both methods using $p_{\text{ifSmall}}$ and $p_{\text{sPlateau}}$, which
counts the proportion of times when our method (BASE) achieves a smaller value
in $\ell_{\text{Start}}$ and $\ell_{\text{Plateau}}$, respectively than the
direct combination based on the posterior samples. Similarly, we let
$p_{\text{ifPrefer}}$ denote the frequency that the correct half of the
external data is favored for DGP 3 and 4, i.e., the number of selected
trajectories for the correct half is no less than that for the wrong half. The
results are presented in Table 1 and 2.

Table 1: Median (The interquartile range; Q3 - Q1) of the three statistics over 100 Monte Carlo replications for each of the four DGPs given different sample sizes $K$. | Data selection (BASE) | Direct combination
---|---|---
$K$ | $\rho$ | DGP | $\ell_{\text{Start}}$ | $\ell_{\text{Plateau}}$ | $p_{Ext}$ | $\ell_{\text{Start}}$ | $\ell_{\text{Plateau}}$
20 | 0.0 | 1 | 0.11 (0.04) | 0.18 (0.17) | 0.85 (0.11) | 0.11 (0.05) | 0.20 (0.21)
| 0.5 |  | 0.12 (0.04) | 0.18 (0.17) | 0.70 (0.15) | 0.12 (0.06) | 0.21 (0.22)
40 | 0.0 |  | 0.10 (0.04) | 0.21 (0.26) | 0.85 (0.10) | 0.10 (0.06) | 0.22 (0.30)
| 0.5 |  | 0.11 (0.05) | 0.20 (0.23) | 0.69 (0.18) | 0.11 (0.06) | 0.22 (0.31)
20 | 0.0 | 2 | 0.12 (0.04) | 0.37 (0.36) | 0.85 (0.11) | 0.11 (0.05) | 0.38 (0.39)
| 0.5 |  | 0.12 (0.04) | 0.35 (0.25) | 0.70 (0.20) | 0.12 (0.06) | 0.40 (0.35)
40 | 0.0 |  | 0.10 (0.05) | 0.44 (0.42) | 0.84 (0.08) | 0.10 (0.05) | 0.45 (0.42)
| 0.5 |  | 0.11 (0.04) | 0.43 (0.41) | 0.70 (0.15) | 0.11 (0.06) | 0.45 (0.44)
20 | 0.0 | 3 | 0.20 (0.11) | 0.21 (0.21) | 0.80 (0.15) | 0.24 (0.16) | 0.22 (0.28)
| 0.5 |  | 0.17 (0.09) | 0.23 (0.21) | 0.60 (0.15) | 0.25 (0.16) | 0.22 (0.28)
40 | 0.0 |  | 0.24 (0.16) | 0.24 (0.22) | 0.73 (0.30) | 0.33 (0.18) | 0.23 (0.26)
| 0.5 |  | 0.23 (0.11) | 0.23 (0.19) | 0.60 (0.15) | 0.33 (0.17) | 0.25 (0.26)
20 | 0.0 | 4 | 0.20 (0.11) | 0.22 (0.23) | 0.80 (0.15) | 0.30 (0.18) | 0.43 (0.25)
| 0.5 |  | 0.17 (0.08) | 0.24 (0.20) | 0.60 (0.11) | 0.30 (0.17) | 0.43 (0.25)
40 | 0.0 |  | 0.25 (0.15) | 0.26 (0.26) | 0.73 (0.25) | 0.40 (0.15) | 0.54 (0.20)
| 0.5 |  | 0.22 (0.09) | 0.26 (0.24) | 0.60 (0.13) | 0.40 (0.16) | 0.53 (0.21)
Table 2: The three frequency values that compare the data selection procedure to the direct combination. $K$ | $\rho$ | DGP | $p_{\text{ifSmall}}$ | $p_{\text{sPlateau}}$ | DGP | $p_{\text{ifSmall}}$ | $p_{\text{ifPrefer}}$ | $p_{\text{sPlateau}}$
---|---|---|---|---|---|---|---|---
20 | 0.0 | 1 | 0.44 | 0.63 | 3 | 0.89 | 0.74 | 0.44
| 0.5 |  | 0.39 | 0.61 |  | 0.90 | 0.97 | 0.5
40 | 0.0 |  | 0.44 | 0.73 |  | 0.97 | 0.74 | 0.62
| 0.5 |  | 0.39 | 0.74 |  | 1.00 | 0.90 | 0.67
20 | 0.0 | 2 | 0.49 | 0.62 | 4 | 0.96 | 0.71 | 0.91
| 0.5 |  | 0.42 | 0.65 |  | 0.93 | 0.97 | 0.82
40 | 0.0 |  | 0.45 | 0.66 |  | 0.99 | 0.74 | 0.97
| 0.5 |  | 0.39 | 0.64 |  | 1.00 | 0.86 | 0.97

We first discuss Table 1. When data preceding $T^{*}$ are generated under the
true model (DGP 1 and 2), both methods achieve comparable accuracy in
estimating the starting parameters, indicated by the close
$\ell_{\text{Start}}$ values, which is of our expectation. On the other hand,
when data are generated incorrectly (DGP 3 and 4), the data selection
procedure demonstrates a remarkable reduction in $\ell_{\text{Start}}$
compared to the direct combination method. This observation is also supported
by the $p_{\text{ifSmall}}$ columns in Table 2, i.e., $p_{\text{ifSmall}}\geq
89\%$ for DGP 3 and 4. For the plateau value estimation, our method also
achieves a more accurate estimation than direct combination, as evident from
the observation that the $p_{\text{sPlateau}}$ values are consistently above
60% in Table 2. This advantage becomes particularly pronounced under DGP 4,
where a subset of the external data is generated differently from the internal
data, our method achieves significantly lower values in $\ell_{\text{Start}}$
and $\ell_{\text{Plateau}}$, along with high values in $p_{\text{ifPrefer}}$
and $p_{\text{sPlateau}}$. This underscores that our method effectively
selects the correct subset from the external data with high probabilities.

We have conducted additional simulation studies to assess the sensitivity of
our method by varying hyperparameter values and plateau points. Moreover, we
explored the impact of generating random errors from heavy-tailed
t-distributions instead of the normal distribution in our model. The results
are summarized in Section 3 and 4 of the Supplementary File, and they are
consistent with the findings in Table 1 and 2. All these findings confirm the
robustness of our method and its effectiveness in terms of selecting the
correct subset from external data and reducing estimation bias in parameter
estimation.

##  5 Long-term outcome after a hemophilia gene therapy

###  5.1 Data description and clinical question of interest

The data were collected from a three-phase study designed to assess the
effectiveness of gene therapy as a long-term therapeutic solution for
hemophilia B. The study focuses on monitoring the factor ix@ (Fix@) activity
levels of the subjects after receiving treatment with a therapy involving an
adeno associated virus 5 (AAV5) vector. In this work, we treat the 10
trajectories from the phase 1 study, with a 5-year follow-up, as the external
data. The 3 trajectories from the phase 2b study and the 54 trajectories from
the phase 3 study are considered as the internal data since the follow-up is
less than 30 months. The clinical question of interest is to make annual
predictions on the Fix@ activity level and provide statistical evidence to
validate whether the Fix@ activity level increases at the end of the fifth
year compared with the initial level, based on the prediction.

We provide a brief overview of the design and objectives of the three study
phases, with more detailed information available in the cited references for
interested readers. Phase 1 study (NCT02396342;Miesbach et al. (2021)) focused
on assessing the safety and efficacy of AMT-060 gene therapy and involved 10
male participants with severe or moderate-severe hemophilia B. Phase 2b study
(NCT03489291; Gomez et al. (2021)) was an open-label, single-dose, single-arm,
multi-center trial, and it involved 3 adult participants with severe
hemophilia B. The main goal was to assess the sustained efficacy of HEMGENIX
(also called by AMT-061), which is the first FDA approved gene therapy for
hemophilia B. Phase 3 study (NCT03569891; Pipe et al. (2023)) enrolled 54 men
with moderate-to-severe hemophilia B to evaluate the change in the annualized
bleeding rate post-treatment compared to the rate during the lead-in period.
The participants received a single infusion of HEMGENIX. Detailed demographic
and clinical characteristics of the participants in Phase 3 study can be found
in Pipe et al. (2023) and Table 1 therein.

###  5.2 Primary analysis

We implement our data selection method to analyze the Hemophilia data. For the
posterior sampling, we use the same hyperparameter settings as those
introduced in Section 3.2 and consider 1,000 MCMC iterations, with the first
500 iterations discarded as the burn-in samples. To obtain a representative
result, we conducted our selection procedure 100 times and studied the
variation of the posterior median for both the starting value and the plateau
value. The result indicates that the median values for the starting value and
the plateau value are 26.86 (with a standard deviation of 0.09) and 32.08
(with a standard deviation of 0.21), respectively, over the 100 replications.
These findings further indicate that it is trustworthy to study a single chain
because the across-chain variation is tolerable. Based on one of the 100
chains, we visualize the estimated trend in the top panel of Figure 4 using
the posterior median of $\beta^{*}$, with a $95\%$ credible region banded by
the $2.5\%$ and $97.5\%$ quantiles of the posterior samples. In particular, we
mark out the $95\%$ credible interval for each of the five years aligned with
the estimated trend.

Figure 4:  (a) The spaghetti plot of the internal trajectories and the trend
estimated based on the BASE. The green banded region refers to the $95\%$
credible interval; (b) The external trajectories presented with their opacity
proportional to the corresponding posterior frequency of being selected. The
dot-dashed line refers to the estimated trend.

According to the posterior samples of the chain presented, the median values
of the predicted starting and plateau values are 26.93 (with a 95$\%$ credible
interval of $[23.21,30.11]$) and 31.94 (with a 95$\%$ credible interval of
$[22.17,37.17]$), respectively. Moreover, the $95\%$ credible interval at Year
3 is $[27.65,38.77]$ and it contains the median Fix@ activity level (36.0)
reported in Phase 3 study (available at
https://ash.confex.com/ash/2023/webprogram/Paper187624.html). This validates
the accuracy of our proposed method to a good extent. To answer the scientific
question of whether there is statistically significant evidence to assert that
the predicted plateau value is greater than the starting value, we further
examine the posterior probability that the predicted plateau value exceeds the
predicted starting value. Based on the posterior samples, the estimated
probability is $93\%$, which strongly supports our hypothesis. Another finding
is that the bandwidth of the credible interval increases annually with the
trend over the five years, as shown in Figure 4(a). This finding is partly
attributed to the reduction in the effective sample size informing the trend
as it progresses, primarily due to the internal data being censored after the
2.5-year follow-up endpoint.

We also investigate to which extent each external trajectory is preferred by
the internal dataset, quantified by the posterior probability of selection.
These results are presented in Figure 4(b) using opacity, where lighter colors
indicate smaller probabilities, representing the degree to which each external
trajectory is preferred. We highlight three trajectories that are least
preferred as they have the lowest posterior probabilities of being selected
(e.g., 0.69, 0.71, and 0.73). In addition, the posterior median of the
selected external proportion is 0.80, with a standard deviation of 0.25. This
result suggests that the majority of the external trajectories can be utilized
based on the internal information, implying a significant similarity between
the internal and external data. From a scientific perspective, the high
proportion of external data utilization may indicate that the mechanism of
action of hemophilia and durability of gene therapy remains consistent between
AMT-060 (external data) and etranacogene dezaparvovec, which differ by 1
single amino acid resulting in naturally occurring highly active FIX Padua
variant (FIX-R338L) . Furthermore, it supports the hypothesis that this
evolution only influence the baseline and plateau factor levels after the
subjects are treated. Moreover, by examining the study subjects that are
frequently selected, one could identify potential external sub-populations
exhibiting the same trend as the internal population before the follow-up
endpoint. Such insights could be valuable for future research endeavors.

###  5.3 Secondary analysis

We conduct a two-fold secondary analysis. First, we compare the results given
by our proposed method (Selection) with those from the direct combination
(Combination) and no-information-borrowing (No-borrow) from the external
dataset. Second, we evaluate the robustness of our model by considering
different prior hyper-parameter settings.

We present the estimated trends along with their credible intervals obtained
by the three methods under the same MCMC settings in Figure 5(a) and (b).
Additionally, we display the estimated annual factor levels and the turning
point in Table 3. All three methods yield similar results in the early stage,
specifically before the second year, which is due to the availability and
dominance of the internal trajectories before 2.5 years. The curve obtained
without using external information (No-borrow) exhibits high variance starting
from the middle period of the 5-year duration due to the absence of long-term
data. This is also reflected by the large standard deviations in the estimated
factor levels at the fourth and fifth years in Table 3. Both our method and
the direct combination method yield similar estimated trends and plateau
values. This similarity arises because our method selects a high proportion
(around $80\%$) of the external data for inference. However, our method
achieves this estimate in a more principled way.

Figure 5:  (a) The estimated trends (median) together with their $95\%$ credible intervals obtained by the three methods; (b) The zoom-in plot of (a); (c) The sensitivity analysis results under different hyper-parameter settings; (d) The sensitivity analysis results with different plateau time points. Table 3: The posterior median value (standard deviation) of the annual factor level and the turning point $\alpha$. Method | Start | 1st year | 2nd year
---|---|---|---
Selection | 26.93 (1.77) | 39.59 (2.96) | 35.20 (2.64)
Combination | 26.89 (1.68) | 40.31 (2.96) | 35.87 (2.41)
No-borrow | 26.45 (1.67) | 40.80 (3.10) | 35.14 (2.34)
3rd year | 4th year | 5th year | $\alpha$
33.75 (2.87) | 32.60 (3.15) | 31.94 (3.64) | 1.93 (0.09)
33.68 (3.21) | 32.26 (3.19) | 31.40 (3.09) | 2.01 (0.06)
29.58 (3.79) | 20.90 (9.92) | 12.47 (17.80) | 2.07 (0.07)

To investigate the robustness of our results, we perform a sensitivity
analysis by adopting alternative prior hyperparameter values
$\nu_{0},\Psi_{0},a_{0},b_{0}$ defined in (6) and plateau time points $T$.
Specifically, we consider a more informative prior (Informative) and a vaguer
prior (Vague) in comparison to the current setting, and larger plateau time
points ($T=7$ and $8$). More details about these setting are given in Table 1
of the Supplementary File. As shown in Figure 5(c) and (d), the credible bands
exhibit a substantial amount of agreement, which confirms a satisfactory level
of robustness, i.e., the results obtained are not sensitive to the choice of
the hyperparameter values and plateau time points. We have also conducted a
residual diagnosis analysis and explored the use of heavy-tailed residuals in
the model, with results presented in Sections 3 and 4 of the Supplementary
File.

##  6 Discussion

In this paper, we propose a data integration method using Bayesian procedures.
Our work stands out as the first work that selects subsets from external
datasets based on theoretically supported motivations. Compared to the direct
combination method, our approach exhibits extra flexibility and
interpretability, providing inspiring scientific evidence for future clinical
studies. Moreover, our work not only provides a definite answer on the long
term efficacy of HEMGENIX, but also showcases a conceptually novel procedure
that manages to refine the long-term efficacy inference. For instance, the
factor level on the 15-th year post treatment can be predicted if the trials
continue, rendering more available durability data.

Despite these achievements, the current model can be generalized in future
studies. Firstly, the current MCMC sampling scheme relies on the proposal of
the external subset $\mathcal{C}$. Although the current proposal has a
satisfactory empirical performance in simulations and real data analysis,
e.g., the standard deviation of the external proportion is moderate compared
to the point estimate, it remains unexplored whether this proposal is
effective when the external data size is large, for instance, over 100
observations. Intuitively, a more effective proposal should be devised in such
scenarios, given that the number of possible external subsets increases at an
exponential order with the external data size. However, since the size of the
hemophilia external data is manageable (10 observations), this question is
beyond the scope of this work. Additionally, the estimation of the marginal
likelihood can be further investigated. The current strategy (Pajor, 2017) is
proven to be effective under the parametric setting, whilst there is limited
evidence that this method provides efficient Monte Carlo estimation to the
marginal likelihood under the nonparametric setting. Moreover, when the sample
size is sufficiently large, one can consider data-driven strategies that
select models from a family of parametric models. Criteria such as the Akaike
Information Criterion (AIC; Akaike (1974)), Bayesian Information Criterion
(BIC;Schwarz (1978)), and marginal likelihood can be employed to explore if
any alternatives are superior in capturing the mean trend compared to the
Hermite spline specification. Another intriguing future topic is how to embed
the spirit of data selection into non-likelihood problems, e.g., treatment
effects that are given by estimating equations and non-probability samples
(Yang et al., 2020; Chen et al., 2020). In cases where likelihood functions
are not well-defined, it would be impossible to use the marginal likelihood as
a guideline, rendering our current pipeline unsuitable. Furthermore, it is
also challenging to determine whether the internal estimation benefits from
the external data selection.

## References

  * Akaike (1974) Akaike, H. (1974).  A new look at the statistical model identification.  IEEE transactions on automatic control 19(6), 716–723.
  * Athey et al. (2019) Athey, S., R. Chetty, G. W. Imbens, and H. Kang (2019).  The surrogate index: Combining short-term proxies to estimate long-term treatment effects more rapidly and precisely.  Technical report, National Bureau of Economic Research.
  * Cai et al. (2016) Cai, T., T. T. Cai, and A. Zhang (2016).  Structured matrix completion with applications to genomic data integration.  Journal of the American Statistical Association 111(514), 621–633.
  * Chen et al. (2000) Chen, M.-H., J. G. Ibrahim, and Q.-M. Shao (2000).  Power prior distributions for generalized linear models.  Journal of Statistical Planning and Inference 84(1-2), 121–137.
  * Chen et al. (1999) Chen, M.-H., J. G. Ibrahim, and C. Yiannoutsos (1999).  Prior elicitation, variable selection and bayesian computation for logistic regression models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 61(1), 223–242.
  * Chen et al. (2021) Chen, S., B. Zhang, and T. Ye (2021).  Minimax rates and adaptivity in combining experimental and observational data.  arXiv preprint arXiv:2109.10522.
  * Chen et al. (2020) Chen, Y., P. Li, and C. Wu (2020).  Doubly robust inference with nonprobability survey samples.  Journal of the American Statistical Association 115(532), 2011–2021.
  * Dai and Li (2022) Dai, X. and L. Li (2022).  Orthogonalized kernel debiased machine learning for multimodal data analysis.  Journal of the American Statistical Association, 1–15.
  * Farmer et al. (2023) Farmer, C., E. Nikram, L. A. Trigg, S. Robinson, J. Coppell, M. Muthukumar, G. Melendez-Torres, and E. C. Wilson (2023).  Etranacogene dezaparvovec for treating moderately severe or severe haemophilia b.  National Institute for Health and Care Excellence.
  * George et al. (2017) George, L. A., S. K. Sullivan, A. Giermasz, J. E. Rasko, B. J. Samelson-Jones, J. Ducore, A. Cuker, L. M. Sullivan, S. Majumdar, J. Teitel, et al. (2017).  Hemophilia b gene therapy with a high-specific-activity factor ix variant.  New England Journal of Medicine 377(23), 2215–2227.
  * Ghosal et al. (2000) Ghosal, S., J. K. Ghosh, and A. W. Van Der Vaart (2000).  Convergence rates of posterior distributions.  Annals of Statistics, 500–531.
  * Gomez et al. (2021) Gomez, E., A. Giermasz, G. Castaman, N. Key, S. Lattimore, F. Leebeek, W. Miesbach, M. Recht, A. Drygalski, E. Sawyer, and S. Pipe (2021).  Etranacogene dezaparvovec (aav5-padua hfix variant, amt-061), an enhanced vector for gene transfer in adults with severe or moderate-severe hemophilia B: 2.5 year data from a phase 2b trial.  ISTH Congress.
  * Hassler et al. (2023) Hassler, G. W., A. F. Magee, Z. Zhang, G. Baele, P. Lemey, X. Ji, M. Fourment, and M. A. Suchard (2023).  Data integration in bayesian phylogenetics.  Annual Review of Statistics and Its Application 10.
  * Ibrahim et al. (2015) Ibrahim, J. G., M.-H. Chen, Y. Gwon, and F. Chen (2015).  The power prior: theory and applications.  Statistics in medicine 34(28), 3724–3749.
  * Imbens et al. (2022) Imbens, G., N. Kallus, X. Mao, and Y. Wang (2022).  Long-term causal inference under persistent confounding via data combination.  arXiv preprint arXiv:2202.07234.
  * Kee and Maio (2019) Kee, A. and V. Maio (2019).  Value-based contracting: challenges and opportunities.  American Journal of Medical Quality 34(6), 615–617.
  * Leebeek and Miesbach (2021) Leebeek, F. W. and W. Miesbach (2021).  Gene therapy for hemophilia: a review on clinical benefit, limitations, and remaining issues.  Blood 138(11), 923–931.
  * Li and Li (2022) Li, Q. and L. Li (2022).  Integrative factor regression and its inference for multimodal data analysis.  Journal of the American Statistical Association 117(540), 2207–2221.
  * Li and Luedtke (2021) Li, S. and A. Luedtke (2021).  Efficient estimation under data fusion.  arXiv preprint arXiv:2111.14945.
  * Liu et al. (2022) Liu, Y., X. Sun, W. Zhong, and B. Li (2022).  B-scaling: A novel nonparametric data fusion method.  The Annals of Applied Statistics 16(3), 1292–1312.
  * Miesbach et al. (2021) Miesbach, W., K. Meijer, M. Coppens, P. Kampmann, R. Klamroth, R. Schutgens, G. Castaman, S. E.K, and L. F.WG (2021).  Five year data confirms stable fix expression and sustained reductions in bleeding and factor ix use following amt-060 gene therapy in adults with severe or moderate-severe hemophilia B.  ISTH Congress.
  * Nathwani et al. (2018) Nathwani, A. C., U. Reiss, E. Tuddenham, P. Chowdary, J. McIntosh, A. Riddell, J. Pie, J. N. Mahlangu, M. Recht, Y.-M. Shen, K. G. Halka, M. M. Meagher, A. W. Nienhuis, A. M. Davidoff, S. Mangles, C. L. Morton, Z. Junfang, and V. C. Radulescu (2018, 11).  Adeno-associated mediated gene transfer for hemophilia B:8 year follow up and impact of removing ”empty viral particles” on safety and efficacy of gene transfer.  Blood 132(Supplement 1), 491–491.
  * Nathwani et al. (2014) Nathwani, A. C., U. M. Reiss, E. G. Tuddenham, C. Rosales, P. Chowdary, J. McIntosh, M. Della Peruta, E. Lheriteau, N. Patel, D. Raj, et al. (2014).  Long-term safety and efficacy of factor ix gene therapy in hemophilia B.  New England Journal of Medicine 371(21), 1994–2004.
  * Nathwani et al. (2011) Nathwani, A. C., E. G. Tuddenham, S. Rangarajan, C. Rosales, J. McIntosh, D. C. Linch, P. Chowdary, A. Riddell, A. J. Pie, C. Harrington, et al. (2011).  Adenovirus-associated virus vector–mediated gene transfer in hemophilia B.  New England Journal of Medicine 365(25), 2357–2365.
  * Neuenschwander et al. (2009) Neuenschwander, B., M. Branson, and D. J. Spiegelhalter (2009).  A note on the power prior.  Statistics in medicine 28(28), 3562–3566.
  * Oberst et al. (2022) Oberst, M., A. D’Amour, M. Chen, Y. Wang, D. Sontag, and S. Yadlowsky (2022).  Bias-robust integration of observational and experimental estimators.  arXiv preprint arXiv:2205.10467.
  * Pajor (2017) Pajor, A. (2017).  Estimating the marginal likelihood using the arithmetic mean identity.
  * Park et al. (2021) Park, S., H. Xu, and H. Zhao (2021).  Integrating multidimensional data for clustering analysis with applications to cancer patient data.  Journal of the American Statistical Association 116(533), 14–26.
  * Pipe et al. (2023) Pipe, S. W., F. W. Leebeek, M. Recht, N. S. Key, G. Castaman, W. Miesbach, S. Lattimore, K. Peerlinck, P. Van der Valk, M. Coppens, et al. (2023).  Gene therapy with etranacogene dezaparvovec for hemophilia B.  New England Journal of Medicine 388(8), 706–718.
  * Ritchie et al. (2015) Ritchie, M. D., E. R. Holzinger, R. Li, S. A. Pendergrass, and D. Kim (2015).  Methods of integrating data to uncover genotype–phenotype interactions.  Nature Reviews Genetics 16(2), 85–97.
  * Robert et al. (2007) Robert, C. P. et al. (2007).  The Bayesian choice: from decision-theoretic foundations to computational implementation, Volume 2.  Springer.
  * Samelson-Jones et al. (2021) Samelson-Jones, B. J., S. K. Sullivan, J. E. Rasko, A. Giermasz, L. A. George, J. M. Ducore, J. M. Teitel, C. E. McGuinn, A. O’Brien, I. Winburn, et al. (2021).  Follow-up of more than 5 years in a cohort of patients with hemophilia b treated with fidanacogene elaparvovec adeno-associated virus gene therapy.  Blood 138, 3975.
  * Schwarz (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  The annals of statistics, 461–464.
  * Shah et al. (2023) Shah, J., H. Kim, K. Sivamurthy, P. E. Monahan, and M. Fries (2023).  Comprehensive analysis and prediction of long-term durability of factor ix activity following etranacogene dezaparvovec gene therapy in the treatment of hemophilia B.  Current Medical Research and Opinion 39(2), 227–237.
  * Shen et al. (2013) Shen, W., S. T. Tokdar, and S. Ghosal (2013).  Adaptive bayesian multivariate density estimation with dirichlet mixtures.  Biometrika 100(3), 623–640.
  * Tice et al. (2022) Tice, J., S. Walton, B. Herce-Hagiwara, S. Fahim, A. Moradi, J. Sarker, J. Chu, F. Agboola, S. Pearson, and D. Rind (2022, 12).  Gene therapy for hemophilia B and an update on gene therapy for hemophilia A: Effectiveness and value.  Institute for Clinical and Economic Review.
  * Valentinuzzi (2004) Valentinuzzi, M. E. (2004).  Friedman lm, furberg cd, demets dl: Fundamentals of clinical trials 3rd edition: Springer-verlag, new york; 1998. 361 pp, isbn 0-387-98586-7.
  * Xue and Qu (2021) Xue, F. and A. Qu (2021).  Integrating multisource block-wise missing data in model selection.  Journal of the American Statistical Association 116(536), 1914–1927.
  * Yang et al. (2020) Yang, S., C. Gao, D. Zeng, and X. Wang (2020).  Elastic integrative analysis of randomized trial and real-world data for treatment heterogeneity estimation.  arXiv preprint arXiv:2005.10579.
  * Yang et al. (2020) Yang, S., J. K. Kim, and R. Song (2020).  Doubly robust inference when combining probability and non-probability samples with high dimensional data.  Journal of the Royal Statistical Society. Series B, Statistical Methodology 82(2), 445.