# Bayesian Model Averaging for Partial Ordering Continual Reassessment Methods

Luka Kovačević
MRC Biostatistics Unit, University of Cambridge
Thomas Jaki
Department of Machine Learning and Data Science, University of Regensburg
MRC Biostatistics Unit, University of Cambridge
Helen Barnett
Lancaster University
and
Pavel Mozgunov
MRC Biostatistics Unit, University of Cambridge

###### Abstract

Phase I clinical trials are essential to bringing novel therapies from
chemical development to widespread use. Traditional approaches to dose-finding
in Phase I trials, such as the ’3+3’ method and the Continual Reassessment
Method (CRM), provide a principled approach for escalating across dose levels.
However, these methods lack the ability to incorporate uncertainty regarding
the dose-toxicity ordering as found in combination drug trials. Under this
setting, dose-levels vary across multiple drugs simultaneously, leading to
multiple possible dose-toxicity orderings. The Partial Ordering CRM (POCRM)
extends to these settings by allowing for multiple dose-toxicity orderings. In
this work, it is shown that the POCRM is vulnerable to ’estimation
incoherency’ whereby toxicity estimates shift in an illogical way, threatening
patient safety and undermining clinician trust in dose-finding models. To this
end, the Bayesian model averaged POCRM (BMA-POCRM) is proposed. BMA-POCRM uses
Bayesian model averaging to take into account all possible orderings
simultaneously, reducing the frequency of estimation incoherencies. The
effectiveness of BMA-POCRM in drug combination settings is demonstrated
through a specific instance of estimate incoherency of POCRM and simulation
studies. The results highlight the improved safety, accuracy and reduced
occurrence of estimate incoherency in trials applying the BMA-POCRM relative
to the POCRM model.

Keywords: Bayesian inference, Adaptive design, Maximum Tolerable Dose,
Posterior probability, Model averaging

##  1 Introduction

An aim of in-patient Phase I clinical trials is to determine the maximum
tolerable dose (MTD) from a range of dose levels for progression into Phase II
trials. The MTD is the dose level that matches the target probability of
unwanted effects, sometimes called the target toxicity rate (TTR). Depending
on the setting, the MTD can be either a dose level for a single agent or a
combination of dose levels in a multi-agent setting. A common approach to
determining the MTD is to use the number of dose-limiting toxicities (DLTs),
simply called toxicities, on each dose to compute a probability of toxicity.
The definition of a DLT and the desired TTR are set in advance of a trial,
hence, by estimating the risk of DLT for each dose, the dose level with a
toxicity rate closest to the TTR can be selected. Dose-escalation procedures
aim to estimate the location of the true MTD given a TTR and sequential
patient DLT data.

Approaches to dose-escalation in Phase I clinical trials frequently rely on
the assumption of simple dose orderings. In the case of single drug Phase I
trials, a simple ordering can be constructed by assuming a monotonic
relationship between dose level and toxicity. That is, higher doses of a drug
are assumed to be more toxic than lower doses. For example, suppose there are
dose levels $d_{1},\ldots,d_{4}$ where a higher index indicates a higher dose
level. This means that $d_{1}<\ldots<d_{4}$ implies the following simple
ordering assuming dose-toxicity monotonicity,

$d_{1}\rightarrow d_{2}\rightarrow d_{3}\rightarrow d_{4}.$

This also implies that $d_{2}$ is more toxic than $d_{1}$, $d_{3}$ is more
toxic than both $d_{1}$ and $d_{2}$ and so forth. The continual reassessment
method (CRM) is one such approach that is based on the assumption of
monotonicity. It utilises a Bayesian framework to update estimates of the risk
of toxicity to guide dose-escalation based on a given toxicity ordering of
doses (O’Quigley et al., 1990). Several independent studies have shown that
escalation based on the CRM leads to favourable operating characteristics for
finding the true MTD in single agent clinical trials. With the growing need
for combination drug trials, where the toxicity profile of joint
administration of two or more drugs is investigated, methods that allow for
potential uncertainties in the ordering of dose levels are necessary (Mozgunov
et al., 2020).

|  | Drug B
---|---|---
|  | $1$ | $2$
Drug A | $1$ | $\bm{d_{1}}$ | $\bm{d_{2}}$
$2$ | $\bm{d_{3}}$ | $\bm{d_{4}}$
$3$ | $\bm{d_{5}}$ | $\bm{d_{6}}$
Table 1: $3\times 2$ setting resulting in 6 dose levels.

Escalation of multiple drugs concurrently creates uncertainty in the dose-
escalation process as the change in dose toxicity is not obvious for diagonal
transitions in dose level, for example where one drug increases in dose level
and the other decreases. Specifically, it may be reasonable to assume
monotonicity for a single drug, however, this does not extend to multiple
changes in dose level. In Table 1, the 3-by-2 dose level configuration for a
dual drug combination trial with 3 dose levels for drug $A$ and 2 dose levels
for drug $B$ is shown. Here, it is unknown prior to the trial whether toxicity
increases or decreases with a shift between $d_{2}$ and $d_{3}$ and likewise
for $d_{4}$ and $d_{5}$. Based on this dose level matrix, the 5 simple
orderings comprising the partial ordering of dose levels are,

$\displaystyle 1:\;d_{1}\rightarrow d_{2}\rightarrow d_{3}\rightarrow
d_{4}\rightarrow d_{5}\rightarrow d_{6},$ $\displaystyle 2:\;d_{1}\rightarrow
d_{3}\rightarrow d_{5}\rightarrow d_{2}\rightarrow d_{4}\rightarrow d_{6},$
$\displaystyle 3:\;d_{1}\rightarrow d_{3}\rightarrow d_{2}\rightarrow
d_{5}\rightarrow d_{4}\rightarrow d_{6},$ $\displaystyle 4:\;d_{1}\rightarrow
d_{2}\rightarrow d_{3}\rightarrow d_{5}\rightarrow d_{4}\rightarrow d_{6},$
$\displaystyle 5:\;d_{1}\rightarrow d_{3}\rightarrow d_{2}\rightarrow
d_{4}\rightarrow d_{5}\rightarrow d_{6}.$

Several procedures have been developed to handle the problem of uncertain
dose-toxicity orderings, which is a problem that also persists in dose-
schedule (Mozgunov and Jaki, 2019) and combination-schedule settings (Mozgunov
et al., 2022; Riviere et al., 2015). Ivanova and Wang (2004) developed the
_up-and-down design for combinations_ using an isotonic regression combined
with the Narayana design to escalate doses algorithmically (Ivanova et al.,
2003). A further development updated this method to utilise a _T-statistic_
(Ivanova and Kim, 2009) for escalation decisions. Further, Yin and Yuan (2009)
have developed the Bayesian copula regression-based model, which uses a
Bayesian scheme similar to the CRM to update posterior estimates of toxicity
and recommend dose allocation decisions. Most recently, Mozgunov et al. (2020)
proposed a beta distributed surface-free approach to handling drug combination
trials, successful in reducing the average number of toxicities within a
trial. The Continual Reassessment Method for Partial Orderings (POCRM),
developed by Wages et al. (2011), extends the Bayesian framework of the CRM to
allow for several dose-toxicity orderings to be specified.

Further, the POCRM addresses the problem of uncertainty in dose-toxicity
orderings by selecting from the set of proposed dose-toxicity orderings for
each cohort in the trial. Given a set of simple orderings that are pre-
specified, POCRM selects the most likely ordering and uses this to recommend
the next dose to be assigned in the trial via the CRM. This approach is
particularly favourable as it allows for flexibility in the dose-toxicity
ordering used for dose recommendation. However, it also significantly limits
the performance of the model as, particularly in cases with a large number of
dose combinations, not all possible orderings can be considered by the model
since the model can only consider orderings given as candidates. It can be
argued that the true MTD will be selected regardless of whether the correct
dose ordering is present, however, since the proposed orderings alter toxicity
estimates, they also guide dose-escalation. In particular, in the combination
setting, multiple doses that have a toxicity close to the TTR may be present.
Therefore, inaccurate estimation of the risk of toxicity would hide
alternative true MTDs. In this case, not only is escalation important, but
also point estimation. Another instance under which this method is potentially
problematic is where multiple orderings have a similar posterior probability
leading to the uncertainty in dose orderings being disregarded as only the
model with the highest posterior probability is selected at each step. This
article also explores the prevalence of illogical ’jumps’ in dose-toxicity
estimates present in practice, largely due to changes in the preferred simple
ordering.

