# An adaptive enrichment design using Bayesian model averaging for selection
and threshold-identification of tailoring variables

Lara Maleyeff
Department of Epidemiology, Biostatistics, and Occupational Health, McGill
University, Montréal, QC, CA
<EMAIL_ADDRESS>
&Shirin Golchi
Department of Epidemiology, Biostatistics, and Occupational Health, McGill
University, Montréal, QC, CA
&Erica E.M. Moodie
Department of Epidemiology, Biostatistics, and Occupational Health, McGill
University, Montréal, QC, CA
&Marie Hudson
Department of Medicine, McGill University, Montréal, QC, CA
Jewish General Hospital and Lady Davis Institute for Medical Research,
Montréal, QC, CA

###### Abstract

Precision medicine stands as a transformative approach in healthcare, offering
tailored treatments that can enhance patient outcomes and reduce healthcare
costs. As understanding of complex disease improves, clinical trials are being
designed to detect subgroups of patients with enhanced treatment effects.
Biomarker-driven adaptive enrichment designs, which enroll a general
population initially and later restrict accrual to treatment-sensitive
patients, are gaining popularity. Current practice often assumes either pre-
trial knowledge of biomarkers defining treatment-sensitive subpopulations or a
simple, linear relationship between continuous markers and treatment
effectiveness. Motivated by a trial studying rheumatoid arthritis treatment,
we propose a Bayesian adaptive enrichment design which identifies important
tailoring variables out of a larger set of candidate biomarkers. Our proposed
design is equipped with a flexible modelling framework where the effects of
continuous biomarkers are introduced using free knot B-splines. The parameters
of interest are then estimated by marginalizing over the space of all possible
variable combinations using Bayesian model averaging. At interim analyses, we
assess whether a biomarker-defined subgroup has enhanced or reduced treatment
effects, allowing for early termination due to efficacy or futility and
restricting future enrollment to treatment-sensitive patients. We consider
pre-categorized and continuous biomarkers, the latter of which may have
complex, nonlinear relationships to the outcome and treatment effect. Using
simulations, we derive the operating characteristics of our design and compare
its performance to two existing approaches.

##  1 Introduction

Precision medicine, which allows for decisions about medical interventions to
depend on patient characteristics, may be more efficient, less costly, and
result in better patient outcomes, on average, than a traditional one-size-
fits-all approach. This is particularly true in domains where response to
treatment is highly heterogeneous, such as in rheumatoid arthritis (RA) and
oncology. For example, of the several treatment options available for RA, 30%
of patients fail any particular drug (Smolen et al., 2018). The current trial-
and-error approach in RA, in which patients attempt multiple drugs over the
span of months or years until their disease is controlled, results in
unnecessary comorbidity and potential patient harm. In the presence of such
heterogeneous treatment effects, traditional clinical trials, which detect a
population-averaged effect, may be inefficient (Adams and Brantner, 2006). As
molecular phenotyping becomes increasingly available and inexpensive to
obtain, clinical trials to assess treatment efficacy in biomarker-defined
subgroups are becoming more feasible. Such validated, quality information is
crucial to determine appropriate personalized treatments. In addition,
incorporating statistically sound subgroup-identification methods can lessen
the risk of conducting null trials and reduce exposure of patients to
treatments which do not benefit them.

In the clinical trial setting, several biomarker-driven designs exist to
detect heterogeneous effects. Adaptive enrichment designs initially enroll and
randomize patients regardless of their biomarker status, then use interim data
to cease enrollment of certain biomarker-defined subgroups for efficacy or
futility. Wang et al. (2007) developed a two-stage adaptive design to assess
treatment effects across levels of a dichotomous biomarker; multiple similar
extensions have been proposed (Karuri and Simon, 2012; Magnusson and Turnbull,
2013; Rosenblum et al., 2016; Brannath et al., 2009; Jenkins et al., 2011;
Mehta et al., 2014). Basket trials have become more prolific in recent years
due their ability to cost-effectively evaluate the activity of a single
treatment on groups of patients with varying biomarker status, or “baskets"
(Asano and Hirakawa, 2020; Chen and Lee, 2020; Chu and Yuan, 2018; Jin et al.,
2020; Park et al., 2020). In such trials, the goal is to determine the subset
of baskets in which treatment is most effective. Development of these designs
is a rich area of active research; a comprehensive review can be found in
(Simon, 2017). Psioda et al. (2021) developed a Bayesian basket trial design
methodology which uses Bayesian model averaging (BMA). BMA, a coherent
mechanism for accounting for inherent uncertainly in model selection, yields
posterior distributions which are marginalized over a space of submodels. This
approach naturally allows for information borrowing by weighing multiple
models based on their posterior probabilities, leading to more robust and
informed decision-making.

Most biomarker-driven designs, including those described above, focus on pre-
determined patient subgroups. Such subgroups are defined by biomarkers which
are categorized, either as levels of naturally categorized variables or by
using pre-specified thresholds. However, there exist some trial designs which
incorporate biomarker identification. Freidlin and Simon (2005) and Freidlin
et al. (2010) propose a two-stage approach, wherein the first stage machine
learning techniques are used to identify a treatment-sensitive population; in
the second stage, the effectiveness of treatment is assessed for the sensitive
group. There also exist designs which incorporate variable selection at the
interim analysis stage. Gu et al. (2016) propose a multi-stage design, called
BATTLE-2, in which Bayesian lasso is used at interim analyses to refine the
set of biomarkers considered in the remainder of the trial. Park et al. (2022)
propose a Bayesian group sequential enrichment design which incorporates
variable selection to reduce the space of candidate biomarkers and uses an
early response proxy outcome that is used in interim analyses. At each interim
decision, spike-and-slab priors are used to identify a set of tailoring
biomarkers that characterizes treatment-sensitive patients and decisions are
made to stop the trial for efficacy or futility. Subsequent enrollment is
limited to the most recently identified treatment-sensitive subgroup. Their
methods can accommodate both categorical and continuous biomarkers, assuming
continuous biomarkers have a monotone, linear relationship with treatment
effect. This assumption is commonly made in much of the threshold
identification literature (Jiang et al., 2007; Renfro et al., 2014; Ohwada and
Morita, 2016).

Recently, Liu et al. (2022) developed a Bayesian adaptive enrichment design
which allows for a nonlinear and nonmonotone relationship between one or more
continuous biomarkers and treatment. Rather than utilizing algorithmic
dichotomization to determine biomarker thresholds, they flexibly estimate
continuous biomarker relationships directly and quantify the posterior
uncertainty around these estimates. They specify an interim analysis procedure
in which subgroups are identified and assessed for efficacy or futility. If no
subgroup meets the pre-specified criteria, an overall analysis is performed.
In simulation studies, their design had similar or improved operating
characteristics compared with more traditional grid-search based approaches.
While simulation studies showed positive performance, this method does not
incorporate variable selection; that is, the set of tailoring biomarkers is
fixed and assumed to be known a priori.

We propose a Bayesian adaptive enrichment clinical trial which identifies
important tailoring variables out of a larger set of candidate biomarkers,
flexibly models the effects of continuous biomarkers using free knot
B-splines, and marginalizes over the space of all possible variable
combinations, reducing dependence on one single model specification. At
interim analysis, we assess whether a biomarker-defined subgroup, which may be
the entire trial population, has enhanced or reduced treatment effects,
allowing for early stoppage due to efficacy or futility. At the end of the
trial, the set of candidate tailoring variables is refined using posterior
inclusion probabilities and a reduced model is used for treatment
recommendations. The novelty of our work lies in the combination of flexibly
modelling continuous biomarker-treatment effect relationships using free knot
splines, while simultaneously using BMA to reduce the set of candidate
biomarkers and marginalize over the space of all possible submodels. We
compare our proposed method with two state-of-the-art comparators, Park et al.
(2022) and Liu et al. (2022). Both of these methods out-perform existing
enrichment designs. Comparison with Park et al. (2022) emphasizes the
importance of accounting for nonlinear and nonmonotone relationships between
biomarkers and treatment effectiveness and comparison with Liu et al. (2022)
highlights the effect of variable selection.

