# Mixture of Finite Mixtures Model for Basket Trial

Junxian Geng
Boehringer Ingelheim Pharmaceutical Inc.
and
Guanyu Hu
Department of Statistics, University of Missouri - Columbia Corresponding
author, Email<EMAIL_ADDRESS>
###### Abstract

With the recent paradigm shift from cytotoxic drugs to new generation of
target therapy and immuno-oncology therapy during oncology drug developments,
patients with various cancer (sub)types may be eligible to participate in a
basket trial if they have the same molecular target. Bayesian hierarchical
modeling (BHM) are widely used in basket trial data analysis, where they
adaptively borrow information among different cohorts (subtypes) rather than
fully pool the data together or doing stratified analysis based on each
cohort. Those approaches, however, may have the risk of over shrinkage
estimation because of the invalidated exchangeable assumption. We propose a
two-step procedure to find the balance between pooled and stratified analysis.
In the first step, we treat it as a clustering problem by grouping cohorts
into clusters that share the similar treatment effect. In the second step, we
use shrinkage estimator from BHM to estimate treatment effects for cohorts
within each cluster under exchangeable assumption. For clustering part, we
adapt the mixture of finite mixtures (MFM) approach to have consistent
estimate of the number of clusters. We investigate the performance of our
proposed method in simulation studies and apply this method to Vemurafenib
basket trial data analysis.

Keywords— Bayesian Hierarchical Modeling, Clustering Detection, Nonparametric
Bayesian Method, Oncology

##  1 Introduction

With the recent paradigm shift from cytotoxic drugs to new generation of
target therapy and immuno-oncology therapy during oncology drug developments,
clinical trials in oncology no longer solely target for a specific cancer type
based on the anatomic location of the primary tumor (e.g., breast, lung or
GI). Patients with various cancer (sub)types may be eligible to participate in
a clinical trial if they have the same molecular target. The term basket trial
is used to represent those studies where a single targeted or immuno-onocology
therapy is investigated in the context of multiple diseases or disease
subtypes (Woodcock and LaVange, 2017).

For oncology drug development, we usually start patients studies in Phase I.
The phase I study is typically comprised of two parts, Phase Ia and Phase Ib.
The objective in Phase Ia is to find Maximum Tolerated Dose (MTD) or Optimal
Biological Dose (OBD) based mainly on safety data with the supplements of
initial efficacy information, and we usually call this part dose escalation.
After the optimal dose is determined, we enroll more patients to that dose
level to mainly evaluate efficacy while still monitoring safety of the drug,
where we name as expansion part of the study (Phase Ib). During Phase Ib,
basket trial design are widely used to evaluate the efficacy of the new
treatment in multiple cohorts.

Assume single-arm settings for all cohorts in phase Ib without control arm, a
binary endpoint of response (yes/no) is usually used. There are different
criterion for various tumor types. For example, RECIST criteria (Eisenhauer et
al., 2009) is used for solid tumor and RANO criteria (Lin et al., 2015) is
used for brain tumor. Even with slight different definitions among different
criterion, response is typically defined as a tumor shrinkage more than a
certain threshold (e.g. 30$\%$ in RECIST criteria). The reasonable positive
correlation with time-to-event endpoints such as Progression Free Survival
(PFS) and Overall Survival (OS) along with much earlier readout, makes the
binary endpoint of response widely acceptable efficacy endpoint in early phase
oncology studies. Eventually the treatment effect will be reflected in
response rate, which is defined as the number of patients with responses
divided by the total number of patients. Thus without further discussion on
the validity of the endpoint, a binary endpoint is assumed in this paper.

When dealing with multiple cohorts in one study, there are two directions as
we conduct data analysis: pooled analysis or stratified analysis. If the
assumption is the treatment effect is homogeneous among all cohorts, a pooled
analysis may be conducted where we get a response rate by combining the data
from all cohorts. If the assumption is that different cohorts have totally
heterogeneous treatment effect, separate response rate will be estimated in
each cohort. In practice, however, both of these two assumptions might not be
true. It’s very hard to justify a drug has homogeneous effect on different
tumor types, but it’s also difficult to claim that totally different treatment
effects are expected among those cohorts taking the same drug. In addition,
fully pooled analysis might dilute the positive signal if the drug works only
in some of the cohorts and fully stratified analysis might be either lack of
power or costly inefficient.

In order to mitigate the above mentioned problems, several methods were
proposed to find a balance between pooled analysis and stratified analysis,
where they adaptively borrow information among different cohorts rather than
fully pool the data together. Many of those methods are based on Bayesian
hierarchical modeling (BHM). Thall et al. (2003) and Berry et al. (2013)
proposed exchangeable assumptions among the parameters (logit transformation
of response rate) for different cohorts, which led to shrinkage estimators for
each cohort. Those exchangeable assumptions, however, can easily lead to over
shrinkage to a common mean, thus introduce bias to the estimators.
Neuenschwander et al. (2016) and Chu and Yuan (2018) further allow each
cohort‐specific parameter to be exchangeable with other similar strata
parameters or nonexchangeable with any of them in order to mitigate the risk
of over shrinkage. These methods sometimes may suffer from too many pre-
specified parameters.