To address the challenges associated with ordering selection and uncertainty
quantification during trials, the original POCRM is extended by applying
Bayesian model averaging (BMA) in this work. Previously, BMA has been applied
to the original CRM (Yin and Yuan, 2009) where it was highly successful in
improving dose allocation for single-agent trials, particularly for small
sample Phase I trials. The novel BMA-POCRM design aims to incorporate
uncertainty in the toxicity ordering with the aim of making more flexible
dose-toxicity estimates, which are not limited to following a single
predefined ordering, and take into account the additional uncertainty implied
by a partial ordering.

Clinician trust is an essential aspect of running real-world adaptive clinical
trials. Even if a clinical trial design has good statistical properties, as
confirmed via a simulation study, the design will not be implemented in
practice if its recommendations are not aligned with the evidence from the
trial. A lack of trust in the design can lead to the recommendations given by
the model not being followed despite existing evidence that may support the
model’s decision. The concept of estimation coherency, proposed and developed
in this manuscript, aims to provide an additional metric for evaluating design
behaviour in real-world settings. Via a case study and an extensive simulation
study, it will be demonstrated that this property is crucial to ensuring that
the dose-escalation and de-escalation recommendations of a design can be
easily communicated to the clinical team on the trial and, hence, more likely
to be adhered to. This will result in delivering more efficient early phase
dose-escalation trials.

The general framework of the Bayesian CRM is outlined along with the POCRM and
the novel method BMA-POCRM in Section 2. Section 3 provides clear motivation
for the development of BMA-POCRM. Here, the novel concept of estimation
incoherency is introduced, which measures the consistency of dose toxicity
estimate updates with respect to the given dose toxicity orderings.
Furthermore, we provide a case study of the performance of BMA-POCRM and POCRM
on real trial data in Section 4. Together, the evidence from the case study
and further simulation results in Section 5 show that BMA-POCRM improves the
accuracy, safety and operating characteristics with more intuitive escalation
and de-escalation decisions. An examination of the operating characteristics
is carried out in Section 5 leading to a sensitivity analysis of results under
various cohort sizes in Section 6. Finally, a discussion and analysis of the
results are presented in Section 7.

##  2 Methods

###  2.1 General Framework

Consider a setting with a partial ordering corresponding to $M$ simple
orderings and $K$ dose levels, $\\{d_{1},\ldots,d_{K}\\}$. Following a
framework similar to that set by the original CRM (O’Quigley et al., 1990),
let $X_{j}$ be the dose level assigned to the $j$-th patient where
$x_{j}\in\\{d_{1},\ldots,d_{K}\\}$ and let $Y_{j}$ be a binary random variable
for whether patient $j$ experiences a DLT.

For a particular ordering $m\in\\{1,\ldots,M\\}$, the risk of DLT at
$d_{k}\in\\{d_{1},\ldots,d_{K}\\}$ is modelled as,

$\widehat{R}(d_{k})=\Pr[Y_{j}=1|X_{j}=d_{k}]=\psi_{m}(d_{k},a),$

where $\widehat{R}(d_{k})$ is the estimated risk of DLT at $d_{k}$, $\psi_{m}$
is the working model under ordering $m$ (Wages et al., 2011), and $a$ is the
model parameter. A wide range of specifications can be selected for the
working model, each with their own associated assumptions regarding the dose-
toxicity relationship and parameter estimation approaches (Cheung and
Chappell, 2002). A necessary assumption of the working model is that the
relationship between dose level and dose toxicity is monotonic under a given
simple ordering. This allows for multiple models to be specified based on the
defined set of simple orderings, with each model corresponding to a specific
ordering of doses.

The dose level for the next cohort of patients is allocated by minimising the
difference between estimated risk of DLT and the TTR, which can be expressed
as the following criterion,

$x_{j+1}=\operatorname*{arg\,min}_{d_{k}}|\widehat{R}(d_{k})-\theta|,$ (1)

where $\theta$ is the TTR. By repeating this estimation-minimisation process
until the stopping conditions are satisfied, an estimate for the true MTD is
obtained, which is the dose recommended by the model following the final
cohort of patients.

###  2.2 Continual Reassessment Method for Partial Ordering (POCRM)

Under the Bayesian framework of the POCRM, a prior distribution for the model
parameters, $f(a)$, and a prior probability for each ordering $p(m)$ where
$\sum\limits_{m=1}^{M}p(m)=1$ and $p(m)\geq 0\;\forall m$ is required.

Since $Y$ is binary, the likelihood takes the form of a Bernoulli random
variable, where each patient in a cohort either experiences or does not
experience a DLT. The observed data up to patient $j$ is defined as
$\Omega_{j}=\\{x_{1},y_{1},\ldots,x_{j},y_{j}\\}$. This gives the following
likelihood under ordering $m$ after the inclusion of $j$ patients in the
trial,

$L_{m}(a|\Omega_{j})=\prod_{l=1}^{j}\\{\psi_{m}(x_{l},a)\\}^{y_{l}}\\{1-\psi_{m}(x_{l},a)\\}^{1-y_{l}},$
(2)

where $x_{l}$ is the dose allocated to patient $l$, and $y_{l}$ is the binary
variable denoting whether patient $l$ experiences a DLT and $\Omega_{j}$
contains the paired patient data $(x_{l},y_{l})$. Given the likelihood, the
posterior density for the parameter $a$ under the $m$-th model is given by

$f_{m}(a|\Omega_{j})=\frac{L_{m}(a|\Omega_{j})f(a)}{\int_{\mathcal{A}}L_{m}(a|\Omega_{j})f(a)da}.$
(3)

That is, each model will have a unique posterior distribution for $a$, as each
model implies a unique ordering of doses. The posterior probabilities for each
ordering are also obtained via Bayes’ rule as follows,

$p(m|\Omega_{j})=\frac{p(m)\int_{\mathcal{A}}L_{m}(a|\Omega_{j})f(a)da}{\sum\limits_{m^{\prime}=1}^{M}p(m^{\prime})\int_{\mathcal{A}}L_{m^{\prime}}(a|\Omega_{j})f(a)da}.$
(4)

The model used for dose allocation is selected by maximising the posterior
model probabilities,

$m^{*}=\operatorname*{arg\,max}_{m}p(m|\Omega_{j}),$ (5)

which results in a single partial ordering being selected for downstream
estimation. The posterior density corresponding to this model is then used to
estimate a posterior mean for parameter $a$,

$\widehat{a}_{m^{*}}=\int_{\mathcal{A}}af_{m^{*}}(a|\Omega_{j})da,$ (6)

which can be plugged directly into the working model to obtain an estimate for
the risk of DLT for the $k$-th dose,

$\widehat{R}(d_{k})=\psi_{b}(d_{k},\widehat{a}_{m^{*}}),$ (7)

from which the next dose is allocated using the criterion expressed in
Equation (1).

###  2.3 Bayesian Model Averaging POCRM (BMA-POCRM)

Suppose that rather than selecting a single model or ordering for each dose
allocation in the trial, all orderings are taken into account before making
the next decision on dose allocation. Bayesian model averaging (BMA) (Raftery
et al., 1997) would allow for estimates for probability of toxicity under
multiple orderings to be combined for a single probability across several
models. This is done by combining the available posterior information on the
set of parameters $a$ and the $M$ model probabilities. BMA can be applied to
the posterior distribution of the set of parameters $a$. This allows for
incorporating uncertainity in $a$ and partial orderings before obtaining an
estimate of the risk of DLT of each dose level.

Using the previously described model-specific posterior distributions of $a$
and the posterior model probabilities $p(m|\Omega_{j})$, a weighted
combination of posterior distributions is obtained as follows,

$g(a|\Omega_{j})=\sum\limits_{m=1}^{M}p(m|\Omega_{j})f_{m}(a|\Omega_{j}),$ (8)

where $g(a|\Omega_{j})$ is the model averaged posterior of $a$. Hence, rather
than having a model-specific posterior for $a$, a version averaged across all
models, relative to the posterior probabilities $p(m|\Omega_{j})$ is now
obtained.

By applying a change of variables to $g(a|\Omega_{j})$ the working model for
$R(d_{k})$ is expressed as a probability distribution. Although $R(d_{k})$ is
a useful expression for interpretation, $R(d_{k})$ for each $k$ is estimated
independently of the others, hence, $R(d_{k})$ is not a function of $d_{k}$
but a function of $a$ since $R(d_{k})=\psi_{m}(d_{k},a)$. The change in
variable is $R(d_{k})=\psi_{m}(d_{k},a)$ such that from the combined posterior
distribution the following probability distribution function for $R(d_{k})$
can be obtained,

$|f_{m}\left(R(d_{k}|\Omega_{j})\right)dR(d_{k})|=|f_{m}(a|\Omega_{j})da|,$
(9)