The remainder of the paper is organized as follows. Section 2 describes the
notation, models, and interim analysis algorithm associated with the proposed
design. We also describe the computational methods used to derive the
posterior distribution of all unknown parameters. Section 3 describes a
hypothetical trial to detect treatment-sensitive RA patients, motivating the
simulation study in Section 4, in which we assess the operating
characteristics of the proposed design under a variety of settings and compare
its performance to two existing methods. We finish with a discussion in
Section 5.

##  2 Methodology

In this section, we propose analytic methods, appropriate for both Bayesian
adaptive enrichment designs and broader applications, to determine a precision
medicine strategy in the presence of multiple candidate biomarkers, which may
be dichtomous or continuous. To this end, we first define a class of
appropriate models, that flexibly relate the treatment, predictive and
tailoring variables, and their interactions to the study outcome of interest.
We then use BMA to derive a posterior distribution which accounts for the
heterogeneity inherent in selecting an appropriate model. Using the derived
posterior distribution, we can then identify the range of values where
treatment has a positive, clinically important effect.

###  2.1 Notation

For $i=1,\dots,n$, let $Y_{i}$ be the study outcome for individual $i$,
$T_{i}$ be a treatment indicator for individual $i$ ($T_{i}=0$ if assigned to
control, $T_{i}=1$ if assigned to treatment),
$\mathbf{x}_{i}=(x_{1i},\dots,x_{Ji})^{\top}$ be a set of continuous
predictive variables for $Y_{i}$,
$\mathbf{z}_{i}=(z_{1i},\dots,z_{Ri})^{\top}$ be a set of candidate binary
predictive variables for $Y_{i}$, and $\mathbf{\tilde{x}}_{i}$,
$\mathbf{\tilde{z}}_{i}$ be sets of candidate continuous and binary tailoring
variables, respectively. Throughout, tailoring variables refer to those which
interact with treatment to modify treatment effects within specific subgroups,
while a predictive variable directly influences the outcome regardless of
interactions. For both binary and continuous variables, the set of candidate
tailoring variables is constrained to be a subset of the set of candidate
predictive variables. We denote
$\mathcal{D}=\\{(Y_{i},T_{i},\mathbf{x}_{i},\mathbf{z}_{i})\text{,
}i=1,\dots,n\\}$ as the complete set of information available for all
individuals in the study.

###  2.2 A fully Bayesian model specification

In order to identify a subgroup of individuals where treatment is most
effective, we first postulate a fully Bayesian model,

$g(\mathbb{E}\left[Y_{i}\mid\mathbf{x}_{i},\mathbf{z}_{i},T_{i}\right])=\mu+\sum_{j=1}^{J}h_{1,j}(x_{ji})+\boldsymbol{\beta}_{1}^{T}\mathbf{z}_{i}+\left(\phi+\sum_{j=1}^{\tilde{J}}h_{2,j}(\tilde{x}_{ji})+\boldsymbol{\beta}_{2}^{T}\mathbf{\tilde{z}}_{i}\right)T_{i},$
(1)

where $g$ is a pre-specified link function,
$\mathbb{E}\left[Y_{i}\mid\mathbf{x}_{i},\mathbf{z}_{i},T_{i}\right]$ is the
conditional mean of the outcome, $\mu$ is the intercept, $\phi$ is the main
effect of treatment, $h_{1,s}(\cdot)$, $s=1,\dots,J$ and $h_{2,s+J}(\cdot)$,
$s=1,\dots,\tilde{J}$ are free knot splines with pre-specified candidate knots
$K_{s}$ corresponding to the predictive and tailoring effects of the candidate
continuous biomarkers, respectively, and $\boldsymbol{\beta}_{1}$ and
$\boldsymbol{\beta}_{2}$ are the predictive and tailoring effects of the
candidate binary biomarkers, respectively. Additionally, we denote the
coefficients associated with spline $s$ as $\boldsymbol{\theta}_{s}$ and the
complete set of model coefficients as
$\boldsymbol{\zeta}=(\mu,\phi,\boldsymbol{\theta}_{1}^{\top},\dots,\boldsymbol{\theta}_{J}^{\top},\boldsymbol{\beta}_{1}^{\top},\boldsymbol{\theta}_{J+1}^{\top},\dots,\boldsymbol{\theta}_{J+\tilde{J}}^{\top},\boldsymbol{\beta}_{2}^{\top})^{\top}$.
We note that the model formulation described in (1) contains all candidate
tailoring variables, representing the postulated “full model" in the proposed
procedure. Throughout the paper, we refer to the proposed procedure as
“fitting Model (1)"; this refers to the computation of the posterior
distribution averaged over the space of candidate submodels and sub-splines.

We differentiate between model coefficients, terms, and variables. Variables
refer to the biomarkers included in the model. These variables may be
predictive (included in the main effect), tailoring (included in the
interaction effect), or both. Within a model, a variable may be included in up
to two terms: namely, a main effect term and a tailoring effect term. During
the model fitting procedure, terms are iteratively added or removed from the
model. For example, the interaction effect associated with age can be removed,
while the main effect of age remains. Conversely, the main effect of age
cannot be removed without first removing its interaction effect to maintain
well-formulated model hierarchy. A term may consist of one or more
coefficients. While terms associated with binary variables contain a single
coefficient, terms associated with continuous variables generally contain
multiple coefficients that describe the spline function. For example, suppose
we have continuous age and binary sex as candidate predictive and tailoring
variables. Using a 3-degree B-spline with no knots and excluding the intercept
and main effect of treatment, this model has 2 variables (age and sex), 4
associated terms (the main and interaction effect of each), and 8 coefficients
(3 for each spline term, 1 for each binary term).

In order to detect heterogeneous treatment effects, a common goal is to
identify the effective subspace. First, we let
$\boldsymbol{\gamma}(\mathbf{\tilde{x}},\mathbf{\tilde{z}}\mid\phi,\boldsymbol{\theta}_{J+1},\dots,\boldsymbol{\theta}_{J+\tilde{J}},\boldsymbol{\beta}_{2})=\phi+\sum_{j=1}^{\tilde{J}}h_{2,j}(\tilde{x}_{ji})+\boldsymbol{\beta}_{2}^{T}\mathbf{\tilde{z}}$
be the variable-specific treatment effect or the blip effect of treatment
(Robins, 1997); we omit the conditioning and use the shorthand
$\boldsymbol{\gamma}(\mathbf{\tilde{x}},\mathbf{\tilde{z}})$ throughout. Then
the effective subspace can be defined as

$\mathcal{X}^{*}=\left\\{\mathbf{\tilde{x}},\mathbf{\tilde{z}}:P\left(\boldsymbol{\gamma}(\mathbf{\tilde{x}},\mathbf{\tilde{z}})>e_{1}\mid\mathcal{D}\right)>1-\alpha\right\\},$
(2)

where thresholds $\alpha$ and $e_{1}$ are predefined; $\alpha$ can be tuned
through simulation, while $e_{1}$ is typically determined based on clinical
expertise. Given the effective subspace, the hypothesis of a trial designed to
detect a treatment effect in the effective subspace can be written as

$H_{0}:\gamma(\mathbf{\tilde{x}},\mathbf{\tilde{z}})\leq e_{1}\text{
}\forall\text{ }\mathbf{\tilde{x}},\mathbf{\tilde{z}}\in\mathcal{X}\text{ vs.
}H_{A}:\exists\text{
}\mathbf{\tilde{x}},\mathbf{\tilde{z}}\in\mathcal{X}\text{ s.t.
}\gamma(\mathbf{\tilde{x}},\mathbf{\tilde{z}})>e_{1},$