In this article, we propose a two-step procedure to find the balance between
pooled and stratified analysis in basket trial design. In the first step, we
treat it as a clustering problem by grouping cohorts into clusters that share
the similar treatment effect. In the second step, we use shrinkage estimator
proposed by Berry et al. (2013) to estimate treatment effects for cohorts
within each cluster under exchangeable assumption. For clustering, we adapt
the mixture of finite mixtures (MFM) approach (Miller and Harrison, 2018; Geng
et al., 2019; Geng and Slate, 2020), which admits a clustering scheme similar
to the famous Chinese restaurant process (CRP) but alleviates the drawback of
CRP by automatic model-based pruning of the tiny extraneous clusters leading
to consistent estimate of the number of clusters. The contribution of this
paper is three-fold. First, a full Bayesian framework is developed and the
clustering results yield useful probabilistic interpretation. In addition, we
establish consistency results for the estimation of clusters for the binary
data. Thirdly, high probability of selecting the correct number of clusters
and the better estimation performance than benchmark methods when
heterogeneity exits between cohorts are empirically demonstrated.

The rest of the paper is organized as follows. We start with a brief review of
the existing popular models for basket trial design and MFM model, as well as
propose our two-step procedures in Section 2. The Bayesian methods for
simultaneous inference on the number of clusters and cluster-specific
estimates are discussed in Section 3 and the MCMC algorithm is provided in
Section 3.1. Simulation studies and comparisons with existing methods are
provided in Section 4 and illustrations of our method by two case studies are
presented in section 5. Finally, we have some discussions in Section 6.

##  2 Methodology

###  2.1 Existing models for basket design

Suppose there are $N$ cohorts in the proposed basket design. In cohort $i$,
there are $n_{i}$ patients enrolled with response rate $p_{i}$. If the number
of responses in cohort $i$ is denoted as $r_{i}$, we have:

$\displaystyle r_{i}\mid
n_{i},p_{i}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\mbox{Binomial}(n_{i},p_{i}),\quad
i=1,\ldots,N$ (1)

In a single arm setting in early phase oncology trials, the goal is to get an
estimation of $p_{i}$ for each cohort $i$ and compare that estimate to the
benchmark response rate threshold. If the estimated response rate is greater
than the pre-specified threshold in a cohort, the sponsor will consider
continuing the development of the tested drug for that indication.

A popular Bayesian hierarchical modeling used to get estimations of $p_{i}$’s
is proposed by Berry et al. (2013), where they use shrinkage estimators for
each cohort under exchangeable assumptions and we call this approach Berry’s
BHM in the rest of the paper. Then model and prior can be briefly expressed
as:

$\displaystyle\theta_{i}=\text{log}(\frac{p_{i}}{1-p_{i}})-\text{log}(\frac{p_{Ti}}{1-p_{Ti}})$
(2) $\displaystyle\theta_{i}\sim\mbox{N}(\mu,\sigma^{2}),\quad i=1,\ldots,N,$
$\displaystyle\mu\sim\mbox{N}(0,2^{2}),\quad\sigma^{2}\sim\mbox{HN}(0,1),$
$\displaystyle r_{i}\mid
n_{i},p_{i}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\mbox{Binomial}(n_{i},p_{i}),\quad
i=1,\ldots,N$

Where $p_{Ti}$ is usually set as the benchmark response rate or control rate
for cohort $i$. $\theta_{i}$ can be understood as the logit transformed
treatment effect relative to the control rate. The values of the
hyperparameters in aforementioned model for $\mu$ and $\sigma^{2}$ are popular
used and they are generally not too sensitive to different choices.

The performance of this model, however, is highly dependent on the validity of
the exchangeable assumptions as well as the choices of control rate, which
sometimes is hard to justify in real studies. If the exchangeable assumptions
are not valid for $\theta$’s, there is high risk of over shrinkage on the
estimations of cohort-specific response rate.

Neuenschwander et al. (2016) proposed a method that further allow each
cohort‐specific parameter to be exchangeable with other similar strata
parameters or nonexchangeable with any of them, which is a data-drive approach
to mitigate the risk of over shrinkage and we call this approach EXNEX in the
rest of the paper. The model and prior can be briefly expressed as:

$\displaystyle\theta_{i}=\text{log}(\frac{p_{i}}{1-p_{i}})$ (3)
$\displaystyle\mbox{EX: with probability }\pi_{i};\quad$
$\displaystyle\theta_{i}\sim\mbox{N}(\mu,\sigma^{2}),\quad i=1,\ldots,N,$
$\displaystyle\mu\sim\mbox{N}(a,b^{2}),\quad\sigma^{2}\sim\mbox{HN}(c,d),$
$\displaystyle\mbox{NEX: with probability }1-\pi_{i};\quad$
$\displaystyle\theta_{i}\sim\mbox{N}(\mu_{i},\sigma^{2}_{i}),\quad
i=1,\ldots,N,$ $\displaystyle r_{i}\mid
n_{i},p_{i}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\mbox{Binomial}(n_{i},p_{i}),\quad
i=1,\ldots,N$

Where $\pi_{i}$ is the weight of exchangeable component for cohort $i$, which
is usually set the same for all cohorts by default; the hyperparameters
$a,b,c,d$ for the exchangeable part is usually from the method for setting
weakly informative priors given in the on-line Appendix of Neuenschwander et
al. (2016) unless there are specific prior information; the hyperparameters
$\mu_{i}$ and $\sigma^{2}_{i}$ for each cohort are usually set so that
$\mu_{i}$ equals to the log-odds of a plausible guess (or benchmark response
rate) for cohort $i$ and $\sigma^{2}_{i}$ makes approximately one observation.

This model sometimes may suffer from too many pre-specified parameters or
parameters that need to be tuned.

###  2.2 Mixture of Finite Mixtures model for binary data

To mitigate the possible issues mentioned above, we first cluster the cohorts
into groups based on the observed response rates and then get shrinkage
estimators within each group by assuming the exchangeable assumption.

Bayesian models offer a natural solution to simultaneously estimate the number
of clusters and cluster assignments. The Chinese restaurant process (CRP;
Neal, 2000) offers choices to allow for uncertainty in the number of clusters
by assigning a prior distribution on the cluster assignments parameters
$(z_{1},z_{2},\ldots,z_{N})$, where $z_{i}$ denotes the cluster assignment for
cohort $i$. In the CRP, $z_{i},i=2,\ldots,N$ have the following conditional
distribution (i.e., a Pólya urn scheme, Blackwell et al., 1973)

$\displaystyle P(z_{i}\mid
z_{1},\ldots,z_{i-1})\propto\begin{cases}\absolutevalue{c},\quad\text{at an
existing cluster labeled}\,c\\\ \alpha,\quad\quad\quad\,\text{at a new
cluster}.\end{cases}$ (4)

Here $\absolutevalue{c}$ refers to the size of cluster labeled $c$, and
$\alpha$ is the concentration parameter of the underlying Dirichlet process.
The prior distribution on the cluster assignments induces a prior distribution
on the sizes of the clusters in the partition. Let $\mathcal{C}_{N}$ denotes a
partition of the set $\\{1,2,3,\ldots,N\\}$ and
$t=\absolutevalue{\mathcal{C}_{N}}$ denote the number of blocks in the
partition $\mathcal{C}_{N}$.

Under (4), one can obtain the probability of block-sizes
$\bm{b}=(b_{1},b_{2},\ldots,b_{t})$ of a partition $\mathcal{C}_{N}$ as

$\displaystyle p_{\mathrm{DP}}(\bm{b})\propto\prod_{j=1}^{t}b_{j}^{-1}.$ (5)

It is clear from (5) that CRP tends to assign large probabilities to highly
imbalanced cluster sizes in which, necessarily, some clusters will be quite
small. This results in producing extraneous clusters in the posterior that
leads inconsistent estimation on the number of clusters even when the sample
size goes to infinity.

Miller and Harrison (2018) proposed a modification to the CRP, which is called
a mixture of finite mixtures (MFM) model, to mitigate this issue:

$\displaystyle\begin{split}k&\sim p(\cdot),\\\
\pi=(\pi_{1},\ldots,\pi_{k})\mid k&\sim\mbox{Dir}(\gamma,\ldots,\gamma),\\\
z_{i}\mid k,\pi&\sim\sum_{s=1}^{k}\pi_{s}\delta_{s},\quad
i=1,\ldots,N,\end{split}$ (6)

where $p(\cdot)$ is a proper probability mass function on $\\{1,2,\ldots\\}$,
and $\delta_{s}$ is a point-mass at $s$. The joint distribution of
$(z_{1},\ldots,z_{n})$ under (6) admits a Pólya urn scheme akin to the CRP:

$\displaystyle P(z_{i}\mid
z_{1},\ldots,z_{i-1})\propto\begin{cases}\absolutevalue{c}+\gamma,\quad\text{at
an existing cluster labeled}\,c\\\
\frac{V_{i}(t+1)}{V_{i}(t)}\gamma,\quad\quad\quad\,\text{at a new
cluster},\end{cases}$ (7)

where the coefficients $V_{i}(t)$ are given by,

$\displaystyle\begin{split}V_{i}(t)&=\sum_{k=1}^{+\infty}\dfrac{k_{(t)}}{(\gamma
k)^{(i)}}p(k)\end{split}$