since the probability contained under a differential area must be remain
constant under a change of variables. Thus,

$f_{m}(R(d_{k})|\Omega_{j})=\left|\frac{d\psi^{-1}_{m}(d_{k},R(d_{k}))}{dR(d_{k})}\right|f_{m}(a|\Omega_{j}).$
(10)

Reapplying BMA the following combined posterior distribution for the risk of
toxicity is determined,

$g(R(d_{k})|\Omega_{j})=\sum_{m=1}^{M}p(m|\Omega_{j})f_{m}(a|\Omega_{j}),$
(11)

which is independent of ordering. The expectation for the risk of toxicity
under dose $d_{k}$ is then,

$\mathbb{E}[R(d_{k})]=\int_{0}^{1}R\left(d_{k})\;g_{\psi(d_{k},a)}(R(d_{k})|\Omega_{j}\right)\;d\\{R(d_{k})\\}$
(12)

which under the law of the unconscious statistician (Papoulis and Unnikrishna
Pillai, 2002) is equivalent to

$\mathbb{E}[R(d_{k})]=\int_{\mathcal{A}}\psi(d_{k},a)g(a|\Omega_{j})da,$ (13)

where $\mathcal{A}$ is the domain of $a$. Finally, dose allocation is carried
out under this framework by setting $\widehat{R}(d_{k})=\mathbb{E}[R(d_{k})]$
and applying the criterion described in Section 2.1.

##  3 Coherency in the Presence of Partial Ordering

###  3.1 Defining Estimation Coherency

Coherence in Phase I clinical trials is a useful concept for assessing the
theoretical qualities of trial methodology. Throughout this article,
coherence, which is defined based on escalation and de-escalation behaviour of
a dose-finding model, is referred to as _escalation coherence_. In practice,
an escalation coherent design will benefit patient safety as it reduces the
likelihood of assigning an overly toxic dose to a patient or cohort but will
also ensure that the maximal dose level considered safe is administered.
Expanding to the drug combination setting, Park and Liu (2020) introduces
definitions of strong and weak coherency, both of which rely on evaluating the
characteristics of escalation and de-escalation to define coherency. A
definition of escalation coherency which is equally applicable to both single-
agent and dual-agent combination trials is also presented.

There are two sets of doses that are relevant to escalation coherency. Let
$\mathcal{E}_{n}$ and $\mathcal{D}_{n}$ contain the candidate dose levels for
escalation and de-escalation, respectively, for dose allocation $X_{n}$.

######  Definition 3.1 (Escalation Coherency).

A design is coherent in dose escalation if
$\Pr[X_{n+1}\in\mathcal{E}_{n}|Y_{n}=1]=0$ for $n=1,\ldots,d-1$ and is
coherent in dose de-escalation if $\Pr[X_{n+1}\in\mathcal{D}_{n}|Y_{n}=0]=0$
for $n=1,\ldots,d-1$. A design is escalation coherent if dose escalation and
de-escalation are coherent.

This definition emphasises the sole concept of dose selection. Whereas, during
the administration of a real-world trial, there is frequent interface between
domain expert and dose-escalation model. Domain experts use both the given
toxicity estimates and the recommended next dose to make a final escalation
decision. This is exceedingly relevant to the combination setting where there
could be several combinations to choose from at any given point in the trial.
To address the need for both toxicity estimates and dose escalation
recommendations to be coherent, the following specification of _estimation
coherency_ is introduced.

######  Definition 3.2 (Estimation Coherency).

Suppose there is a given partial ordering, from which $M$ simple orderings are
derived. Let $I_{m}(d_{i})$ be the index of dose $i$ within ordering $m$. If
dose $i$ is more toxic than dose $j$ under ordering $m$ then
$I_{m}(d_{j})<I_{m}(d_{i})$ and vice versa. Based on the $M$ orderings, each
dose $i$ will have an a priori known set of less toxic doses,

$\nu_{i}=\\{d_{j}:I_{m}(d_{j})<I_{m}(d_{i})\;\forall m\\},$

and similarly a known set of more toxic doses

$\xi_{i}=\\{d_{j}:I_{m}(d_{j})>I_{m}(d_{i})\;\forall m\\}.$

A design is estimation coherent if it satisfies both (i) following no DLT at
$d_{i}$ the estimated risk of toxicity for dose $j$, $\widehat{R}(d_{j})$,
decreases for all $d_{j}\in\nu_{i}\cup\xi_{i}$, and (ii) following DLT at
$d_{i}$, the estimated risk of toxicity for $d_{j}$, $\widehat{R}(d_{j})$,
increases for all $d_{j}\in\nu_{i}\cup\xi_{i}$.

This definition of estimation coherency acts as a ’sanity check’ for the dose-
finding model. By ensuring the model coincides with both the prior knowledge
implied by the partial ordering and most recently gained information, it is
able to detect changes which are potentially illogical and may endanger the
integrity of the trial. This type of coherency is not relevant to single dose-
finding methods as, under monotonicity, this condition is always satisfied.

###  3.2 Illustrative Example

The POCRM is particularly vulnerable to estimation incoherencies. Here, a
specific example of estimation incoherency under the POCRM is explored.

Consider the $3\times 2$ drug combination design shown in Table 1. Under this
combination setting, the chosen cohort size is 1 with a TTR of 0.4. Applying
the partial ordering specification recommended by Wages and Conaway (2013),
the following 6 simple orderings are obtained, of which 5 are unique,

$\displaystyle m=1:\;d_{1}\rightarrow d_{2}\rightarrow d_{3}\rightarrow
d_{4}\rightarrow d_{5}\rightarrow d_{6},$ $\displaystyle
m=2:\;d_{1}\rightarrow d_{3}\rightarrow d_{5}\rightarrow d_{2}\rightarrow
d_{4}\rightarrow d_{6},$ $\displaystyle m=3:\;d_{1}\rightarrow
d_{3}\rightarrow d_{2}\rightarrow d_{5}\rightarrow d_{4}\rightarrow d_{6},$
$\displaystyle m=4:\;d_{1}\rightarrow d_{2}\rightarrow d_{3}\rightarrow
d_{4}\rightarrow d_{5}\rightarrow d_{6},$ $\displaystyle
m=5:\;d_{1}\rightarrow d_{2}\rightarrow d_{3}\rightarrow d_{5}\rightarrow
d_{4}\rightarrow d_{6},$ $\displaystyle m=6:\;d_{1}\rightarrow
d_{3}\rightarrow d_{2}\rightarrow d_{4}\rightarrow d_{5}\rightarrow d_{6}.$

This set of orderings is the complete set of possible orderings, assuming that
dose-toxicities increase monotonically only where the dose level increases in
only one drug of the combination. At the start of the trial, the _a priori_
probability for each ordering is equal.

From these simple orderings one can obtain the sets of interest $\nu_{i}$ and
$\xi_{i}$ for each dose as shown for each dose level in Table 2. The sets used
to check for estimation incoherencies are composed as follows. For $d_{2}$,
$\nu_{2}=\\{d_{1}\\}$ and $\xi_{2}=\\{d_{4},d_{6}\\}$. This is obtained by
considering that $d_{1}$ is less toxic than $d_{2}$ under every simple
ordering, and both $d_{4}$ and $d_{6}$ are more toxic than $d_{2}$ under every
simple ordering. These sets exist for each dose level and can be used
following each Bayesian probability update to detect any estimation
incoherencies.

Dose Level, $i$ | $\nu_{i}$ | $\xi_{i}$
---|---|---
1 | $\emptyset$ | $\\{d_{2},d_{3},d_{4},d_{5},d_{6}\\}$
2 | $\\{d_{1}\\}$ | $\\{d_{4},d_{6}\\}$
3 | $\\{d_{1}\\}$ | $\\{d_{4},d_{5},d_{6}\\}$
4 | $\\{d_{1},d_{2},d_{3}\\}$ | $\\{d_{6}\\}$
5 | $\\{d_{1},d_{3}\\}$ | $\\{d_{6}\\}$
6 | $\\{d_{1},d_{2},d_{3},d_{4},d_{5}\\}$ | $\emptyset$
Table 2: Sets for detecting estimation incoherencies where $\nu_{i}$ is the
set of doses always less toxic than dose $i$ and $\xi_{i}$ is the set of doses
always more toxic than dose $i$.

###  3.3 Applying POCRM with Model Selection