where $\mathcal{X}$ is the complete space of candidate tailoring variables. In
other words, the null hypothesis is that the treatment is ineffective for all
variable combinations and the alternative hypothesis is that there exists at
least one variable combination for which the treatment is effective.

We then assign a prior probability to each submodel. A submodel containing $m$
biomarker terms is assumed to follow a right-truncated Poisson distribution at
$p$, where $p$ is the total number of terms in the full model, excluding the
intercept and the main effect of treatment. This prior probability is given
as:

$p(m\mid\lambda_{1})=\frac{e^{-\lambda_{1}}\lambda_{1}^{m}}{C_{1}m!},$

where $C_{1}$ is a normalization constant. We let $\boldsymbol{\omega}$ be a
vector of indicators of length $p$ denoting which biomarker terms are included
in a given submodel. In the absence of prior information, one can assume that
all biomarker terms have equal probability of being included in a given
submodel, such that $p(\boldsymbol{\omega}\mid m)={m\choose p}^{-1}$.
Similarly, we assign a prior probability to the number of knots in each
spline, $k_{s}$:

$p(k_{s}\mid\lambda_{2})=\frac{e^{-\lambda_{2}}\lambda_{2}^{k_{s}}}{C_{2}k_{s}!},$

where $C_{2}$ is a normalization constant. All knot locations from a set of
candidate knots of size $K_{s}$ have equal probability such that
$p(\boldsymbol{\kappa}_{s}\mid k_{s})={K_{s}\choose k_{s}}^{-1}$, where
$\boldsymbol{\kappa}_{s}$ is a vector of length $K_{s}$ of indicator variables
denoting which knots are used for a given spline. Both $\lambda_{1}$ and
$\lambda_{2}$ must be specified a priori; larger $\lambda_{1}$ favors larger
models and larger $\lambda_{2}$ favors more complex splines. The model
coefficients, $\boldsymbol{\zeta}$, are assumed to have independent Normal
priors, centered at $0$ with variance $\sigma_{B}^{2}$, where $\sigma_{B}^{2}$
is selected such that the prior is weakly informative.

In our simulation study, $Y_{i}$ is Normally distributed such that
$Y_{i}=\mathbb{E}\left[Y_{i}\mid\mathbf{x}_{i},\mathbf{z}_{i},T_{i}\right]+\tau_{i}$,
where $\tau_{i}\sim\mathcal{N}(0,\sigma_{\tau}^{2})$. We postulate a prior
distribution on $\sigma_{\tau}$:

$\sigma_{\tau}\sim\mathcal{IG}(a_{0},b_{0}),$

where $\mathcal{IG}$ is an Inverse-Gamma distribution.

###  2.3 Computation

In this section, we propose computational methods to derive the posterior
distribution of all unknown parameters in Model (1). Due to the complicated
nature of the model, a closed-form analytic expression for the joint posterior
distribution is unavailable and we must use a Markov Chain Monte Carlo (MCMC)
sampler. To compute the posterior distribution, we simultaneously average over
the space of knot configurations for each spline term and average over the
space of all combinations of appropriate biomarker terms. Both of these
elements are independently challenging as the MCMC sampler must move across
spaces of varying dimension. Reversible jump (rj) MCMC, a commonly used
sampler in BMA, facilitates efficient exploration of model spaces with varying
dimensions (Chen et al., 2011).

The proposed rjMCMC sampling procedure proceeds as follows (see Web Appendix A
for more detail). Unlike MH, which operates within a fixed parameter space and
model structure, rjMCMC dynamically explores different model spaces by
proposing changes in both model parameters and model dimensionality. First,
the procedure initializes parameter states. It then alternates between
proposing changes in model structure using reversible jump moves (e.g.,
changing the knots for each spline term and adding or removing biomarker terms
from the model), updating model parameters and spline coefficients within the
current model using MH steps, and updating the individual-level heterogeneity
parameter using a Gibbs step. In each step, a proposed state in drawn from a
proposal distribution and accepted with a given probability. For reversible
jump steps, the acceptance probability depends on the proposed state, the
current state, and the determinant of the Jacobian matrix of the
transformation to account for the change in dimensionality and ensure detailed
balance. If the proposed state is accepted, the algorithm moves to the new
state; otherwise, it stays in the current state. The process iterates through
these steps for a large number of iterations, generating a sequence of states
that characterize the posterior distribution.

More specifically, updating the set of knots for each spline term currently in
the model proceeds as follows. To facilitate the description, we use the term
active set and inactive set to refer to the set of knots or biomarker terms
currently in and not in the model, respectively. The first step is a knot
position change in which we propose a new location for a randomly selected
knot, while keeping the number of knots constant. As this does not require a
dimension change, we can use a traditional MH acceptance probability. Next, a
stochastic decision is made to either propose the addition of a knot (birth
step) or the removal of a knot (death step). This is a reversible jump step,
in which the dimension of the proposed parameter space differs from the
dimension of the current parameter space. For the birth step, a randomly-
selected knot from the inactive set is proposed to be be added. This involves
(1) stochastically generating a proposed model coefficient for the new knot,
and (2) augmenting the proposal for other model coefficients associated with
the spline term to account for this addition. The other coefficients must be
augmented as changing the knot locations will change the interpretation of all
coefficients in that spline term. Various methods have been proposed to deal
with this challenging problem; we use an adaptation by Moore et al. (2020) of
methods first suggested by Gamerman (1997). Similarly, the death step is
simply a reverse of the birth step.

The next reversible jump step of the proposed procedure facilitates the
averaging of posterior distributions of parameters across the complete space
of candidate biomarkers. Here, we randomly decide between proposing to add a
biomarker term (birth step) or to remove a biomarker term (death step). Rather
than adding or removing one model coefficient at a time as in traditional BMA,
we perform operations on the terms themselves; namely, the removal or addition
of a spline term will generally correspond to a change in dimension greater
than one. Further, the terms to be added or removed must be carefully selected
to ensure well-formulated hierarchy is maintained. For example, if the current
submodel does not contain the main effect of age, its interaction effect is
not eligible to be added. It is only after the main effect is in the active
set of a submodel that we can add the corresponding interaction effect.
Similarly, we cannot remove a main effect term without removing its
interaction effect.

For the birth step, we select an eligible term from the inactive set and
propose its addition to the model. Once the eligible term is selected, we
propose updates to the associated model coefficients. This involves (1)
stochastically generating one or more proposed coefficients associated with
the added term, and (2) proposing adjustments to the other coefficients to
account for this addition. While (2) is not necessary in general BMA, our
findings suggest that such adjustments significantly enhance acceptance rates,
particularly when adding a term may correspond to the addition of multiple
coefficients. The death step is simply the reverse of the birth step.

###  2.4 Adaptive and stopping decision rules

At each interim analysis, we propose using the posterior distribution of the
treatment effect averaged over the enriched population to define our decision
criteria. We denote by $\mathcal{D}_{\ell}$ the data accumulated up to interim
analysis $\ell$; we will use this information to make decisions to stop the
trial early for superiority or futility and to adaptively enrich the trial
population. We propose the following procedure at each interim analysis
$\ell=1,\dots,L-1$:

  1. 1.

Identify the effective subspace: Fit (or re-fit) Model (1) using the
reversible jump sampling procedure specified in Section 2.3. The model fitting
procedure is based on the data observed up to the given interim analysis
point, $\mathcal{D}_{\ell}$. Identify the effective subspace,
$\mathcal{X}_{\ell}^{*}$, using the definition provided in (2), with
$\mathcal{D}$ replaced by $\mathcal{D}_{\ell}$. Recalling that
$\boldsymbol{\gamma}(\mathbf{\tilde{x}},\mathbf{\tilde{z}})=\phi+\sum_{j=1}^{\tilde{J}}h_{2,j}(\tilde{x}_{ji})+\boldsymbol{\beta}_{2}^{T}\mathbf{\tilde{z}}$
represents the subgroup-specific treatment effect, the effective subspace at
the $\ell$-th interim analysis is