where $k_{(t)}=k(k-1)...(k-t+1)$, and $(\gamma k)^{(i)}={\gamma k}(\gamma
k+1)...(\gamma k+i-1)$. (By convention, $x^{(0)}=1$ and $x_{(0)}=1$). While
this restaurant process bears close resemblance to the CRP, the introduction
of new clusters is slowed down by a factor
$V_{i}(\absolutevalue{\mathcal{C}_{i-1}}+1)/V_{n}(\absolutevalue{\mathcal{C}_{i-1}})$,
thereby pruning the tiny extraneous clusters.

An alternative way to understand the natural pruning of extraneous clusters is
through the probability distribution induced on the cluster sizes. Again, let
$\mathcal{C}_{N}$ denotes a partition with block-sizes
$\bm{b}=(b_{1},b_{2},\ldots,b_{t})$ and $t=|\mathcal{C}_{N}|$ under MFM. In
contrast to (5), the probability of the cluster sizes $(b_{1},\ldots,b_{t})$
under the MFM is

$\displaystyle
p_{\mathrm{MFM}}(\bm{b})\propto\prod_{j=1}^{t}b_{j}^{\gamma-1}.$ (8)

From (5) and (8), the comparison easily reveals that MFM assigns comparatively
smaller probabilities to highly imbalanced cluster sizes.

We adapt MFM to our model setting in order to group different cohorts into
clusters, then the model and prior can be expressed hierarchically as:

$\displaystyle k\sim q(\cdot),\text{where $q(\cdot)$ is a p.m.f on \\{1,2,
\ldots\\} }$ $\displaystyle
P_{s}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\mbox{Beta}(\alpha,\beta),\quad
s=1,\ldots,k,$ $\displaystyle\mbox{pr}(z_{i}=s\mid\bm{\pi},k)=\pi_{s},\quad
s=1,\ldots,k,\,i=1,\ldots,N,$ $\displaystyle\bm{\pi}\mid
k\sim\mbox{Dirichlet}(\gamma,\ldots,\gamma),$ $\displaystyle r_{i}\mid
n_{i},\bm{z},\bm{P},k\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\mbox{Binomial}(n_{i},p_{i}),\quad
p_{i}=P_{z_{i}},\,i=1,\ldots,N,$ (9)

We assume $q(\cdot)$ is a $\mbox{Poisson}(1)$ distribution truncated to be
positive through the rest of the paper, which has been proved by Miller and
Harrison (2018) and Geng et al. (2019) to guarantee consistency for the mixing
distribution and the number of clusters. Here, we define
$G=\sum_{s=1}^{k}\pi_{s}\delta_{P_{s}}$, where $\delta$ is the point mass
measure, and $P_{s}$ is the collection of parameters for the binomial
distribution in cluster $s$ for $s=1,\ldots,k$ and $P$ is Beta probability
measure.

Let $k_{0}$, $G_{0}$, $P_{0}$ be the true number of clusters, the true mixing
measure, and the corresponding probability measure, respectively. Then the
following proposition establishes the posterior consistency and contraction
rate for $k$ and $G$.

######  Proposition 1.

Let $\Pi_{N}(\cdot\mid r_{1},\ldots,r_{N})$ be the posterior distribution
obtained from (2.2) given a random sample $r_{1},\ldots,r_{N}$. Assume that
all the parameters are restricted to a compact space $\Theta^{*}$. Then we
have

$\displaystyle\Pi_{N}(k=k_{0}\mid r_{1},\ldots,r_{N})\rightarrow
1,~{}~{}\Pi_{N}(W(G,G_{0})\lesssim(\log N/N)^{-1/2}\mid
r_{1},\ldots,r_{N})\rightarrow 1,$

almost surely under $P_{0}$ as $N\rightarrow\infty$.

In order to prove the Proposition 1, we need to verify the conditions
(P.1)-(P.4) in Guha et al. (2019) hold. Condition (P.1) is satisfied since we
restrict our parameters of interest to a compact space $\bm{\Theta^{*}}$ and
beta distribution is first-order identifiable. Condition (P.2) also holds
since we assign a non-zero continuous distribution on the parameters within a
bounded support. Beta base distribution is sufficient for Condition (P.3).
Condition (P.4) holds since we choose a truncated Poisson distribution on
$q(\cdot)$. The proof can be finished by using the results in Yin et al.
(2020).

###  2.3 Proposed two-step procedure

We now propose a two-step procedure to find the balance between pooled and
stratified analysis in basket trial design. In the first step, we treat it as
a clustering problem by grouping cohorts into clusters that share the similar
treatment effect using the model proposed in (2.2). In the second step, we
further address heterogeneity for cohorts within each cluster by fitting the
Bayesian Hierarchical model by Berry et al. (2013) in each cluster. For
cluster $s$, the model and prior can be expressed as:

$\displaystyle\theta_{j}=\text{log}(\frac{p_{j}}{1-p_{j}})-\text{log}(\frac{p_{T}}{1-p_{T}})$
(10) $\displaystyle\theta_{j}\sim\mbox{N}(\mu,\sigma^{2}),\quad
j=1,\ldots,n_{s},$
$\displaystyle\mu\sim\mbox{N}(0,2^{2}),\quad\sigma^{2}\sim\mbox{HN}(0,1),$
$\displaystyle r_{j}\mid
n_{j},p_{j}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\mbox{Binomial}(n_{j},p_{j}),\quad
j=1,\ldots,n_{s},$

where $n_{s}$ denotes the number of cohorts in cluster $s$.

##  3 Bayesian Inference

###  3.1 MCMC Algorithm

For the step one of our proposed methods, we design an efficient Gibbs sampler
algorithm without reversible jumps algorithm which is easier to get the
posterior samples for clustering labels. The Gibbs sampler for step one is
presented in Algorithm 1.

Algorithm 1 Collapsed sampler for MFM-BD

1:procedure c-MFM-BD

2:Initialize $\bm{z}=(z_{1},\ldots,z_{N})$ and $\bm{P}=(P_{s})$.

3: for each iter $=1$ to M do

4:Update $\bm{P}=(P_{s})$ conditional on $z$ in a closed form as
$\displaystyle\begin{split}P_{s}\mid\bm{r},\bm{n},z&\sim\mbox{Beta}(\alpha+r_{[s]},\beta+N_{s}-r_{[s]})\end{split}$
Where ${r}_{[s]}=\sum_{z_{i}=s}r_{i}$ and ${N}_{s}=\sum_{z_{i}=s}n_{i}$,
$r=1,\ldots,k$. Here $k$ is the number of clusters formed by current $z$.