The parametrisation of POCRM used here matches that which is used in the full
simulation study as described in Section SM.2 of the Supplementary Materials
with a cohort size of 1. 1,000 simulations were initiated with the first
estimation incoherency identified in the second simulated trial at the
induction of cohort 12. At this point in the trial, the posterior ordering
probabilities are shown in Table 3. The dose allocations and DLTs observed up
to and including cohort 11 are as follows,

$\displaystyle\mathbf{n}^{[11]}$ $\displaystyle=(1,0,1,6,2,1),$
$\displaystyle\mathbf{y}^{[11]}$ $\displaystyle=(0,0,0,3,1,1),$

where the $i$-th entry in $\mathbf{n}\in\mathbb{Z}^{6}$ is the number of
patients assigned to the $i$-th dose level, $d_{i}$, and
$\mathbf{y}\in\mathbb{Z}^{6}$ is the number of patients that experienced a DLT
after being assigned dose level $d_{i}$. As seen in Table 3, following the
induction of cohort 11, $d_{2}$ is recommended as the next dose by POCRM.
Hence, cohort 12 is inducted and administered $d_{2}$, which yields the
following allocation vectors,

$\displaystyle\mathbf{n}^{[12]}$ $\displaystyle=(1,1,1,6,2,1),$
$\displaystyle\mathbf{y}^{[12]}$ $\displaystyle=(0,0,0,3,1,1),$

where no new DLT is observed for $d_{2}$.

Recalling the set of doses with known toxicity relative to $d_{2}$ in Section
3.2, since no DLT is observed at $d_{2}$ the dose-toxicity estimates for all
dose levels in $\nu_{2}\cup\xi_{2}=\\{d_{1},d_{4},d_{6}\\}$ are expected to
decrease. However, in Table 3 an increase in the toxicity estimate for $d_{4}$
from 0.49 to 0.56 is observed. This is a change of +0.07 despite there being
no information gained from the previous cohort that indicates a greater
toxicity of $d_{4}$. Throughout the coming analyses, the occurrence of
estimation incoherencies is considered as a key operating characteristic of
the methods being studied.

Method | Cohort |  $p(\cdot|x):$ (Posterior Model Probability)
---|---|---
$m=1$ | $m=2$ | $m=3$ | $m=4$ | $m=5$ | $m=6$
POCRM | 11 | 0.1568 | 0.1497 | 0.1878 | 0.1568 | 0.1582 |  0.1906
12 | 0.1743 | 0.1091 | 0.1840 | 0.1743 |  0.1840 | 0.1743
BMA-POCRM | 11 | 0.1568 | 0.1497 | 0.1878 | 0.1568 | 0.1582 | 0.1906
12 | 0.1743 | 0.1091 | 0.1840 | 0.1743 | 0.1840 | 0.1743
Method | Cohort |  $R(\cdot):$ (Estimated Probability of Toxicity)
$d_{1}$ | $d_{2}$ | $d_{3}$ | $d_{4}$ | $d_{5}$ | $d_{6}$
POCRM | 11 | 0.0672 |  0.3261 | 0.1756 | 0.4859 | 0.6282 | 0.7412
12 | 0.0331 | 0.1114 | 0.2432 | 0.5562 |  0.4023 | 0.6853
BMA-POCRM | 11 | 0.0802 | 0.2671 | 0.2371 | 0.5247 |  0.5019 | 0.7111
12 | 0.0654 | 0.2281 | 0.2210 | 0.4975 |  0.4860 | 0.6933
Table 3: Model probabilities and dose-toxicity estimates for POCRM and BMA-
POCRM dose-escalation frameworks in a setting with a TTR of 0.4. The
recommended dose levels and selected orderings following each cohort are shown
with an underline. Incoherencies are shown in red. Note: Despite, the
recommended dose-level given by BMA-POCRM being $d_{5}$ following cohort 11,
$d_{2}$ is treated as the recommended dose to maintain consistency for
comparison of POCRM and BMA-POCRM.

This example illustrates the often illogical changes in toxicity estimates
observed under the POCRM. Since toxicity estimates guide clinicians and affect
the next allocated dose, it is crucial that changes in these toxicity
estimates are robust to scrutiny.

###  3.4 Applying BMA-POCRM

Again, under this setting, the dose-toxicity estimates for all dose levels in
$\nu_{2}\cup\xi_{2}=\\{d_{1},d_{4},d_{6}\\}$ are expected to decrease or
remain the same as there was no DLT observed for $d_{2}$. Returning to the
results presented in Table 3 for BMA-POCRM, for all doses in
$\nu_{2}\cup\xi_{2}$ the toxicity estimate decreases. Conversely to POCRM,
BMA-POCRM remains coherent with respect to estimation coherency in this case.

##  4 Case Study

To consider the performance of BMA-POCRM in comparison to POCRM in a real-
world setting, we apply these methods to a Phase I study dosing patients
combinations of neratinib and temsirolimus (Gandhi et al., 2014). This trial
involved fifty-two patients treated on 12 doses in a 4-by-4 grid of possible
neratinib-temsirolimus combinations with a TTR of $\theta=1/3$. The remaining
4 doses in the grid were never assigned to patients in the original trial. A
DLT is defined as an inability to maintain the prescribed dose for the first
28 days of treatment due to treatment-related toxicity. The two initial
cohorts consisted of two patients each were enrolled simultaneously with: (1)
160 mg of neratinib/15 mg of temsirolimus; and (2) 120 mg of neratinib/25mg of
temsirolimus.

We apply BMA-POCRM and POCRM the parameters selected in the calibration study
in Section SM.2 and follow the approach of Wages and Conaway (2013) to derive
candidate orderings based on the 12 available doses in the trial. The full
trial data used for these simulations can be found in Table S3 in the
Supplementary Materials.

The aim here is to identify incoherencies and large changes in toxicity
estimates. By showing that these occur in real-world trials we aim to further
motivate the use of BMA-POCRM in practice.

###  4.1 Data Generation

To allow for a fair evaluation of model-guided dose escalation, we use the
original trial data to conduct a case study according to the scheme outlined
by Barnett et al. (2024). We define a fixed set of fifty-two patient dose
responses for each dose. Here, we denote the number of patients assigned to
dose $j$ by $n_{j}$ and the number of observed DLTs under dose $j$ by $y_{j}$.

To define a fixed set of fifty-two patient dose responses for each dose we
take the first $n_{j}$ responses to be a random permutation of responses from
the original study. The remaining $52-n_{j}$ responses are generated from a
Beta$(1+y_{j},1+n_{j}-y_{j})$ distribution by first sampling a probability of
DLT and sampling a binary response from a Bernoulli distribution with the
given probability of DLT. Where no patients are assigned to a dose combination
in the real study, probabilities are generated from a Beta$(3,3)$ distribution
instead. As each method allocates a dose to each cohort, the ordered set of
responses corresponding to each dose is used to determine treatment response
sequentially. This ensures that the $n$th patient allocated to each dose,
regardless of method of dose escalation used, will have the same response to
treatment.

###  4.2 Results

Figure 1 shows that the range of magnitude of changes in toxicity estimate for
POCRM is larger than for BMA-POCRM. The range is shown to be [-0.37, 0.35] and
[-0.17, 0.18] for POCRM and BMA-POCRM, respectively.

Table 4 shows an example of a cohort where the large change in toxicity
estimate is particularly great. After assigning $d_{8}$ to cohort 3, we
observe a DLT, which leads to several incoherencies and large changes in
toxicity estimates for the estimates that are to be used to assign doses to
cohort 4. For POCRM, we observe two changes that are $>25$, for $d_{4}$ and
the administered dose $d_{8}$. We also observe two estimate incoherencies for
doses $d_{5}$ and $d_{6}$, which we know are less toxic than $d_{8}$ according
to the partial orderings corresponding to this trial setting. For this cohort,
the toxicity estimates given by BMA-POCRM do not yield any incoherencies or
large changes as is true for the remainder of the motivating trial, as shown
in Figure 2.

Figure 1: Distribution of the magnitude of changes in toxicity estimates in
the motivating trial simulation. Whiskers indicate range of maximum and
minimum values observed.

$R(\cdot):$ (Estimated Probability of Toxicity) $d_{\text{sel}}$
$y_{\text{sel}}$ Method Cohort $d_{1}$ $d_{2}$ $d_{3}$ $d_{4}$ $d_{5}$ $d_{6}$
$d_{7}$ $d_{8}$ $d_{9}$ $d_{10}$ $d_{11}$ $d_{12}$ 3 0.00 0.00 0.01 0.04 0.08
0.15 0.23 0.33 0.43 0.53 0.61 0.69 $d_{8}$ 1 POCRM 4 0.00 0.01 0.07 0.32 0.03
0.14 0.42 0.68 0.22 0.51 0.75 0.60 – – 3 0.01 0.04 0.13 0.25 0.06 0.13 0.29
0.48 0.18 0.36 0.52 0.41 $d_{10}$ 0 BMA-POCRM 4 0.01 0.02 0.09 0.20 0.03 0.08
0.22 0.41 0.11 0.27 0.44 0.32 – –