$\mathcal{X}_{\ell}^{*}=\left\\{\mathbf{\tilde{x}},\mathbf{\tilde{z}}:P\left(\boldsymbol{\gamma}(\mathbf{\tilde{x}},\mathbf{\tilde{z}})>e_{1}\mid\mathcal{D}_{\ell}\right)>1-\alpha\right\\}.$

Then compute the prevalence of the effective subspace, defined by the
percentage of individuals in the trial who have biomarker values in
$\mathcal{X}^{*}$. If this prevalence is less than some cutoff, $\pi$, stop
the trial for futility. We then compute the enriched treatment effect as

$\Delta_{\ell}=\int_{\mathcal{X}_{\ell}^{*}}\boldsymbol{\gamma}(\mathbf{\tilde{x}},\mathbf{\tilde{z}})\hat{f}(\mathbf{\tilde{x}},\mathbf{\tilde{z}})\partial\mathbf{\tilde{x}}\partial\mathbf{\tilde{z}},$
(3)

where $\hat{f}$ is the empirical probability distribution function of the
tailoring variables.

  2. 2.

Stop for efficacy: Using the posterior predictive probability of trial success
in the effective subspace, assess if the trial should be stopped for efficacy
using the following criterion:

$P\left(\Delta_{\ell}>b_{1}\mid\mathcal{D}_{\ell}\right)>B_{1}.$ (4)

  3. 3.

Stop for futility: Using the posterior predictive probability of trial failure
in the effective subspace, assess if the trial should be stopped for futility
using the following criterion:

$P\left(\Delta_{\ell}<b_{2}\mid\mathcal{D}_{\ell}\right)>B_{2}.$ (5)

Figure 1 summarizes the procedure. If the trial is not stopped early, at the
last interim analysis, we conclude that the treatment is effective in the
sensitive population using the efficacy criterion (4) for $\ell=L$. Otherwise,
we conclude that treatment is not superior in the sensitive population. This
design accounts for the possibility that the entire trial population is in the
effective subspace. In this case, we assess the efficacy of treatment in the
entire trial population, as in a traditional randomized clinical trial, and
$\Delta_{\ell}$ becomes the classic average treatment effect.

Figure 1: A visual representation of the interim decision rules

We additionally incorporate adaptive enrichment after each interim analysis
step. We use $\mathcal{X}_{\ell}^{*}$, defined in (3), to define the inclusion
criteria for individuals enrolled after the $\ell$-th interim analysis. That
is, after interim analysis $\ell$ only patients with variables
$\mathbf{\tilde{z}},\mathbf{\tilde{x}}\in\mathcal{X}_{\ell}^{*}$ will be
enrolled. In addition to improving the probability of enrolling patients who
will benefit from treatment, we anticipate that this adaptive enrichment will
improve study power. The procedure is as follows. Before the first interim
analysis, there is no adaptive enrichment. After the first interim analysis,
we use the accumulated data to restrict entry of patients to those expected to
be in the effective subspace, updating the criteria at each interim time
period.

To ensure that the effective subgroup is statistically meaningful, we may
impose restrictions on the prevalence of the effective subgroup in the study.
If only a small percentage of patients benefit from the treatment, it may not
be useful in practice. To prevent this, we can additionally impose that the
trial should be stopped for futility if the prevalence of the effective
subspace falls below some proportion $\pi$ between 0 and 1. Park et al. (2022)
recommend using values of $\pi$ between $0.01-0.10$; we use $\pi=0.10$ in the
simulation study.

Throughout the trial, we avoid explicitly pruning the list of candidate
variables as BMA will not necessarily completely drop variables which have
little to no relation to the outcome. Instead, irrelevant variables will be
excluded from most of the models in the sampler and have a low posterior
inclusion probability. We propose to prune variables with a low posterior
inclusion probability (say, $<10\%$) at the last interim analysis, or at the
end of the trial. This approach of relying on fewer biomarkers for treatment
recommendations not only simplifies clinical implementation but may also
reduce costs, particularly if some biomarkers are costly to measure. It is
entirely possible that all tailoring variables will be trimmed, which simply
corresponds to a traditional, overall analysis. Finally, Model (1) can be re-
fit excluding the irrelevant variables to remove unnecessary noise in
preparation for clinical recommendations.

In practice, the cutoffs ($b_{1}$, $b_{2}$, $B_{1}$, $B_{2}$) and the number
and timing of interim analyses must be specified. The values of $b_{1}$ and
$b_{2}$ are typically pre-specified by clinicians and correspond to clinically
meaningful effect sizes. The values of $B_{1}$ and $B_{2}$ are generally tuned
in design simulation studies and will depend on factors such as the number of
candidate variables and their anticipated effect sizes, desired type I error
rate, power, and the number of planned interim analyses. Generally, logistical
considerations prevent more than three interim analyses to be conducted. A
reasonable rule-of-thumb is to conduct the first interim analysis after $1/3$
or $1/2$ of the maximum sample size has been enrolled (Park et al., 2022).

##  3 Application to rheumatoid arthritis

We implement the proposed trial design to identify the subset of RA patients
who are expected to benefit from anti-tumor necrosis factor (anti-TNF)
treatment. The hypothetical trial will enroll a diverse cohort of RA patients,
including those who have previously failed alternative anti-TNF therapies in
the past. Patients will undergo comprehensive molecular profiling using multi-
omics technologies before initiating anti-TNF treatment and at specified
intervals during the treatment course. Clinical outcomes, including disease
activity scores, radiographic scores, and patient-reported outcomes, will be
closely monitored throughout the trial. Statistical methods, such as fitting
proposed Model (1), will be used to identify treatment-sensitive subgroups and
the associated robust biomarkers. By integrating molecular data with clinical
outcomes, this trial design aims to pave the way for personalized medicine
approaches in managing RA, ensuring that patients receive the most effective
therapies tailored to their individual biological profiles.

A meta-analysis of 41 studies investigating 99 different laboratory biomarkers
revealed inconsistent associations with anti-TNF treatment response in RA
patients (Wientjes et al., 2022). Notably, among the five biomarkers assessed
in more than two studies—anti-cyclic citrullinated peptide (anti-CCP),
rheumatoid factor (RF), -308 polymorphism in the TNF-$\alpha$ gene, presence
of one or two shared epitope (SE) copies in the HLA-DRB1 gene, and FcGR2A
polymorphism (rs1081274)—mixed results were observed. A weak association was
observed for three markers studied less frequently: granulocyte-macrophage
colony-stimulating factor (GM-CSF), interleukin (IL)-34 concentration, and the
combined biomarker of serum IL-6 and survivin level. Building on these
insights, our proposed trial design includes several clinical biomarkers, such
as sex, body mass index, smoking status, alongside candidate binary and
continuous molecular biomarkers. Section 4 describes simulation studies to
determine the operating characteristics of such a trial design. In Simulation
Study I, we consider five binary biomarkers (representing, for example, sex,
smoking status, RF, anti-CCP, SE copies) and one continuous marker
(representing, for example, IL-34 concentration). In Simulation Study II, we
consider two continuous biomarkers, such as level of GM-CSF and IL-34
concentration.

##  4 Simulation study

In this section, we summarize the results of two simulation studies to assess
the operating characteristics of the proposed Free Knot-Bayesian Model
Averaging (FK-BMA) compared with two existing methods within the context of
the RA trial. We assumed a maximum sample size of 500 with one interim
analyses at 300 patients and considered candidate biomarkers which are both
continuous and dichotomous. In order to define the effective subspace, we set
various subsets of the variables to have a non-zero interaction with
treatment; we refer to biomarkers with a non-zero interaction effect as the
true tailoring variables.