5:Update $\bm{z}=(z_{1},\ldots,z_{N})$ conditional on $\bm{P}=(P_{s})$, for
each $i$ in $(1,\ldots,N)$, we can get a closed form expression for
$P(z_{i}=c\mid z_{-i},\bm{r},\bm{n},\bm{P})$:
$\propto\left\\{\begin{array}[]{ll}[\absolutevalue{c}+\gamma]\mbox{
}\text{dbinom}(r_{i},n_{i},P_{z_{i}})&\text{at an existing table $c$}\\\
\frac{V_{n}(\absolutevalue{\mathcal{C}_{-i}}+1)}{V_{n}(\absolutevalue{\mathcal{C}_{-i}}}\gamma
m(\mathscr{S}_{i})&\text{if $c$ is a new table}\\\ \end{array}\right.$ where
$\mathcal{C}_{-i}$ denotes the partition obtained by removing $z_{i}$ and
$\displaystyle
m(\mathscr{S}_{i})=\binom{n_{i}}{r_{i}}\frac{B(r_{i}+\alpha,n_{i}-r_{i}+\beta)}{B(\alpha,\beta)}$
Where $B(.)$ is beta function.

6: end for

7:end procedure

For step two, we don’t derive MCMC algorithm ourselves since the full
conditional distributions are not in closed form. The gMAP function in Rbest
package in R is used to handle Berry’s BHM.

###  3.2 Post MCMC inference

Another important task for our proposed method is the inference of MCMC
results.

We first discussion the inference on cluster assignment parameter $\bm{z}$.
The posterior mean or median of $\bm{z}$ is not suitable. Dahl’s method (Dahl,
2006) provides a remedy for posterior inference of the clustering
configurations $\bm{z}$. It chooses the iteration in the posterior sample that
optimizes a least squares criterion as the estimate for $\bm{z}$. We firstly
define an $n\times n$ membership matrix

$\displaystyle B=(B(i,j))=\big{(}1(z_{i}=z_{j})\big{)},\quad
i,j\in\\{1,2,\dots,n\\},$ (11)

where $1(\cdot)$ is the identical function. The $(i,j)$-th element of
membership matrix is 1 if the $i$-th and $j$-th observations belong to the
same component (or have the same baseline intensity); it is 0 otherwise. For
the MCMC samples, we calculate the $B^{(i)}$ for $i$-th iterations. And then
calculate the element-wise mean of the membership matrices
$\bar{B}=\frac{1}{L}\sum_{l=1}^{L}{B}^{(l)}$, where the summation is element-
wise. Finally, we identify the most _close_ posterior draw as the one that is
closest to $\bar{B}$ with respect to the element-wise Euclidean distance
$\sum_{i=1}^{n}\sum_{j=1}^{n}(B^{(l)}(i,j)-\bar{B}(i,j))^{2}$ among the
retained $l=1,\ldots,L$ posterior draws. The posterior estimates of cluster
memberships $z_{1},\ldots,z_{n}$ can be summarized based on the draw
identified by Dahl’s method.

Within each cluster, the cohort specific parameter $\theta_{j}$’s will be
summarized by posterior mean and the Bayesian estimate for response rate
$p_{j}$ will be directly from solving the equation in (10).

##  4 Simulation

In this section, we investigate the performance of our proposed method from a
variety of measures.

###  4.1 Simulation Settings and Evaluation Metrics

In the simulation study, we have five different scenarios with two different
sample sizes (20 or 30) in each cohort. The response rates for each scenarios
are shown in Table 1.

Scenario | Response Rate Design
---|---
Scenario 1 | $P_{1}=P_{2}=\ldots=P_{10}=0.4$
Scenario 2 |  $P_{1}=\ldots=P_{5}=0.2$, $P_{6}=\ldots=P_{10}=0.6$
Scenario 3 |  $P_{1}=\ldots=P_{5}=0.2$, $P_{6}=\ldots=P_{10}=0.5$
Scenario 4 |  $P_{1}=P_{2}=P_{3}=0.1$, $P_{4}=P_{5}=P_{6}=0.4$, $P_{7}=\ldots=P_{10}=0.7$
Scenario 5 | $P_{1}=P_{2}=\ldots=P_{10}=0.2$
Table 1: Simulation designs for different scenarios

In all the simulation scenarios considered below, we employed Algorithm 1 with
$\gamma$ = 1 and $\alpha=\beta=1$ to fit the MFM-BD model. We arbitrarily
initialized our algorithm with 5 clusters and randomly allocated the cluster
assignments. We experimented with various other choices and did not find any
evidence of sensitivity to the initialization. Results from our method is
based on 5000 MCMC iterations leaving out a burn-in of 2000 in step one and
8000 MCMC iterations leaving out a burn-in of 2000 in step two.

For Berry’s BHM and EXNEX approach, the choice of hyper-parameters are from
the default settings as stated in section 2.1. Results from those two
approaches are based on 10000 MCMC iterations leaving out a burn-in of 2000.

For our proposed method, the estimated number of clusters $\hat{K}$ for each
replicate is summarized from the MCMC iteration picked by Dahl’s method which
is introduced in Section 3.2. The performance of the posterior estimates of
parameters for all methods were evaluated by the average absolute bias (AAB)
and the mean square error (AMSE) in the following ways, take $\theta_{i}$ as
an example:

$\begin{split}\text{AAB}=\frac{1}{10}\sum_{i=1}^{10}\left|\frac{1}{500}\sum_{r=1}^{500}(\theta^{r}_{i}-\theta^{0}_{i})\right|,\\\
\text{MMSE}=\frac{1}{10}\sum_{i=1}^{10}\sqrt{\frac{1}{500}\sum_{r=1}^{500}(\theta^{r}_{i}-\theta^{0}_{i})^{2}},\end{split}$

where $\theta^{r}_{i}$ is the posterior mean of rth replicate for ith cohort
and $\theta^{0}_{i}$ is the true value of ith cohort.

###  4.2 Simulation Results

First, we present the cluster estimates of our proposed methods in Table 2.
From the results shown in Table 2, we see that our proposed methods
successfully recover the number of clusters within a reasonable range for all
five different scenarios under two different sample size. In addition, the
large sample size will cause better estimates of the number of clusters by
comparing results between $n=20$ and $N=30$.

Table 2: Clustering Performance for Simulation Studies Scenario | n | $\hat{K}$ |  S.D. of $\hat{K}$
---|---|---|---
Scenario 1 | $n=20$ | 1.046 | 0.210
| $n=30$ | 1.036 | 0.186
Scenario 2 | $n=20$ | 2.168 | 0.417
| $n=30$ | 2.131 | 0.346
Scenario 3 | $n=20$ | 1.945 | 0.516
| $n=30$ | 2.105 | 0.376
Scenario 4 | $n=20$ | 2.661 | 0.612
| $n=30$ | 2.929 | 0.559
Scenario 5 | $n=20$ | 1.026 | 0.159
| $n=30$ | 1.020 | 0.140

Furthermore, we compare our proposed methods with Berry’s BHM and EXNEX under
the criteria we proposed in Section 4.1. The simulation results are presented
in Table 3 and Table 4. All simulations results are based on 500 replicates.
From the results shown in Table 3 and Table 4, we see that our proposed
methods have better estimation performance based on AAB when heterogeneity
exits among different cohorts and the AMSE is comparable with benchmark
methods.

Table 3: Simulation Results ($n=20$)  |  | Berry’s BHM | EXNEX | MFM-BD
---|---|---|---|---
Scenario 1 | AAB | 0.0015 | 0.0026 | 0.0025
| AMSE | 0.0024 | 0.0056 | 0.0026
Scenario 2 | AAB | 0.0430 | 0.0274 | 0.0131
| AMSE | 0.0091 | 0.0095 | 0.0069
Scenario 3 | AAB | 0.0485 | 0.0304 | 0.0266
| AMSE | 0.0087 | 0.0093 | 0.0096
Scenario 4 | AAB | 0.0305 | 0.0256 | 0.0261
| AMSE | 0.0086 | 0.0086 | 0.0102
Scenario 5 | AAB | 0.0023 | 0.0064 | 0.0119
| AMSE | 0.0016 | 0.037 | 0.0016
Table 4: Simulation Results ($n=30$)  |  | Berry’s BHM | EXNEX | MFM-BD
---|---|---|---|---
Scenario 1 | AAB | 0.0012 | 0.0019 | 0.0015
| AMSE | 0.0016 | 0.0036 | 0.0017
Scenario 2 | AAB | 0.0303 | 0.0180 | 0.0056
| AMSE | 0.0063 | 0.0064 | 0.0035
Scenario 3 | AAB | 0.0349 | 0.0223 | 0.0107
| AMSE | 0.0061 | 0.0065 | 0.0053
Scenario 4 | AAB | 0.0207 | 0.0179 | 0.0136
| AMSE | 0.0059 | 0.0060 | 0.0061
Scenario 5 | AAB | 0.0015 | 0.0039 | 0.0068
| AMSE | 0.0011 | 0.0024 | 0.0010

In order to have a closer look of estimation performance, we have the boxplots
of posterior estimates of different cohorts in Figure 1.

(a) $n=20$

(b) $n=30$

Figure 1: Posterior Estimates Boxplots for Different Cohorts under Scenario 3

From the results shown in Figure 1, we see that in each cohort our proposed
method performs better than other two methods when there exit heterogeneities
among different cohorts.

##  5 Real Data Analysis

In this section, we apply the proposed MFM-BD to the analysis of a famous
basket trial (Hyman et al., 2015). This trial is a histology-independent
“basket” study of vemurafenib in BRAF V600 mutation–positive nonmelanoma
cancers, where patients were enrolled in six prespecified cancer cohorts
including anaplastic thyroid cancer (ATC), Erdheim-Chester disease or
Langerhans’-cell histiocytosis (ECD/LCH), cholangiocarcinoma (CCA), colorectal
cancer treated by vemurafenib (CRC-V), colorectal cancer treated by
vemurafenib and cetuximab (CRCVC), and non small cell lung cancer(NSCLC). The
primary end point was the response rate. A total of 84 evaluable patients are
included in the analysis. The sample size and observed responses for each
cohort are presented in Table 5.

Results from our method (MFM-BD) is based on 5000 MCMC iterations leaving out
a burn-in of 2000 in step one and 8000 MCMC iterations leaving out a burn-in
of 2000 in step two; initialized at a randomly generated configuration with 5
clusters. From Table 5, 6 subtypes are grouped into two clusters. The first
cluster includes ATC, ECD/LCH and NSCLC and the the second cluster includes
CCA, CRC-V and CRC-VC.

For comparison, we also conduct analysis based on Berry’s BHM and EXNEX
approach. The parameter settings for both methods are consistent from those
used in simulation studies. From the results in Table 5, we can see that when
comparing to the other approaches, our proposed MFM-BD give the estimations
more shrinkage within each cluster and more divergence between clusters.

Table 5: Response estimations Results Subtype | $r/n$ | $\%$ | Posterior mean (CI)
---|---|---|---
MFM-BD | Berry’s BHM | EXNEX
Cluster 1 | 1\. ATC | 2/7 | 28.6 | 31.8 (6.9,62.3) | 26.7 (6.3,57.6) | 27.5 (5.1,60.1)
2\. ECD/LCH | 6/14 | 42.9 | 41.6 (20.3,65.7) | 38.3 (17.7,63.2) | 40.7 (18.1,65.3)
6\. NSCLC | 8/19 | 42.1 | 41.2 (21.9,61.9) | 38.6 (19.8,60.5) | 40.5 (20.6,62.3)
Cluster 2 | 3\. CCA | 1/8 | 12.5 | 12.6 (1.7,35.7) | 16.2 (2.1,39.9) | 15.4 (1.5,42.5)
4\. CRC-V | 1/26 | 3.8 | 6.6 (0.7,16.7) | 7.5 (0.9,20.3) | 5.9 (0.6,16.7)
5\. CRC-VC | 0/10 | 0.0 | 6.5 (0.1,20.4) | 8.3 (0.2,26.5) | 5.8 (0.1,23.3)

##  6 Discussion

In this article, we propose a two-step procedure to evaluate the efficacy of
the new treatment (binary response endpoint) in multiple cohorts in a basket
trial, which has excellent performance in both simulation and real data
examples. The article provides a new way of borrowing information across
different cohorts by firstly clustering into groups without pre-specifying the
number of groups, and then borrow information within each group under
exchangeable assumption. The clustering of those cohorts can give further
directions for development of a new drug. The approach is also proved to yield
consistent detection of the number of clusters.

While we recommend using Berry’s BHM in the second step of our approach, other
popular approaches such as EXNEX model can also be implemented. This proposed
approach is more appropriate for basket trials where the number of cohorts is
relatively large (e.g. 6 cohorts or more).

Further extension of this work includes: analysing time-to-event endpoint such
as Progression Free Survival (PFS) and Overall Survival (OS) (Xu et al., 2019)
under the same framework; adding another layer of heterogeneity so that it can
also be used in platform trials where there are multiple treatment options and
multiple cancer types.

## References

  * Berry et al. (2013) Berry, S. M., K. R. Broglio, S. Groshen, and D. A. Berry (2013).  Bayesian hierarchical modeling of patient subpopulations: efficient designs of phase ii oncology clinical trials.  Clinical Trials 10(5), 720–734.
  * Blackwell et al. (1973) Blackwell, D., J. B. MacQueen, et al. (1973).  Ferguson distributions via Pólya urn schemes.  The Annals of Statistics 1(2), 353–355.
  * Chu and Yuan (2018) Chu, Y. and Y. Yuan (2018).  A bayesian basket trial design using a calibrated bayesian hierarchical model.  Clinical Trials 15(2), 149–158.
  * Dahl (2006) Dahl, D. B. (2006).  Model-based clustering for expression data via a dirichlet process mixture model.  Bayesian inference for gene expression and proteomics 4, 201–218.
  * Eisenhauer et al. (2009) Eisenhauer, E. A., P. Therasse, J. Bogaerts, L. H. Schwartz, D. Sargent, R. Ford, J. Dancey, S. Arbuck, S. Gwyther, M. Mooney, et al. (2009).  New response evaluation criteria in solid tumours: revised recist guideline (version 1.1).  European journal of cancer 45(2), 228–247.
  * Geng et al. (2019) Geng, J., A. Bhattacharya, and D. Pati (2019).  Probabilistic community detection with unknown number of communities.  Journal of the American Statistical Association 114(526), 893–905.
  * Geng and Slate (2020) Geng, J. and E. H. Slate (2020).  Discovery among binary biomarkers in heterogeneous populations.  In Statistical Modeling in Biomedical Research, pp. 213–232. Springer.
  * Guha et al. (2019) Guha, A., N. Ho, and X. Nguyen (2019).  On posterior contraction of parameters and interpretability in bayesian mixture modeling.  arXiv preprint arXiv:1901.05078.
  * Hyman et al. (2015) Hyman, D. M., I. Puzanov, V. Subbiah, J. E. Faris, I. Chau, J.-Y. Blay, J. Wolf, N. S. Raje, E. L. Diamond, A. Hollebecque, et al. (2015).  Vemurafenib in multiple nonmelanoma cancers with braf v600 mutations.  New England Journal of Medicine 373(8), 726–736.
  * Lin et al. (2015) Lin, N. U., E. Q. Lee, H. Aoyama, I. J. Barani, D. P. Barboriak, B. G. Baumert, M. Bendszus, P. D. Brown, D. R. Camidge, S. M. Chang, et al. (2015).  Response assessment criteria for brain metastases: proposal from the rano group.  The lancet oncology 16(6), e270–e278.
  * Miller and Harrison (2018) Miller, J. W. and M. T. Harrison (2018).  Mixture models with a prior on the number of components.  Journal of the American Statistical Association 113(521), 340–356.
  * Neal (2000) Neal, R. M. (2000).  Markov chain sampling methods for Dirichlet process mixture models.  Journal of Computational and Graphical Statistics 9(2), 249–265.
  * Neuenschwander et al. (2016) Neuenschwander, B., S. Wandel, S. Roychoudhury, and S. Bailey (2016).  Robust exchangeability designs for early phase clinical trials with multiple strata.  Pharmaceutical statistics 15(2), 123–134.
  * Thall et al. (2003) Thall, P. F., J. K. Wathen, B. N. Bekele, R. E. Champlin, L. H. Baker, and R. S. Benjamin (2003).  Hierarchical bayesian approaches to phase ii trials in diseases with multiple subtypes.  Statistics in medicine 22(5), 763–780.
  * Woodcock and LaVange (2017) Woodcock, J. and L. M. LaVange (2017).  Master protocols to study multiple therapies, multiple diseases, or both.  New England Journal of Medicine 377(1), 62–70.
  * Xu et al. (2019) Xu, Y., P. Müller, A. M. Tsimberidou, and D. Berry (2019).  A nonparametric Bayesian basket trial design.  Biometrical Journal 61(5), 1160–1174.
  * Yin et al. (2020) Yin, F., G. Hu, and W. Shen (2020).  Analysis of professional basketball field goal attempts via a Bayesian matrix clustering approach.  arXiv preprint arXiv:2010.08495.