Table 4: Example of a sudden change in toxicity estimates under POCRM in the
motivating trial. Any sudden changes (i.e. those greater than 0.25) are shown
in  red and any incoherencies are shown in  blue.

Figure 2 shows that incoherencies only occurred for POCRM in the motivating
trial. The incoherencies shown for step 3-4 correspond to the results in Table
4. The full dose allocations can be found in the Appendix in Tables S4 and S5.

Figure 2: Incoherencies that occurred in motivating trial.

##  5 Simulation Study

###  5.1 Specification

A simulation study is necessary to evaluate and compare the operating
characteristics for POCRM and BMA-POCRM in a large-scale setting. For this
study, a 4-by-4 structure is selected in a 2 drug combination study. For each
scenario, $10,000$ trial simulations are run with a TTR of $\theta=0.3$ and a
cohort size of $1$ for $60$ cohorts. In Section 6, cohort sizes of 3 and 6 are
considered with 20 and 10 cohorts, respectively. We select model parameters
based on the calibration study presented in Section SM.2.

There are several candidates for a working model $\psi_{m}(d_{k},a)$. The
parametrisation of the power model given in Wages et al. (2011) is used here
and takes the following form,

$\psi_{m}(d_{k},a)=\alpha_{mk}^{a},\;\;\;\;k=1,\ldots,K$

where $a\in[0,\infty)$ and $0<\alpha_{m1}<\ldots,<\alpha_{mK}$ is the
probability skeleton, which represents the prior estimates of dose toxicity at
each dose level under ordering $m$. The prior distribution of $a$ is Normal
with mean 0 and variance 1.34 as suggested by Wages et al. (2011). The
selected variance of the normal prior has no impact on model performance as
shown in SM.3. This working model is applied throughout the paper for
implementations of both POCRM and BMA-POCRM.

Five metrics are used to assess models in simulation trials according to
recommendation accuracy, assignment accuracy and patient safety. In this case,
an overly toxic dose is chosen to be any dose that has a true probability of
toxicity greater than $110\%$ of the TTR.

  1. (i)

The proportion of correct selections (PCS) is the proportion of trials that
recommended doses with a true probability of toxicity equal to the TTR.

  2. (ii)

The proportion of acceptable selections (PAS) is the proportion of trials that
recommended a dose with a true probability of toxicity within
$[\theta-0.1,\theta]$ where $\theta$ is the TTR.

  3. (iii)

Proportion of trials that give overly toxic selections (POTS).

  4. (iv)

The number of patients treated at overly toxic doses (NPTOT).

  5. (v)

Finally, the estimation coherency of model estimates at the induction of each
trial as defined in Definition 3.2.

The full specification for all 24 scenarios can be found in Table S1 of the
Supplementary Materials. Scenarios with both symmetric and asymmetric true MTD
locations are included. The number of accurate orderings present in the set of
candidate orderings is also varied across scenarios. Those found in the left-
most column of Table S1 have no correct orderings, implying that none of the
orderings exactly match the specified probabilities of toxicity. Scenarios
found in the middle column are satisfied by a single ordering in the list of
candidate orderings. Finally, those found in the right-most column are equally
satisfied by 2 candidate orderings.

###  5.2 Results

Figure 3 illustrates model performance under all $24$ scenarios, and provides
an arithmetic mean performance across scenarios. BMA-POCRM exhibits better PCS
and PAS across all scenarios but scenario 15, indicating that it more
consistently selects desirable dose levels following the Phase I trial
simulations.

Figure 3: Panel A. Arithmetic mean PCS across 10,000 repeated simulations with
a cohort size 1 for each scenario. The mean performance across scenarios is
also shown. Panel B. Arithmetic mean PAS across 10,000 repeated simulations
with a cohort size 1 for each scenario. The mean performance across scenarios
is also shown. Panel C. Arithmetic mean POTS across 10,000 repeated
simulations with a cohort size 1 for each scenario. The mean performance
across scenarios is also shown. Panel D. Arithmetic mean PTOTD across 10,000
repeated simulations with a cohort size 1 for each scenario. The mean
performance across scenarios is also shown.

BMA-POCRM leads to an average increase of 5.2% in PCS when compared to that of
POCRM. Across all scenarios with non-zero PCS, there is an improvement in
performance when applying BMA-POCRM. Considering scenarios 1-13, no scenario
shows a lower difference than 2.77% in favour of BMA-POCRM. Scenario 12 shows
the greatest discrepancy between the methods with BMA-POCRM exceeding POCRM by
10.7%. Likewise, investigating Figure 3.B, BMA-POCRM leads to an improvement
by at least 3.34% and an average of 5.5%. For this measure, scenario 12 leads
to the greatest difference in performance with 10.5%. Differences in the
number of correct orderings between scenarios does not clearly favour either
method and there is no consistent effect on performance associated with
changes in the number of correct orderings.

The design matrices for scenarios 14 and 15 are such that there are only
overly toxic doses and doses with a lower toxicity than the TTR, respectively.
Under scenario 15, POCRM has a mean PAS of $85.61\%$ whilst BMA-POCRM has one
of $84.67\%$. Nearly half of the doses in scenario 15 are considered
acceptable by definition. Noting that these doses are in the lower triangular
of the combination matrix, these are the most toxic doses under both drugs.
This suggests that BMA-POCRM is the more conservative approach to dose
escalation.

The safety of BMA-POCRM is further supported by Figures 3.C and 3.D, where it
recommends overly toxic doses in a smaller proportion of trials and allocates
overly toxic doses to fewer patients, respectively. This trend is maintained
throughout all scenarios as the mean POTS is 4.89% lower for BMA-POCRM than
POCRM. Scenario 2 exhibits the smallest difference between methods in terms of
POTS with a 1.84% decrease in the number of overly toxic selections with BMA-
POCRM. In scenario 12, BMA-POCRM leads to a 10% reduction in the number of
overly toxic doses. Figure 3.D shows that on average, at least one cohort less
will be assigned an overly toxic dose.

Figure 4.A shows the magnitude and direction of changes in toxicity estimates
across all trial steps in the simulated trials. From 4.A, BMA-POCRM has less
dispersed both in terms of the box and whiskers for each scenario. Conversely,
for changes in model probabilities shown in Figure S3, since the scheme for
obtaining posterior model probabilities is shared between BMA-POCRM and POCRM,
the distributions are similar.

The frequency of incoherencies within simulated trials is shown in Figure 4.B.
The previous definition of estimation incoherency is here applied to evaluate
the estimation operating characteristics of POCRM and BMA-POCRM. In all
scenarios but 13, 14 and 15, POCRM leads to at least $90\%$ of trials with at
least one estimation incoherent change of toxicity estimates during the trial.
Meanwhile, BMA-POCRM exhibits estimation incoherencies in a worst-case
$3.13\%$ of trials under scenario 11 and no incoherencies in most scenarios.
Further, examining Figure 4.C under scenario 11, the magnitude of estimation
incoherent toxicity estimates is significantly smaller for BMA-POCRM with a
maximum incoherent change in toxicity estimate of $0.060$ compared to POCRM
for which this is $0.288$. Figure 4.D shows that in trials where BMA-POCRM
does exhibit estimation incoherencies, this occurs for a smaller number of
cohorts. The distribution of the number of cohorts where BMA-POCRM exhibits an
incoherency is consistently biased towards zero when compared to those of
POCRM.

Figure 4: Panel A. Distribution of changes in dose toxicity estimates
following the induction of each cohort for simulated trials with cohort size
1. Panel B. Proportion of trials where at least one incoherent estimate is
observed. Panel C. Magnitude of incoherent changes in toxicity estimates.
Panel D. Mean number of cohorts with at least one observed incoherency.

##  6 Sensitivity Analysis

The operating characteristics of BMA-POCRM illustrated in the previous section
show that for a cohort size of 1, this method outperforms the existing POCRM.
In this section, the analysis conducted for cohorts of 1, are expanded to
cohort sizes of 3 and 6. The simulation settings remain the same as those
described in Section 5. Moreover, the calibrated designs selected in Figure S2
are also used. With a cohort size of 3 and 6, the number of cohorts falls to
20 and 10, respectively, to maintain a constant number of patients at trial
completion.