We compared the proposed FK-BMA with two existing methods: Park, Liu, Thall,
and Yuan [PLTY (Park et al., 2022)] in Simulation Study I and Liu, Kairalla,
and Renfro [LKR (Liu et al., 2022)] in Simulation Study II. The simulations
were divided into two studies, one with mixed biomarker type and another with
solely continuous biomarkers, to facilitate evaluation of LKR, which can only
handle continuous biomarkers. While PLTY can be used when tailoring variables
are continuous, performance is poor in such settings and thus we excluded its
evaluation from Study I.

We used the fitting procedures from either method, adapted to continuous
outcomes, with the interim decision rules described in Section 2.4. In PLTY,
spike-and-slab priors are used for variable selection and each continuous
biomarker is assumed to have a linear predictive relationship with outcome and
a linear tailoring relationship with treatment. Additionally, well-formulated
model hierarchy is maintained by first determining if the main effect should
be included; if so, the interaction effect follows suit; otherwise, the
interaction effect is excluded. This implies that the inclusion probability of
interaction effects directly depends on their corresponding main effects’
inclusion probability. Thus, the comparison between our method and PLTY
emphasizes the benefits of accounting for nonlinear predictive and tailoring
effects for continuous biomarkers and provides insights into the relative
effectiveness of BMA and spike-and-slab priors as approaches for variable
selection. To facilitate a fair comparison, we forced inclusion of the main
effect of treatment in PLTY. In LKR, penalized splines are used to flexibly
model the relationships between biomarkers and treatment and no variable
selection is performed. The comparison between our method and LKR highlights
the differences between the two spline estimation techniques, namely free-knot
(FK-BMA) and penalized splines (LKR), and the benefits of variable selection
in FK-BMA.

We assessed eight simulation scenarios for each study, inspired by the
hypothetical RA trial, described in Table 1. In Simulation Study I, we
considered six candidate biomarkers as both predictive and tailoring
variables: one continuous biomarker, $x_{i}$, and five binary biomarkers,
$\mathbf{z}_{i}=(z_{1i},\dots,z_{5i})$ with prevalence 0.35, 0.50, 0.65, 0.20,
and 0.35, respectively. Data were generated from

$Y_{i}=h_{1}(x_{i})+\boldsymbol{\beta}_{1}^{\top}\mathbf{z}_{i}+\left(\phi+h_{2}(x_{i})+\boldsymbol{\beta}_{2}^{\top}\mathbf{z}_{i}\right)T_{i}+\tau_{i},$
(6)

for $i=1,\dots,n$ and $\tau_{i}\sim N(0,1)$ and
$\gamma(x,\mathbf{z})=\phi+h_{2}(x)+\boldsymbol{\beta}_{2}^{\top}\mathbf{z}$.
In Simulation Study II, we considered two continuous biomarkers
$\mathbf{x}_{i}=(x_{1i},x_{2i})$. For $i=1,\dots,n$, the simulated individual-
level data were generated from:

$Y_{i}=0.5x_{1i}+\left(\phi+h_{2,1}(x_{1i})+h_{2,2}(x_{2i})\right)T_{i}+\tau_{i},$

where $\tau_{i}\sim N(0,1)$ and
$\gamma(x_{1},x_{2})=\phi+h_{2,1}(x_{1})+h_{2,2}(x_{2})$. All continuous
biomarkers were uniformly distributed over the interval $(0,1)$. Additional
details on simulation parameters and plots are available in Web Appendix B:
the true values for the predictive effects are shown in Table S1 and visual
representations of scenarios with true continuous tailoring biomarkers are
shown in Figures S1 and S2.

Table 1: True model parameters used for Simulation Study I and II. $\Delta$ is the true average treatment effect in the effective subspace. | Tailoring Variables | Effective Subspace
---|---|---
Scenario | Type | Dimension | $\gamma(\mathbf{x},\mathbf{z})$ | Prevalence | $\Delta$
Simulation Study I
1 | - | 0 | 0 | 0 | 0
2 | - | 0 | 0.28 | 1 | 0.28
3 | Binary | 1 | $z_{1}-0.3$ | 0.35 | 0.70
4 | Binary | 1 | $0.7z_{2}-0.14$ | 0.50 | 0.56
5 | Binary | 1 | $0.8z_{3}-0.3$ | 0.65 | 0.50
6 | Binary | 2 | $0.9z_{4}+0.9z_{1}-0.2$ | 0.48 | 0.81
7 | Continuous | 1 | $2.3x-1.15$ | 0.50 | 0.58
8 | Continuous | 1 | $\text{cos}(2\pi x)$ | 0.50 | 0.64
Simulation Study II
1 | - | 0 | 0 | 0 | 0
2 | - | 0 | 0.35 | 1 | 0.35
3 | Continuous | 1 | $2.3x_{1}-1.15$ | 0.50 | 0.58
4 | Continuous | 1 | $\text{cos}(2\pi x_{1})$ | 0.50 | 0.64
5 | Continuous | 1 | $\frac{1.4\text{exp}\\{25(x_{1}-0.5)\\}}{1+\text{exp}\\{25(x_{1}-0.5)\\}}-0.6$ | 0.51 | 0.71
6 | Continuous | 1 | If $x_{1}\leq 0.5$, $\frac{2\text{exp}\\{30(x_{1}-0.3)\\}}{1+\text{exp}\\{30(x_{1}-0.3)\\}}-1$ | 0.40 | 0.77
|  |  | If $x_{1}>0.5$, $\frac{2}{1+\text{exp}\\{30(x_{1}-0.7)\\}}-1$ |  |
7 | Continuous | 1 | If $x_{1}\leq 0.5$, $\frac{1.5}{1+\text{exp}\\{30(x_{1}-0.3)\\}}-0.75$ | 0.60 | 0.64
|  |  | If $x_{1}>0.5$, $\frac{1.5\text{exp}\\{30(x_{1}-0.7)\\}}{1+\text{exp}\\{30(x_{1}-0.7)\\}}-0.75$ |  |
8 | Continuous | 2 | $2.3x_{1}+\text{cos}(2\pi x_{2})-1.15$ | 0.50 | 0.79

For both simulation studies, we set $e_{1}=b_{1}=b_{2}=0$ and $\pi=0.10$ so
that the trial is stopped if less than $10\%$ of patients are sensitive to
treatment. We set the overall type I error rate for FK-BMA to be $0.05$ and
tuned the design parameters accordingly ($B_{1}=0.975$; $B_{2}=0.8$;
$\alpha=0.3$ for I and $\alpha=0.2$ for II). It is important to note that the
trial parameters were specifically tuned for FK-BMA to facilitate a direct
comparison of the analysis models. In practical applications, one should
adjust the trial parameters according to the specific analysis model employed.
We additionally tuned the MCMC parameters in each sampler for optimal
performance. In Simulation Study I, we set the parameters $\lambda_{1}=0.1$,
$\lambda_{2}=1$, and $\sigma_{B}=10$ for the FK-BMA method, tuned based on
simulation. For PLTY, we tuned via simulation, as per their established
notation, sparsity parameters $\mu$ and $\tau$ (Web Appendix C, Table S2). In
Simulation Study II, we set $\lambda_{1}=0.01$, $\lambda_{2}=1$, and
$\sigma_{B}=10$ for the FK-BMA approach and used 5 internal knots placed at
the $(q/6)$-quantiles ($q=1,\dots,5$) for each continuous candidate biomarker
in LKR method, tuned by simulation. The candidate internal knots for FK-BMA
were placed at the $(q/10)$-quantiles ($q=1,\dots,9$) for each continuous
candidate biomarker. We generated 1000 simulated trials in each scenario
considered. Web Appendix D describes the MCMC configuration settings used in
the simulation study and diagnostics to assess convergence in our proposed
procedure.

As operating characteristics, we considered traditional type I error, power,
generalized power to detect subgroup-specific effects, success rate in
correctly identifying the appropriate tailoring variables, and accuracy of
treatment recommendations using a large, external testing dataset. More
specifically, we define generalized power, similar to that of Chapple and
Thall (2018), as the empirical probability that the following two events
occur: (1) the efficacy criteria is met (4) in either an interim analysis step
or the final analysis, and (2) the correct tailoring variables are identified.
To assess accuracy of treatment recommendations, we generated an external
dataset with 10,000 individuals and their biomarker information at the start
of the simulation study. In each simulation iteration, we used the simulated
trial results to recommend a treatment plan for individuals in the external
dataset. Since the true functional relationship between biomarkers and the
effect of treatment is known, we then calculated the percentage of individuals
for which this treatment recommendation was optimal.

Table 2: Results from 1000 simulated datasets in Simulation Study I for the
proposed FK-BMA and the PLTY approach.

| Power | Generalized Power | Correct Marker | Accuracy | Size
---|---|---|---|---|---
Scenario | FK-BMA | PLTY | FK-BMA | PLTY | FK-BMA | PLTY | FK-BMA | PLTY | FK-BMA | PLTY
|  $\mu=100$, $\tau=0.05$
1 | 0.04 | 0.05 | - | - | 1.00 | 0.91 | 0.83 | 0.81 | 349.20 | 351.00
2 | 0.86 | 0.89 | 0.86 | 0.79 | 1.00 | 0.90 | 0.96 | 0.97 | 361.60 | 360.60
3 | 0.88 | 0.70 | 0.87 | 0.66 | 0.96 | 0.92 | 0.95 | 0.88 | 375.00 | 413.60
4 | 0.89 | 0.93 | 0.73 | 0.84 | 0.79 | 0.91 | 0.80 | 0.83 | 369.20 | 388.20
5 | 0.92 | 0.93 | 0.80 | 0.86 | 0.86 | 0.92 | 0.89 | 0.92 | 367.40 | 383.20
6 | 1.00 | 0.93 | 0.73 | 0.65 | 0.73 | 0.70 | 0.82 | 0.56 | 324.40 | 370.80
7 | 0.98 | 0.37 | 0.98 | 0.36 | 1.00 | 0.94 | 0.92 | 0.66 | 322.60 | 389.00
8 | 0.98 | 0.04 | 0.98 | 0.02 | 1.00 | 0.21 | 0.89 | 0.50 | 322.82 | 359.00
|  $\mu=50$, $\tau=0.1$
1 | 0.03 | 0.04 | - | - | 1.00 | 0.98 | 0.80 | 0.79 | 352.80 | 351.80
2 | 0.89 | 0.90 | 0.89 | 0.88 | 1.00 | 0.98 | 0.97 | 0.97 | 356.80 | 359.60
3 | 0.86 | 0.58 | 0.86 | 0.57 | 0.96 | 0.98 | 0.95 | 0.82 | 372.00 | 419.40
4 | 0.88 | 0.81 | 0.72 | 0.79 | 0.78 | 0.98 | 0.80 | 0.70 | 372.80 | 411.80
5 | 0.92 | 0.84 | 0.80 | 0.82 | 0.86 | 0.98 | 0.90 | 0.82 | 364.00 | 407.60
6 | 1.00 | 0.87 | 0.70 | 0.75 | 0.71 | 0.86 | 0.82 | 0.52 | 322.20 | 382.40
7 | 0.98 | 0.33 | 0.98 | 0.33 | 1.00 | 0.97 | 0.92 | 0.65 | 322.82 | 393.00
8 | 0.98 | 0.03 | 0.98 | 0.00 | 0.99 | 0.07 | 0.88 | 0.50 | 318.02 | 358.60

Results from Simulation Study I are presented in Table 2 for two sets of PLTY
sparsity parameters. The tuning parameters for FK-BMA are held constant; any
observed discrepancies can be attributed to simulation variation. In Scenario
1, both methods control type I error; in Scenario 2 (constant, positive
treatment effect), they perform similarly across metrics. In scenarios with
one or more binary tailoring variables (3-6), PLTY has similar or higher rates
of correct marker detection, while FK-BMA has similar or higher accuracy
rates. This is emphasized in Scenario 6, where the accuracy of PLTY is
0.52-0.56, compared with 0.82 for FK-BMA. Generalized power results differ in
settings with a single binary tailoring variable. Specifically, in Scenario 3
with lower prevalence (0.35), FK-BMA demonstrates superior generalized power;
the reverse is true in Scenarios 4-5 with higher prevalence (ranging from 0.5
to 0.65). In Scenario 6 with two binary tailoring variables, PLTY exhibits
higher generalized power for the first set of sparsity parameters, while the
second set favors FK-BMA. In Scenarios 7-8, there is a single continuous
tailoring variable with either a linear (7) or nonlinear (8) relationship with
treatment effect. With either set of sparsity parameters, the proposed FK-BMA
out-performs PLTY on all metrics considered. Despite high rates of correct
biomarker detection in Scenario 7, the power and accuracy of PLTY remain low.
As expected, when trends are nonlinear (8), performance is poor across all
metrics for PLTY as its assumptions are not met.

Table 3: Results from 1000 simulated datasets in Simulation Study I (with no
predictive biomarkers, i.e., $h_{1}(x)=0$ and
$\boldsymbol{\beta}_{1}=\mathbf{0}$) for the proposed FK-BMA and the PLTY
approach.

| Power | Generalized Power | Correct Marker | Accuracy | Size
---|---|---|---|---|---
Scenario | FK-BMA | PLTY | FK-BMA | PLTY | FK-BMA | PLTY | FK-BMA | PLTY | FK-BMA | PLTY
1 | 0.04 | 0.04 | - | - | 1.00 | 0.89 | 0.81 | 0.76 | 352.80 | 365.80
2 | 0.88 | 0.89 | 0.88 | 0.79 | 1.00 | 0.88 | 0.97 | 0.97 | 362.20 | 360.60
3 | 0.82 | 0.15 | 0.81 | 0.04 | 0.92 | 0.14 | 0.93 | 0.55 | 374.00 | 381.40
4 | 0.84 | 0.66 | 0.59 | 0.06 | 0.67 | 0.10 | 0.75 | 0.50 | 385.20 | 398.40
5 | 0.90 | 0.68 | 0.71 | 0.08 | 0.76 | 0.12 | 0.86 | 0.62 | 373.20 | 394.60
6 | 0.99 | 0.86 | 0.55 | 0.01 | 0.56 | 0.01 | 0.77 | 0.48 | 339.40 | 367.20
7 | 0.98 | 0.17 | 0.98 | 0.15 | 1.00 | 0.77 | 0.92 | 0.55 | 323.60 | 372.80
8 | 0.98 | 0.04 | 0.98 | 0.00 | 1.00 | 0.03 | 0.88 | 0.50 | 317.20 | 355.00

Recall that the comparator PLTY establishes well-formulated model hierarchy by
linking the inclusion probability of interaction terms with that of the
corresponding main effect terms. Specifically, if a main effect term is absent
from the submodel, the interaction effect is automatically excluded. To
evaluate the robustness of this approach, we conducted simulations where there
are no biomarker effects observed in the control group, i.e., $h_{1}(x)=0$ and
$\boldsymbol{\beta}_{1}=\mathbf{0}$ in (6), using PLTY sparsity parameters
$(\mu,\tau)=(100,0.05)$ (Table 3). In FK-BMA, the performance declined
compared with Table 2 in the case of true binary tailoring variables
(Scenarios 3-6); otherwise, the results are similar. For most of the scenarios
with subgroup differences (3-6, 8), the rates of detecting the correct
biomarker are low, corresponding to low generalized power and low accuracy in
PLTY. In Scenario 7, the rate of detecting the correct biomarker is high
(0.77), while power is low (0.17). As expected, performance of PLTY remains
adequate when there are no subgroup specific effects, as the set of tailoring
variables (none) matches the set of predictive variables (none). Hence,
ensuring a correspondence between predictive and tailoring variables is
essential to achieve satisfactory performance in PLTY.