Figure 5 shows that there is no significant change in performance when
increasing cohort size to 3. Specifically, in Figure 5.A-B, BMA-POCRM has mean
PCS of 50.42% compared to that of POCRM of 45.41%, a difference of 5.01%. For
PAS, mean performance was 64.34% and 59.12% for BMA-POCRM and POCRM
respectively, a difference of 5.25%. These discrepancies are consistent with
those observed for a cohort size of 1. Furthermore, when further increasing
cohort size to 6, there is a gap of 6.44% and 6.59% in mean PCS and mean PAS,
respectively. For both cohort sizes 3 and 6 under scenario 15, where POCRM
previously had superior PAS, BMA-POCRM is now dominant.

Figure 5: Panel A. Geometric mean of PCS for simulated trials with cohort size
3 and 20 cohorts. Panel B. Geometric mean of PCS for simulated trials with
cohort size 6 and 10 cohorts. Panel C. Geometric mean of PAS for simulated
trials with cohort size 3 and 20 cohorts. Panel D. Geometric mean of PAS for
simulated trials with cohort size 6 and 10 cohorts. Panel E. Geometric mean of
PTOS for simulated trials with cohort size 3 and 20 cohorts. Panel F.
Geometric mean of PTOS for simulated trials with cohort size 6 and 10 cohorts.
Panel G. Geometric mean of NPTOT for simulated trials with cohort size 3 and
20 cohorts. Panel H. Geometric mean of NPTOT for simulated trials with cohort
size 6 and 10 cohorts.

The gains in trial safety also remain consistent with previous simulations as
shown in Figures 5.E-G. BMA-POCRM recommends fewer overly toxic doses and a
smaller number of patients are exposed to overly toxic doses relative to
POCRM. For a cohort size of 3, the greatest improvement in safety is seen for
scenario 11 where 8.66% of trials result in an overly toxic dose being
selected. Across scenarios, the average difference in POTS between methods was
4.3%. For scenario 22, 5.6 fewer patients were exposed to overly toxic doses
as a result of BMA-POCRM. Under a cohort size of 6, a mean decrease in POTS of
5.4% is observed with BMA. The greatest difference in POTS is observed for
scenario 8A where BMA-POCRM selects an overly toxic dose in 10.1% fewer
trials. The average number of patients exposed to overly toxic doses also
decreases from 21.9 to 18.4, with the sharpest difference being for scenario
9A with 6.5 fewer patients exposed to overly toxic doses. The rate of overly
toxic dose assignment was most similar under scenario 1 where there was only a
difference of 0.1 patients.

Figure 6.A shows the root mean square error (RMSE) of toxicity estimates at
the conclusion of trials for a cohort size of 3. In 19 out of 24 scenarios,
BMA-POCRM has toxicity estimates closer to the true values than POCRM. This is
most exagerated for scenario 13 where the mean RMSE of BMA-POCRM across
simulated trials is 0.0388 lower than for POCRM. In the 5 scenarios where
POCRM is a better estimator, the largest difference is of only 0.0075, a five-
fold smaller difference. For the average RMSE across all scenarios, shown in
Figure 6.D, BMA-POCRM is at 0.0871 whilst POCRM is at just 0.0969.

Figure 6: Panel A. Distribution of RMSE at the end of a simulated trial for
cohort size of 1. Panel B. Distribution of RMSE at the end of a simulated
trial for cohort size of 3. Panel C. Distribution of RMSE at the end of a
simulated trial for cohort size of 6. Panel D. Mean RMSE across scenarios for
all cohort sizes.

As with trials of cohort size 1, Figure 6.B shows that for trials with cohort
size of 3 and 6, BMA-POCRM has a lower RMSE for 19 and 20 scenarios
respectively. Scenario 13 led to the greatest improvement in RMSE with BMA-
POCRM having a 0.0391 and 0.0390 lower RMSE than POCRM for cohort sizes 3 and
6, respectively. The greatest difference in RMSE in favour of POCRM was for
scenario 11, where POCRM had a 0.00584 lower RMSE for a cohort size of 3.
Similarly, for cohort size of 6, POCRM had best performance relative to BMA-
POCRM for scenario 10 where it had an advantage of 0.0058. The negligible
differences in RMSE with changes to cohort size is due to the constant total
number of patients in each trial, with dose allocations varying depending on
cohort size. Hence, leading to some variability in this figure.

##  7 Discussion

In this paper, the concept of estimation coherency in Phase I clinical trials
for drug combinations is introduced. The application of the POCRM and its
operating characteristics in this combination setting are explored with the
proposition of a novel method, BMA-POCRM. This approach modifies POCRM to take
into account uncertainty in dose-toxicity ordering prior to making dose
escalations by applying Bayesian model averaging. To evaluate and compare the
POCRM and BMA-POCRM, a thorough study of their operating characteristics
spanning model accuracy, safety and robustness to estimation incoherency is
carried out. The novel method improved on its predecessor in its ability to
select correct and acceptable dose levels for progression to Phase II trials.
Moreover, the BMA-POCRM also led to fewer overly toxic dose recommendations.
The occurrence of estimation incoherencies significantly reduced with the BMA-
based method across all scenarios. For example, whilst for POCRM more than 90%
of trials exhibit incoherencies for 21 of 24 scenarios, the worst-case for
BMA-POCRM led to only $3.13\%$ of trials with incoherencies.

The study of operating characteristics shows a clear improvement in
performance achieved by BMA-POCRM. Combining estimates of the probability of
toxicity across several candidate orderings leads to greater flexibility in
potential predictions by the model. POCRM relies on selecting a single most
probable ordering on which it bases dose allocations for the current cohort.
Thus, BMA-POCRM has a distinct advantage in cases where there is no strong
prior knowledge of the underlying dose-toxicity orderings. In particular, the
BMA component of the model allows for the specification of intermediary
orderings not explicitly included in the original POCRM as shown by the
improved RMSE of toxicity estimates across cohort sizes. Despite this, even in
scenarios where the correct toxicity ordering is included as a candidate, BMA-
POCRM outperforms POCRM.

Despite the strong performance of the BMA-POCRM in the considered simulations,
there are several aspects of this model that have not yet been explored. The
novel approach does not completely eliminate the occurrence of estimation
incoherencies and the cause of these incoherencies is not clear. Furthermore,
a comparison study of the impact of prior information on BMA-POCRM and POCRM
performance is also necessary. This is particularly relevant when a correct
ordering is included as a candidate, prior information may be available to
inform these approaches, potentially leading to favourable performance for
POCRM.

Awknowledgements

This report is independent research supported by the National Institute for
Health Research (NIHR Advanced Fellowship, Dr Pavel Mozgunov, NIHR300576). The
views expressed in this publication are those of the authors and not
necessarily those of the NHS, the National Institute for Health Research or
the Department of Health and Social Care (DHSC). T. Jaki and P. Mozgunov
received funding from UK Medical Research Council (MC_UU_0000214 and
MC_UU_0000219, respectively). For the purpose of open access, the author has
applied a Creative Commons Attribution (CC BY) licence to any Author Accepted
Manuscript version arising.