Table 4: Results from 1000 simulated datasets in Simulation Study II for the
proposed FK-BMA and the LKR approach.

| Power | Generalized Power | Correct Marker | Accuracy | Size
---|---|---|---|---|---
Scenario | FK-BMA | LKR | FK-BMA | LKR | FK-BMA | LKR | FK-BMA | LKR | FK-BMA | LKR
1 | 0.04 | 0.24 | - | - | 1.00 | 0.00 | 0.87 | 0.86 | 334.20 | 375.83
2 | 0.96 | 0.99 | 0.96 | 0.00 | 1.00 | 0.00 | 0.99 | 0.79 | 329.20 | 310.28
3 | 0.93 | 1.00 | 0.93 | 0.00 | 0.99 | 0.00 | 0.89 | 0.91 | 325.80 | 301.84
4 | 0.94 | 1.00 | 0.93 | 0.00 | 0.99 | 0.00 | 0.85 | 0.89 | 327.43 | 300.60
5 | 0.94 | 1.00 | 0.93 | 0.00 | 0.98 | 0.00 | 0.92 | 0.94 | 334.20 | 300.62
6 | 0.95 | 1.00 | 0.95 | 0.00 | 0.99 | 0.00 | 0.93 | 0.94 | 310.03 | 300.41
7 | 0.88 | 1.00 | 0.84 | 0.00 | 0.90 | 0.00 | 0.81 | 0.88 | 353.16 | 301.21
8 | 0.99 | 1.00 | 0.98 | 1.00 | 0.99 | 1.00 | 0.85 | 0.89 | 308.32 | 300.00

Simulation Study II results, depicted in Table 4, reveal LKR’s notably
inflated type I error (0.24) and near-maximal power in non-null scenarios.
Generalized power is zero for LKR in all but Scenario 8, as it does not
perform tailoring variable selection. For FK-BMA, the type I error rate is
0.04, with lower power. Further, FK-BMA selects the correct biomarker over 90%
of the time in all scenarios considered. LKR is more sensitive to subgroup
differences due to its use of a saturated, flexible model, meaning it often
detects subgroup differences where they don’t exist, such as in Scenario 2.
Here, LKR concludes there are subgroup differences 90% of the time (not
shown), leading to a decrease in accuracy compared with FK-BMA (0.79 vs 0.99).
In scenarios with subgroup-specific effects, FK-BMA demonstrates slightly
lower but comparable accuracy rates to LKR, indicating a minimal loss of
predictive performance from FK-BMA’s model simplification.

##  5 Discussion

In this work, we propose Bayesian adaptive enrichment methods that effectively
identify important tailoring variables from a candidate pool, flexibly model
continuous biomarker effects using free knot B-splines, and reduce dependence
on a single model specification by marginalizing over all possible variable
combinations. Our proposed approach is the first of its kind to combine BMA
with free-knot splines, leveraging the strength of both approaches. BMA
integrates model uncertainty, leading to more robust inference and mitigating
overfitting risks. Free knot B-splines automatically determine knot number and
positions, offering a data-driven, flexible modeling approach without
subjective bias. Additionally, we propose a design that incorporates adaptive
enrichment and early stopping of treatment arms, potentially reducing trial
costs and improving patient outcomes.

Compared with existing approaches, our proposed Bayesian adaptive enrichment
design has desirable operating characteristics. Our proposed approach
demonstrates comparable performance to Park et al. (2022) in detecting binary
tailoring variables. Notably, our method exhibits enhanced operating
characteristics in scenarios involving continuous tailoring variables,
including cases with nonlinear and nonmonotone effects where the assumptions
of Park et al. (2022) are not met, as well as in instances with simple linear
effects. Additionally, our method maintains good performance in settings where
a biomarker shows no association with the outcome in the control group. Thus,
either our method or that of Park et al. (2022) is appropriate in settings
where (1) all candidate biomarkers are binary and (2) the set of true
predictive biomarkers is expected to exactly match the set of true tailoring
biomarkers. If one or more of these conditions is not met, our approach is
more appropriate. We also compared our method with Liu et al. (2022), which
uses penalized splines without variable selection. Liu et al. (2022) is more
sensitive to minor subgroup differences, leading to higher type I error rates
and increased power. Tuning trial parameters to Liu et al. (2022) would
mitigate the type I error rate concern, making it suitable for settings with
predetermined tailoring variables or a single candidate tailoring variable.
However, our method performs well across all settings, evidenced by its
comparable accuracy and high rates of correct biomarker detection.

A practical limitation of our proposed design is that it is highly
computationally intensive, particularly when the number of candidate
biomarkers is large. Specifically, when using binary biomarkers, computational
time is reduced due to simplified model fitting. However, adding continuous
markers with free-knot splines, increases computational complexity,
underscoring the trade-offs between biomarker types and computational
efficiency. Additionally, estimating quantities such as the empirical
probability density function can be challenging within this context. Future
work will explore the consequences of high dimensional variable spaces and
investigate alternative sampling algorithms, such as reversible jump piece-
wise deterministic processes, which may be more computationally efficient
(Chevallier et al., 2023). Moreover, the effective subspace may not lead to
straightforward thresholds, as we did not enforce monotonicity constraints on
spline terms. To address this, methods like Bayesian additive regression trees
can be employed to derive simplified cutoffs suitable for clinical
applications. Additionally, we assume that the primary endpoint of the trial
can be observed relatively quickly compared to the enrollment rate, ensuring
that interim analyses and adaptations are well-informed by early outcomes.
However, if the primary endpoint is not immediately observable, one can use
early response outcomes as in Park et al. (2022). Future work will focus on
adapting our methods to incorporate this proxy outcome approach. Additionally,
this design should only be used in settings with a large enough sample size to
detect clinically meaningful benefits. Conducting simulation studies specific
to the disease context in which the trial will be conducted is crucial to
determine the feasibility of such an approach.

Although our design is for a single treatment comparison, it can be extended
to compare multiple treatments. Following the methods proposed by Chapple and
Thall (2018), one option is to introduce a latent treatment variable with
strategic priors to allow treatments to be flexibly combined throughout the
MCMC procedure, allowing for efficiency gains without sacrificing flexibility.
Future work will formalize the model specification and assess the performance
of the proposed approach in this more complex setting.

In sum, personalized medicine offers a promising avenue to improve patient
outcomes and reduce healthcare costs. This dynamic field necessitates new ways
of approaching patient recruitment and trial design, particularly in domains
like RA where tailoring biomarkers and their functional forms are largely
unknown but are the focus of intense investigation. To this end, we offer
innovative analysis methods, computational tools, and a trial design
framework, presenting an important pathway to address these evolving needs.

## Acknowledgements

Shirin Golchi is a Fonds de Recherche du Québec, Santé, Chercheuse-boursiere
(Junior 1) and acknowledges support from the Natural Sciences and Engineering
Council of Canada (NSERC-Discovery Grant), Canadian Statistical Sciences
Institute and Fonds de Recherche du Québec, Nature et technologies (FRQNT-
NSERC NOVA). Erica E.M. Moodie is a Canada Research Chair (Tier 1) in
Statistical Methods for Precision Medicine and acknowledges the support of a
Chercheur de Mérite Career Award from the Fonds de Recherche du Québec, Santé.
Research reported in this publication was supported in part by the National
Institute of Mental Health of the National Institutes of Health under Award
Number R01 MH114873 and the McGill Interdisciplinary Initiative in Infection
and Immunity.

## Supplementary Materials

Web Appendix A-D, referenced in Sections 2.3 and 4, are available with this
paper at the Biometrics website on Wiley Online Library. The R code to
implement the proposed methods can be found at
https://github.com/laramaleyeff1/Bayesian_free_knot_model_average_adaptive_design.

## References

  * Adams and Brantner (2006) Adams, C. P. and Brantner, V. V. (2006).  Estimating the cost of new drug development: is it really $802 million?  Health Affairs 25, 420–428.
  * Asano and Hirakawa (2020) Asano, J. and Hirakawa, A. (2020).  A Bayesian basket trial design accounting for uncertainties of homogeneity and heterogeneity of treatment effect among subpopulations.  Pharmaceutical Statistics 19, 975–1000.
  * Brannath et al. (2009) Brannath, W., Zuber, E., Branson, M., Bretz, F., Gallo, P., Posch, M., and Racine-Poon, A. (2009).  Confirmatory adaptive designs with Bayesian decision tools for a targeted therapy in oncology.  Statistics in Medicine 28, 1445–1463.
  * Chapple and Thall (2018) Chapple, A. G. and Thall, P. F. (2018).  Subgroup-specific dose finding in phase I clinical trials based on time to toxicity allowing adaptive subgroup combination.  Pharmaceutical Statistics 17, 734–749.
  * Chen and Lee (2020) Chen, N. and Lee, J. J. (2020).  Bayesian cluster hierarchical model for subgroup borrowing in the design and analysis of basket trials with binary endpoints.  Statistical Methods in Medical Research 29, 2717–2732.
  * Chen et al. (2011) Chen, X., Wang, Z. J., and McKeown, M. J. (2011).  A bayesian lasso via reversible-jump mcmc.  Signal Processing 91, 1920–1932.
  * Chevallier et al. (2023) Chevallier, A., Fearnhead, P., and Sutton, M. (2023).  Reversible jump PDMP samplers for variable selection.  Journal of the American Statistical Association 118, 2915–2927.
  * Chu and Yuan (2018) Chu, Y. and Yuan, Y. (2018).  A Bayesian basket trial design using a calibrated Bayesian hierarchical model.  Clinical Trials 15, 149–158.
  * Freidlin et al. (2010) Freidlin, B., Jiang, W., and Simon, R. (2010).  The cross-validated adaptive signature design.  Clinical Cancer Research 16, 691–698.
  * Freidlin and Simon (2005) Freidlin, B. and Simon, R. (2005).  Adaptive signature design: an adaptive clinical trial design for generating and prospectively testing a gene expression signature for sensitive patients.  Clinical Cancer Research 11, 7872–7878.
  * Gamerman (1997) Gamerman, D. (1997).  Sampling from the posterior distribution in generalized linear mixed models.  Statistics and Computing 7, 57–68.
  * Gu et al. (2016) Gu, X., Chen, N., Wei, C., Liu, S., Papadimitrakopoulou, V. A., Herbst, R. S., and Lee, J. J. (2016).  Bayesian two-stage biomarker-based adaptive design for targeted therapy development.  Statistics in Biosciences 8, 99–128.
  * Jenkins et al. (2011) Jenkins, M., Stone, A., and Jennison, C. (2011).  An adaptive seamless phase II/III design for oncology trials with subpopulation selection using correlated survival endpoints.  Pharmaceutical Statistics 10, 347–356.
  * Jiang et al. (2007) Jiang, W., Freidlin, B., and Simon, R. (2007).  Biomarker-adaptive threshold design: a procedure for evaluating treatment with possible biomarker-defined subset effect.  Journal of the National Cancer Institute 99, 1036–1043.
  * Jin et al. (2020) Jin, J., Riviere, M.-K., Luo, X., and Dong, Y. (2020).  Bayesian methods for the analysis of early-phase oncology basket trials with information borrowing across cancer types.  Statistics in Medicine 39, 3459–3475.
  * Karuri and Simon (2012) Karuri, S. W. and Simon, R. (2012).  A two-stage Bayesian design for co-development of new drugs and companion diagnostics.  Statistics in Medicine 31, 901–914.
  * Liu et al. (2022) Liu, Y., Kairalla, J. A., and Renfro, L. A. (2022).  Bayesian adaptive trial design for a continuous biomarker with possibly nonlinear or nonmonotone prognostic or predictive effects.  Biometrics 78, 1441–1453.
  * Magnusson and Turnbull (2013) Magnusson, B. P. and Turnbull, B. W. (2013).  Group sequential enrichment design incorporating subgroup selection.  Statistics in Medicine 32, 2695–2714.
  * Mehta et al. (2014) Mehta, C., Schäfer, H., Daniel, H., and Irle, S. (2014).  Biomarker driven population enrichment for adaptive oncology trials with time to event endpoints.  Statistics in Medicine 33, 4515–4531.
  * Moore et al. (2020) Moore, C. M., MaWhinney, S., Carlson, N. E., and Kreidler, S. (2020).  A Bayesian natural cubic B-spline varying coefficient method for non-ignorable dropout.  BMC Medical Research Methodology 20, 1–14.
  * Ohwada and Morita (2016) Ohwada, S. and Morita, S. (2016).  Bayesian adaptive patient enrollment restriction to identify a sensitive subpopulation using a continuous biomarker in a randomized phase 2 trial.  Pharmaceutical Statistics 15, 420–429.
  * Park et al. (2020) Park, J. J., Hsu, G., Siden, E. G., Thorlund, K., and Mills, E. J. (2020).  An overview of precision oncology basket and umbrella trials for clinicians.  CA: A Cancer Journal for Clinicians 70, 125–137.
  * Park et al. (2022) Park, Y., Liu, S., Thall, P. F., and Yuan, Y. (2022).  Bayesian group sequential enrichment designs based on adaptive regression of response and survival time on baseline biomarkers.  Biometrics 78, 60–71.
  * Psioda et al. (2021) Psioda, M. A., Xu, J., Jiang, Q., Ke, C., Yang, Z., and Ibrahim, J. G. (2021).  Bayesian adaptive basket trial design using model averaging.  Biostatistics 22, 19–34.
  * Renfro et al. (2014) Renfro, L. A., Coughlin, C. M., Grothey, A. M., and Sargent, D. J. (2014).  Adaptive randomized phase II design for biomarker threshold selection and independent evaluation.  Chinese Clinical Oncology 3, 3489\.
  * Robins (1997) Robins, J. M. (1997).  Causal inference from complex longitudinal data.  In Latent variable modeling and applications to causality, pages 69–117. Springer.
  * Rosenblum et al. (2016) Rosenblum, M., Luber, B., Thompson, R. E., and Hanley, D. (2016).  Group sequential designs with prospectively planned rules for subpopulation enrichment.  Statistics in Medicine 35, 3776–3791.
  * Simon (2017) Simon, R. (2017).  Critical review of umbrella, basket, and platform designs for oncology clinical trials.  Clinical Pharmacology & Therapeutics 102, 934–941.
  * Smolen et al. (2018) Smolen, J. S., Aletaha, D., Anne Barton, G. R. B., Emery, P., Firestein, G. S., Arthur Kavanaugh, I. B. M., Solomon, D. H., Strand, V., and Yamamoto, K. (2018).  Rheumatoid arthritis.  Nature Reviews Disease Primers 4, 18001\.
  * Wang et al. (2007) Wang, S.-J., O’Neill, R. T., and Hung, H. J. (2007).  Approaches to evaluation of treatment effect in randomized clinical trials with genomic subset.  Pharmaceutical Statistics: The Journal of Applied Statistics in the Pharmaceutical Industry 6, 227–244.
  * Wientjes et al. (2022) Wientjes, M. H., den Broeder, A. A., Welsing, P. M., Verhoef, L. M., and van den Bemt, B. J. (2022).  Prediction of response to anti-tnf treatment using laboratory biomarkers in patients with rheumatoid arthritis: a systematic review.  RMD open 8, e002570.