## References

  * Barnett et al. (2024) Barnett, H., M. George, D. Skanji, G. Saint-Hilary, T. Jaki, and P. Mozgunov (2024).  A comparison of model-free phase i dose escalation designs for dual-agent combination therapies.  Statistical Methods in Medical Research, 09622802231220497.
  * Cheung (2019) Cheung, K. (2019).  dfcrm: Dose-Finding by the Continual Reassessment Method.  R package version 0.2-2.1.
  * Cheung and Chappell (2002) Cheung, Y. K. and R. Chappell (2002).  A simple technique to evaluate model sensitivity in the continual reassessment method.  Biometrics 58(3), 671–674.
  * Gandhi et al. (2014) Gandhi, L., R. Bahleda, S. M. Tolaney, E. L. Kwak, J. M. Cleary, S. S. Pandya, A. Hollebecque, R. Abbas, R. Ananthakrishnan, A. Berkenblit, et al. (2014).  Phase i study of neratinib in combination with temsirolimus in patients with human epidermal growth factor receptor 2–dependent and other solid tumors.  Journal of Clinical Oncology 32(2), 68–75.
  * Ivanova and Kim (2009) Ivanova, A. and S. H. Kim (2009).  Dose finding for continuous and ordinal outcomes with a monotone objective function: a unified approach.  Biometrics 65(1), 307–315.
  * Ivanova et al. (2003) Ivanova, A., A. Montazer-Haghighi, S. G. Mohanty, and S. D. Durham (2003).  Improved up-and-down designs for phase i trials.  Statistics in medicine 22(1), 69–82.
  * Ivanova and Wang (2004) Ivanova, A. and K. Wang (2004).  A non-parametric approach to the design and analysis of two-dimensional dose-finding trials.  Statistics in Medicine 23(12), 1861–1870.
  * Lee and Cheung (2009) Lee, S. M. and Y. K. Cheung (2009).  Model calibration in the continual reassessment method.  Clinical Trials 6(3), 227–238.
  * Mozgunov et al. (2020) Mozgunov, P., M. Gasparini, and T. Jaki (2020).  A surface-free design for phase i dual-agent combination trials.  Statistical methods in medical research 29(10), 3093–3109.
  * Mozgunov and Jaki (2019) Mozgunov, P. and T. Jaki (2019).  An information theoretic phase i–ii design for molecularly targeted agents that does not require an assumption of monotonicity.  Journal of the Royal Statistical Society. Series C, Applied Statistics 68(2), 347.
  * Mozgunov et al. (2022) Mozgunov, P., T. Jaki, I. Gounaris, T. Goddemeier, A. Victor, and M. Grinberg (2022).  Practical implementation of the partial ordering continual reassessment method in a phase i combination-schedule dose-finding trial.  Statistics in Medicine 41(30), 5789–5809.
  * O’Quigley et al. (1990) O’Quigley, J., M. Pepe, and L. Fisher (1990).  Continual reassessment method: a practical design for phase 1 clinical trials in cancer.  Biometrics, 33–48.
  * Papoulis and Unnikrishna Pillai (2002) Papoulis, A. and S. Unnikrishna Pillai (2002).  Probability, random variables and stochastic processes.
  * Park and Liu (2020) Park, Y. and S. Liu (2020).  On the coherence of model-based dose-finding designs for drug combination trials.  Plos one 15(11), e0242561.
  * Raftery et al. (1997) Raftery, A. E., D. Madigan, and J. A. Hoeting (1997).  Bayesian model averaging for linear regression models.  Journal of the American Statistical Association 92(437), 179–191.
  * Riviere et al. (2015) Riviere, M.-K., F. Dubois, and S. Zohar (2015).  Competing designs for drug combination in phase i dose-finding clinical trials.  Statistics in medicine 34(1), 1–12.
  * Wages and Conaway (2013) Wages, N. A. and M. R. Conaway (2013).  Specifications of a continual reassessment method design for phase i trials of combined drugs.  Pharmaceutical statistics 12(4), 217–224.
  * Wages et al. (2011) Wages, N. A., M. R. Conaway, and J. O’Quigley (2011).  Continual reassessment method for partial ordering.  Biometrics 67(4), 1555–1563.
  * Yin and Yuan (2009) Yin, G. and Y. Yuan (2009).  Bayesian dose finding in oncology for drug combinations by copula regression.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 58(2), 211–224.

SUPPLEMENTARY MATERIAL

SM.1 SIMULATION SCENARIOS

Drug A Dose Level 1 2 3 4 1 2 3 4 1 2 3 4 Scenario 1: 0 CO, 3 OT Scenario 2: 1
CO, 3 OT Scenario 3: 2 CO, 3 OT 1 0.01 0.04 0.16 0.16 0.01 0.04 0.10 0.16 0.01
0.04 0.10 0.13 2 0.07 0.13 0.20 0.30 0.07 0.13 0.20 0.30 0.07 0.10 0.16 0.30 3
0.10 0.20 0.30 0.45 0.16 0.20 0.30 0.40 0.10 0.20 0.30 0.40 4 0.20 0.30 0.40
0.50 0.20 0.30 0.45 0.50 0.20 0.30 0.45 0.50 Scenario 4: 0 CO, 6 OT Scenario
5: 1 CO, 6 OT Scenario 6: 2 CO, 6 OT 1 0.01 0.10 0.20 0.20 0.01 0.04 0.10 0.20
0.01 0.05 0.15 0.20 2 0.05 0.20 0.20 0.40 0.07 0.15 0.30 0.40 0.10 0.15 0.30
0.40 3 0.15 0.30 0.45 0.60 0.20 0.30 0.45 0.55 0.15 0.30 0.40 0.45 4 0.30 0.50
0.55 0.65 0.30 0.50 0.60 0.65 0.30 0.40 0.50 0.60 Scenario 7: 0 CO, 10 OT
Scenario 8: 1 CO, 10 OT Scenario 9: 2 CO, 10 OT 1 0.01 0.10 0.30 0.40 0.01
0.10 0.30 0.40 0.10 0.05 0.30 0.40 2 0.05 0.20 0.45 0.60 0.05 0.20 0.45 0.60
0.10 0.30 0.45 0.60 3 0.15 0.50 0.65 0.80 0.15 0.50 0.65 0.75 0.30 0.50 0.60
0.65 4 0.55 0.70 0.75 0.85 0.55 0.70 0.80 0.85 0.55 0.60 0.70 0.75 Scenario
10: 0 CO, 13 OT Scenario 11: 1 CO, 13 OT Scenario 12: 2 CO, 13 OT 1 0.15 0.30
0.50 0.55 0.15 0.20 0.40 0.55 0.15 0.20 0.40 0.45 2 0.20 0.45 0.60 0.75 0.30
0.45 0.60 0.75 0.30 0.40 0.50 0.65 3 0.40 0.65 0.80 0.96 0.50 0.65 0.80 0.90
0.40 0.55 0.65 0.70 4 0.70 0.85 0.90 1.00 0.70 0.85 0.95 1.00 0.60 0.65 0.75
0.80 Scenario 1A: 0 CO Scenario 2A: 1 CO Scenario 3A: 2 CO 1 0.01 0.10 0.13
0.15 0.01 0.05 0.10 0.18 0.01 0.05 0.14 0.17 2 0.05 0.13 0.15 0.30 0.10 0.15
0.20 0.30 0.10 0.14 0.20 0.30 3 0.08 0.15 0.20 0.45 0.18 0.20 0.40 0.44 0.14
0.20 0.40 0.45 4 0.30 0.40 0.50 0.55 0.30 0.40 0.47 0.60 0.30 0.40 0.50 0.55
Scenario 4A: 0 CO Scenario 5A: 1 CO Scenario 6A: 2 CO 1 0.01 0.10 0.15 0.30
0.05 0.15 0.30 0.50 0.05 0.15 0.30 0.45 2 0.05 0.20 0.40 0.55 0.30 0.40 0.55
0.75 0.30 0.40 0.50 0.65 3 0.30 0.45 0.60 0.75 0.45 0.55 0.75 0.80 0.40 0.55
0.65 0.70 4 0.50 0.65 0.70 0.80 0.60 0.80 0.85 0.90 0.60 0.65 0.75 0.80
Scenario 7A: 0 CO Scenario 8A: 1 CO Scenario 9A: 2 CO 1 0.05 0.30 0.40 0.50
0.05 0.15 0.30 0.45 0.05 0.30 0.40 0.45 2 0.15 0.40 0.55 0.65 0.30 0.40 0.50
0.60 0.15 0.40 0.50 0.60 3 0.30 0.55 0.65 0.80 0.50 0.50 0.65 0.70 0.30 0.50
0.60 0.70 4 0.60 0.70 0.75 0.85 0.55 0.70 0.75 0.80 0.55 0.60 0.75 0.80
Scenario 13: 15 OT Scenario 14 Scenario 15 1 0.30 0.40 0.50 0.57 0.37 0.42
0.50 0.57 0.05 0.07 0.17 0.19 2 0.40 0.48 0.57 0.68 0.40 0.48 0.57 0.68 0.07
0.12 0.19 0.24 3 0.45 0.61 0.72 0.87 0.43 0.61 0.72 0.87 0.10 0.19 0.23 0.25
Drug B 4 0.64 0.75 0.80 0.95 0.64 0.75 0.85 0.95 0.17 0.23 0.25 0.25

Table S1: scenario settings for all scenarios. Scenarios 1-12 are symmetric,
meaning the MTDs are on the same diagonal. Scenarios 1A-9A have asymmetric
toxicity matrices. Scenario 14 has all dose levels overly toxic. Scenario 15
has all dose levels below the TTR. Some scenarios have specific design
intentions indicated here by the number of correct orderings (CO) and the
number of overly toxic doses (OT). The CO number indicates the number of model
choices that would have a true ordering for each corresponding setting. The OT
number indicates the number of doses that have probability of toxicity greater
than $110\%$ of the TTR.

SM.2 CALIBRATION

The skeleton is optimised independently for POCRM and BMA-POCRM to ensure a
fair comparison. This is done by evaluating performance under scenarios $2$,
$6$, $5$A and $13$, which characterise a diverse set of scenarios. Scenario 2
and 6 have 1 and 2 correct orderings, respectively. Scenario 5A has 1 correct
ordering and an asymmetric MTD layout. Finally, scenario 13 is where all but
one dose level are overly toxic. The geometric mean across simulated trials is
considered as it more accurately accounts for poor performance relative to the
arithmetic mean. For each scenario, $5,000$ trials are simulated for which the
results are aggregated by taking the geometric mean of the PCS across all
$20,000$ trials. Trial simulations have a total of $60$ patients with $60$
cohorts of size $1$. The skeleton with the highest geometric mean PCS and
lowest standard deviation of PCS are selected for the complete simulation
study.

Two approaches for generating the underlying probability skeleton are
considered. First, the indifference interval protocol defined in Lee and
Cheung (2009) and implemented by the `dfcrm` `R` package (Cheung, 2019). This
approach has two parameters, the prior MTD and the $\delta$ parameter, which
controls the distance between the prior probability of toxicity for adjacent
doses in the skeleton. Grid optimisation across these parameters is carried
out to select the optimal indifference interval skeleton.

A linear skeleton is also implemented, whereby the prior probability of
toxicity at the lowest dose-level is given, and the skeleton increases at a
fixed interval. Recursively, this is defined as,

$\widehat{R}^{(0)}(d_{i})=\begin{cases}p_{0}&\text{if $i=1$,}\\\
\widehat{R}^{(0)}(d_{i-1})+\eta&\text{otherwise,}\end{cases}$

where $\widehat{R}^{(0)}(d_{i})$ is the prior probability of toxicity for dose
$i$, $p_{0}$ is the prior probability of toxicity for the lowest dose, or
_start probability_ and $\eta$ is the fixed increment, or _spacing_. Again,
grid optimisation is carried out to find the best linear skeleton.

The calibration results for POCRM under an indifference interval are presented
in Figure S1. These results show that for many prior MTDs, selecting delta as
0.03 leads to a misspecified skeleton. The same trend was empirically observed
for greater values of delta. In particular, when delta is set to 0.04 a
diverging marginal likelihood interval is observed for all settings of the
prior MTD in at least one scenario. Since the indifference interval protocol
begins at the prior MTD and obtains each adjacent skeleton value sequentially,
if $\delta$ is too high, the protocol will reach the bounds of 0 and 1
quickly, resulting in several adjacent values which are arbitrarily close to 0
and 1. For example, selecting $d_{15}$ as the prior MTD and $\delta=0.04$ the
following skeleton is obtained,

$\\{1.94e^{-12},4.17e^{-10},3.08e^{-08},\ldots,0.22,0.30,0.38\\}.$

A misspecified skeleton leads to a diverging marginal likelihood integral and
is observed both for BMA-POCRM and POCRM as they share a common approach to
computing the marginal likelihood in Equation (3). The skeletons proposed at
this level are frequently unintuitive and result in small differences between
dose levels at the boundaries near 0 and 1. As a result, $\delta\geq 0.03$ is
excluded for the indifference interval skeleton and
$\eta\;\text{(spacing)}\;=0.04,0.05,0.06$ for the linear skeleton due to the
missing values shown in Figures S1.A and S1.C, respectively.

Under the indifference interval skeleton, POCRM exhibits greatest PCS with a
prior MTD $d_{1}$ and $\delta=0.02$. Switching to a prior MTD $d_{2}$ leads to
a slightly lower PCS but also lower standard deviation in performance. This is
favourable as it implies greater consistency across scenarios. Assuming a
linear skeleton, the design performs optimally under $p_{0}=0.4$ and
$\eta=0.03$, however, this presents significantly greater SD than when
$p_{0}=0.01$ and $\eta=0.03$. Overall, the linear skeleton is much less
consistent across parametrisations when compared to the indifference interval
method for the POCRM.

Similarly, Figure S2 shows that under the indifference interval for BMA-POCRM,
optimal performance is observed according to mean PCS at a prior MTD $d_{2}$
and $\delta=0.02$. This setting also leads to the smallest SD across all valid
parametrisations of the indifference interval (i.e. excluding those where
$\delta>0.03$). This setting of the indifference interval also outperforms all
given parametrisations of the linear skeleton according to mean PCS. The
linear skeleton with $\eta=0.03$ and $p_{0}=0.01$ leads to lower SD, however,
this comes at the expense of a point reduction in PCS.

Moving forward, a matching skeleton is selected for both POCRM and BMA-POCRM.
The indifference interval skeleton with prior MTD $d_{2}$ and $\delta=0.02$ is
used due to the balance between aggregate performance and consistency across
both the methods included (i.e. comparatively low SD).

Figure S1: Results for prior calibration of skeleton parameters for the POCRM
with model selection. Panel A. Geometric mean of PCS for various values of the
indifference interval skeleton under the POCRM model. Panel B. Standard
deviation of PCS across scenarios selected for calibration when applying the
indifference interval skeleton under the POCRM model. Panel C. Geometric mean
of PCS for various values of the linear skeleton under the POCRM model. Panel
D. Standard deviation of PCS across scenarios selected for calibration when
applying the linear skeleton under the POCRM model. Missing values indicate at
least one diverging marginal likelihood integral over simulations. Figure S2:
Results for prior calibration of skeleton parameters for the BMA-POCRM. Panel
A. Geometric mean of PCS for various values of the indifference interval
skeleton under the BMA-POCRM model. Panel B. Standard deviation of PCS across
scenarios selected for calibration when applying the indifference interval
skeleton under the BMA-POCRM model. Panel C. Geometric mean of PCS for various
values of the linear skeleton under the BMA-POCRM model. Panel D. Standard
deviation of PCS across scenarios selected for calibration when applying the
linear skeleton under the BMA-POCRM model. Missing values indicate at least
one diverging marginal likelihood integral over simulations.

SM.3 VARIANCE CALIBRATION

The effect of altering the variance of the prior distribution of the $a$
parameter of the working model is also explored in Supplementary Materials. It
was found that the prior variance does not have any influence on the mean PCS
for neither of the two approaches, and hence the value of 1.34 is used for
further simulations.

Figure S3: Geometric mean of PCS across all scenarios. The optimal skeleton,
$\mathcal{S}$, is used here.

SM.4 FURTHER SIMULATION RESULTS

Figure S4: Distribution of changes in posterior model probabilities following
the induction of each cohort for simulated trials with cohort size 1.

SM.5 TRIAL DATA

Table S2: Real trial data generated by Gandhi et al. (2014) where each cell represents $y_{j}/n_{j}$ the number of observed DLTs and the number of patients assigned to each combination. | Temsirolimus (mg)
---|---
|  | 15 | 25 | 50 | 75
| 120 | 0/2 | 0/4 | 1/5 | 0/4
Neratinib (mg) | 160 | 1/4 | 1/4 | 0/5 | 3/6
| 200 | 0/4 | 1/8 | 1/2 | –
| 240 | 2/4 | – | – | –
Table S3: Key for doses allocated in the Gandhi et al. (2014) trial. | Temsirolimus (mg)
---|---
|  | 15 | 25 | 50 | 75
| 120 | $d_{1}$ | $d_{2}$ | $d_{3}$ | $d_{4}$
Neratinib (mg) | 160 | $d_{5}$ | $d_{6}$ | $d_{7}$ | $d_{8}$
| 200 | $d_{9}$ | $d_{10}$ | $d_{11}$ | –
| 240 | $d_{12}$ | – | – | –
Table S4: Dose allocations with POCRM in the motivating trial (Gandhi et al., 2014) where each cell represents $y_{j}/n_{j}$ the number of observed DLTs and the number of patients assigned to each combination. The starting combination is $d_{5}$. |  | Temsirolimus (mg)
---|---|---
|  | 15 | 25 | 50 | 75
Neratinib (mg) | 120 | 0/0 | 0/0 | 0/1 | 0/2
160 | 0/1 | 0/0 | 0/1 | 1/1
200 | 0/1 | 1/16 | 0/0 | –
240 | 15/30 | – | – | –
Table S5: Dose allocations with BMA-POCRM in the motivating trial (Gandhi et al., 2014) where each cell represents $y_{j}/n_{j}$ the number of observed DLTs and the number of patients assigned to each combination. The starting combination is $d_{5}$. |  | Temsirolimus (mg)
---|---|---
|  | 15 | 25 | 50 | 75
Neratinib (mg) | 120 | 0/0 | 0/0 | 0/0 | 0/2
160 | 0/1 | 0/0 | 0/0 | 1/1
200 | 0/1 | 1/17 | 1/2 | –
240 | 13/28 | – | – | –