crm12comb: Phase I/II adaptive design for drug combinations based on CRM design through R

TL;DR

The paper introduces crm12Comb, an R package that helps design adaptive Phase I/II trials for drug combinations using a Bayesian approach to improve cancer therapy development.

Contribution

The novelty lies in providing a flexible and user-friendly R package for adaptive drug combination trials using CRM methodology.

Findings

crm12Comb supports patient assignment and simulation studies for evaluating trial parameters.

The package handles toxicity and efficacy as binary outcomes with partial orderings for drug combinations.

It offers flexibility through user-specified parameters and various link functions with prior distributions.

Abstract

Adaptive designs for integrated phase I/II trials of drug combinations are increasingly utilized to speed up the drug development process and enhance drug efficacy, particularly in the realm of cancer therapy. The model-based Continual Reassessment Method (CRM) for dose-finding is widely used to leverage accumulated data in guiding patient allocation, drawing on principles from the Bayesian framework. In this paper, we present crm12Comb, an R package we developed to streamline the Phase I/II adaptive design process for drug combinations using the CRM approach. This package supports patients’ assignment guidance in a single trial based on existing data, as well as simulation studies for conducting extensive simulations with multiple trial parameters to evaluate operating characteristics and create visual representations. It accounts for toxicity and efficacy as binary outcomes, applying partial orderings to the dose-toxicity and dose-efficacy relationships for drug combinations. crm12Comb allows for a wide range of user-specified parameters, including maximum number of patients, cohort size, drug combinations, and a variety of link functions with prior distributions, offering flexibility to accommodate diverse clinical scenarios.

Figures

Flowchart of the general algorithm for one trial.

Patient enrolment with toxicity and efficacy outcomes for the first trial in simulation results 𝐬𝐢𝐦𝐑𝐞𝐬.

Patient allocations by dose combinations for the first trial in simulation results 𝐬𝐢𝐦𝐑𝐞𝐬.

Dose combinations identified as ODC among all simulations in 𝐬𝐢𝐦𝐑𝐞𝐬.

Probability of recommending target dose combinations as ODC among 1000 simulation trials.

Average proportion of patients allocated to target ODS(s) among 1000 simulation trials.

Equations

Full text

Introduction

In general, phase I dose-finding designs aim to find the maximum tolerated dose (MTD), which will be used as the recommended dose in the following phase II trials to assess the corresponding efficacy. The MTD is defined as the highest dose level with acceptable toxicity. In most oncology trials, the determination of the MTD involves considering dose-limiting toxicity (DLT). Traditionally, phase I and phase II trials are conducted in a separate and sequential manner. However, the recent trend has seen a rise in the integration of phase I and II trials, primarily due to the significant benefits this approach offers. These benefits include expediting the drug development timeline, minimizing the required number of participants, and ethically optimizing patient allocation to more efficacious dosages. Furthermore, the use of drug combinations in cancer therapy is a critical strategy to avert drug resistance and effectively target cancer cells. To this end, numerous clinical trial designs have been put forward to maximize these advantages.

Traditionally, phase I dose-finding designs have relied on algorithmic approaches, such as 3 + 3 design [1] and accelerated titration design [2], which utilize pre-specified rules to guide dose escalation/de-escalation. Despite their relatively modest performance, these algorithm-based designs are popular in practice [3] due to their simplicity and low DLT rates. To enhance operating effectiveness, model-based designs have been developed. These designs, such as continual reassessment method (CRM) [4] and escalation with overdose control (EWOC) [5], use accumulated data to estimate the dose-toxicity curve, offering more flexible and more accurate dose allocation for each patient cohort by utilizing all available information. However, these phase I dose-finding designs typically focus only on toxicity without considering efficacy and are in the context of single-drug escalations.

Meanwhile, various methods have been developed for dose-finding designs in drug combination trials. For example, linear programming has been employed to identify the optimal combined dose within the boundary of acceptable toxicity [6]. A plethora of Bayesian-based designs have been extensively applied to drug combinations [7–13]. Despite the diversity of these designs, a common limitation is their lack of consideration for efficacy alongside toxicity.

To assess both toxicity and efficacy simultaneously, the integrated phase I/II trials have been proposed, incorporating various definitions of toxicity and efficacy outcomes, such as methods that are based on ordinal outcome [14] and the correlated bivariate binary outcome [15] for single drugs. To accommodate for drug combinations, recent advancements have led to the proposal of various designs specifically tailored to incorporate drug combinations in integrated phase I/II trials, as demonstrated by recent works [16–20]. One practical approach to addressing drug combinations were to introduce the concept of an acceptable set [21,22]. In this approach, a fixed acceptable set is constructed based on phase I data to exclude those overly toxic doses in the phase II trial when assessing efficacy. Wages and Conaway [23] modified this concept to allow for continuous updates of acceptable drug combinations, aligning the efficacy assessment with the CRM design. This dynamic approach has been successfully implemented in real clinical trials [24,25]. However, their design assumed only the empiric link function with a standard normal prior distribution for estimating the toxicity and efficacy. In real-world scenarios, various combinations of link functions and prior distributions could be considered.

The CRM designs have gained widespread applications, bolstered by the development of various software tools that leverage their strong performances. The intricate nature of CRM methodologies underscores the substantial demand for software solutions that offer clinicians insights and simulations across different CRM design configurations. In the R programming environment, several packages have been developed to facilitate Phase I CRM designs, including CRM [26] which provides MTD and operating characteristics from simulations, dfcrm [27], which employs Bayesian inference without using Markov chain Monte Carlo (MCMC) sampling, and crmPack [28] and bcrm [29], which utilize MCMC sampling via JAGS/OpenBUGS/WinBUGS. Additionally, trialr [30] employs MCMC sampling via Stan [31,32], and pocrm [33] implements partial orderings for drug combinations. While these packages have been widely used, CRM/dfcrm/crmPack/bcrm can perform CRM design focusing only on toxicity with single drugs, pocrm utilizes partial ordering to drug combinations but only designed via toxicity, and trialr further includes efficacy but lacks considering the drug combinations. Currently, no R package can comprehensively cover both toxicity and efficacy for drug combinations.

In this paper, we address the existing gap by delving into a CRM-based phase I/II adaptive design for drug combinations that incorporate both toxicity and efficacy and include a wide range of link functions and prior distributions. We present an R package named crm12Comb, designed to facilitate the whole process and is publicly available via the Comprehensive R Archive Network (CRAN) at https://cran.r-project.org/web/packages/crm12Comb/. The crm12Comb package allows for simulation studies to be conducted before the trial commences, providing a general performance overview, or for estimating the next drug combination allocation based on current data.

Method

The initial methodology was proposed by Wages and Conaway [23], which considered the empiric link function with a standard normal prior distribution for the phase I/II adaptive design for drug combinations using CRM design. This approach involves constructing an acceptable set based on the estimated toxicity probabilities and determining the allocation of the next patient or patient cohort based on the estimated efficacy probabilities.

We extended the methodology by deriving various additional link functions and incorporating a wide range of prior distributions beyond the standard normal distribution, with customized parameter values. Additionally, we introduced the scale transformations tailored to the specific characteristics of each link function and prior distribution. (Fig 1) illustrates the general flowchart of this process.

Flowchart of the general algorithm for one trial.

In this section, we delved into the details of this algorithm. Specifically, the algorithm integrated diverse link functions such as hyperbolic tangent function, empiric function, one-parameter logistic function and two-parameter logistic function, making it compatible with a range of widely used prior distributions.

Construction from partial orderings

Partial ordering followed the principle that both dose-toxicity and dose-efficacy curves for a single drug should monotonically increase, which helped determine the possible orderings for the drug combinations with different dosages. More specifically, for any single drug, a higher dose level should correspond to higher DLT and efficacy rates, assuming the dose levels of other drugs remain constant. However, the number of potential orderings can grow exponentially with the addition of more dose levels for a single drug or more drugs in the combination, making it impractical to explore every possible combination. Six typical orderings were suggested to manage this complexity from practical designs [34].

To illustrate this, two drugs A and B with $[eqn]$ possible drug combinations $[eqn]$ were first considered, leading to $[eqn]$ complete toxicity orders and $[eqn]$ complete efficacy orders. We followed the example set by Wages and Conaway [23], which assumed three dose levels for each single drug, as shown in Table 1.

Table 1: 3×3 matrix for two drug combinations with different doses.

By holding the dose of one drug constant, the probabilities of toxicity and efficacy for the other dose increased monotonically across each row or column. For example, if we denoted the DLT probability at dose combination $[eqn]$ be $[eqn]$ , then the DLT probabilities followed $[eqn]$ when dose A was set at level 3. Consequently, six complete toxicity orderings for the 9 dose combinations were outlined below. Similarly, six complete efficacy orderings can be established using analogous reasoning.

Across rows: $[eqn]$ Up columns: $[eqn]$ Up diagonals: $[eqn]$ Down diagonals: $[eqn]$ Alternating down-up diagonals: $[eqn]$ Alternating up-down diagonals: $[eqn]$

Phase I/II CRM dose-finding design for drug combinations

A random variable $[eqn]$ denoted the drug combination for the $[eqn]$ th patient, where $[eqn]$ . Consider toxicity and efficacy as binary endpoints, so that for patient $[eqn]$ ,

$$[eqn]$$

For $[eqn]$ drug combinations, the pre-specified standardized unit, named toxicity skeleton, was used in the dose-toxicity curve to estimate the true toxicity probability $[eqn]$ . The toxicity skeleton had a length of $[eqn]$ , denoted by $[eqn]$ , where the $[eqn]$ th toxicity skeleton for toxicity ordering $[eqn]$ (denoted by $[eqn]$ ) was assigned to the drug combination with $[eqn]$ th toxicity probability arranged in ascending order for the specific toxicity ordering $[eqn]$ , across a total of $[eqn]$ toxicity orderings. It is worth noting that while the $[eqn]$ pre-specified toxicity skeleton remained the same, it was assigned to drug combinations with different orders corresponding to the toxicity orderings. These orderings had pre-determined prior probabilities $[eqn]$ , as outlined in the section about partial ordering, and satisfied the conditions that $[eqn]$ and $[eqn]$ . Similarly, the pre-specified efficacy skeleton for $[eqn]$ drug combinations, $[eqn]$ , had each efficacy ordering $[eqn]$ linked to a skeleton ordering, across $[eqn]$ efficacy orderings. These orderings came with prior probabilities $[eqn]$ , which satisfy the conditions that $[eqn]$ and $[eqn]$ . In case where there was no prior information on orderings, a discrete uniform distribution could be employed.

Starting the trial

The trial initiated by selecting the orderings of toxicity $[eqn]$ and efficacy $[eqn]$ with maximum prior probabilities in $[eqn]$ and $[eqn]$ . If multiple orderings shared the same maximum prior probability, one ordering was chosen at random. With the initial dose combinations $[eqn]$ chosen from $[eqn]$ , and the pre-specified toxicity skeleton $[eqn]$ , the acceptable set was then defined as

$$[eqn]$$

where $[eqn]$ was the pre-specified maximum acceptable toxicity rate (toxicity threshold). Subsequently, the randomization ratio was calculated as

$$[eqn]$$

in which $[eqn]$ represented the pre-specified efficacy skeleton for the efficacy ordering $[eqn]$ . The first patient cohort was randomized to combination $[eqn]$ with probability $[eqn]$ .

Toxicity estimation

Based on the partial orderings discussed before, the DLT probability for each ordering $[eqn]$ was expressed as: $[eqn]$ ,

where $[eqn]$ denoted the dose-toxicity model, often called the link function.

Upon enrolling $[eqn]$ patients, the collected toxicity data were represented by $[eqn]$ . For each ordering $[eqn]$ , the estimate of parameter $[eqn]$ , $[eqn]$ , was determined through Bayesian inference. In crm12Comb, the process for updating the likelihood function drew inspiration from the method used in the R package dfcrm [27] and was extended to include a variety of link functions with different prior distributions. Briefly, the prior distribution of $[eqn]$ was assumed to be $[eqn]$ , and $[eqn]$ represented the prior distribution for the ordering $[eqn]$ . The likelihood function was then formulated as:

$$[eqn]$$

From this, the posterior density for $[eqn]$ was defined as:

$$[eqn]$$

and posterior density for ordering $[eqn]$ became:

$$[eqn]$$

After each new patient or cohort of patients were enrolled, the probabilities of ordering $[eqn]$ were sequentially updated using the posterior density in Equation (3). The ordering $[eqn]$ with the largest posterior probability was chosen according to:

$$[eqn]$$

Then, the posterior probability of DLT was obtained by:

$$[eqn]$$

Based on this, the acceptable set were defined as $[eqn]$ . If at any point during the trial the acceptable set became empty, the next patient would be allocated to the lowest dose level, $[eqn]$ .

Efficacy estimation

The method for efficacy estimation was similar to the approach used for toxicity estimation. Let $[eqn]$ be the number of possible efficacy orderings, the efficacy probability for each ordering $[eqn]$ was then expressed as: $[eqn]$ .

Utilizing Bayesian inference, the likelihood function of order $[eqn]$ , based on the efficacy data collected from currently enrolled patients $[eqn]$ , was defined as:

$[eqn]$

Then the posterior probabilities for parameter $[eqn]$ and ordering $[eqn]$ were determined by:

$$[eqn]$$ $$[eqn]$$

From these, the estimate for $[eqn]$ , $[eqn]$ , was obtained. For subsequent patient or cohort of patients, the order $[eqn]$ with the highest posterior probability was selected, and the efficacy probability for each drug combination under ordering $[eqn]$ was estimated as:

$$[eqn]$$

Then $[eqn]$ was used to guide the allocation of new patient or patient cohort.

Patient allocation

To ensure each patient was allocated to the drug combination with the highest efficacy within the acceptable toxicity set, a two-phase allocation algorithm was implemented: an adaptive randomization phase for the initial stage with limited data, followed by a maximization phase when sufficient data are available.

The adaptive randomization phase aimed to prevent the trial from being stuck to the initially selected drug combinations, which allowed for exploration of a broader range of combinations. The randomization ratio $[eqn]$ was calculated as:

$$[eqn]$$

The next patient or cohort of patients was then randomized to combination $[eqn]$ with probability $[eqn]$ , for $[eqn]$ . This randomization phase was applied to a pre-specified subset of $[eqn]$ patients.

Upon exceeding the $[eqn]$ patient threshold, the trial transitioned to the maximization phase, where patients were allocated to the combination $[eqn]$ that maximized estimated efficacy probability, determined by:

$$[eqn]$$

The trial stopped once the total number of enrolled patients reached the pre-determined maximum of $[eqn]$ , or when pre-defined stopping rules were met. At this point, the optimal dose combination (ODC) was defined as $[eqn]$ after the maximum sample size $[eqn]$ patients have been enrolled. If the trial stopped earlier, no ODC was identified. We will introduce the details about early stopping rules in the following subsection.

Early stopping rules

The trial incorporated stopping rules for both safety concerns and lack of efficacy (futility). At any stage of the trial, the exact binomial 95% confidence interval for toxicity at the least toxic combination $[eqn]$ was computed as $[eqn]$ . If the lower bound $[eqn]$ exceeded the maximum acceptable toxicity rate $[eqn]$ , the trial stopped for safety, and no drug combination was designated as ODC.

Furthermore, the trial could also be stopped for futility if there were evidences suggesting no new dose combination had better effect than the existing ones in terms of efficacy. Upon enrolling $[eqn]$ patients, the exact binomial 95% CI for efficacy at current combination $[eqn]$ was computed as $[eqn]$ . If the upper bound $[eqn]$ fell below the pre-specified futility threshold $[eqn]$ , the trial stopped for futility, and no drug combination was designated as ODC. It is important to note that the futility-based stopping rule was only applicable in the maximization phase, which started after enrolling $[eqn]$ patients, since the upper bound was derived from the maximum estimated efficacy probability during this phase.

The implementation of link functions

The algorithm proposed by Wages and Conaway [23] in prior section utilized an empiric link function with a standard normal prior distribution to represent both the true probabilities of toxicity (under ordering $[eqn]$ ) and efficacy (under ordering $[eqn]$ ) as follows:

$$[eqn]$$

with $[eqn]$ and $[eqn]$ assuming standard normal prior distributions having expectation at $[eqn]$ , i.e., $[eqn]$ and $[eqn]$ . This approach, favoring the exponential form over a linear one, accounted for the presumed monotonic nature of the dose-toxicity and dose-efficacy curves.

Expanding on this framework, this section introduces a broader range of prior distributions for $[eqn]$ and $[eqn]$ , allowing for more versatile selections. The detailed derivation of two-parameter link functions is provided in the S1 File. Additionally, the corresponding scale transformation of the skeleton was also applied, following the methodology outlined in [4], to enhance the model’s adaptability. Since the above prior distribution of empiric link function followed $[eqn]$ , no scale transformation is needed. The details are given in Table 2.

Table 2: Link functions for toxicity and efficacy probabilities with scale transformations and prior distributions.

Using crm12Comb

The R package crm12Comb contains functions that can perform simulations for integrated Phase I and II adaptive design for drug combinations using the CRM dose-finding design. It offers features to generate skeletons and specifications for different link functions, enabling users to conduct comprehensive simulation studies with pre-specified information to assess operating characteristics, or to determine the appropriate dose combination for upcoming patient or patient cohort based on all accumulated data.

To begin, we need to first install the crm12Comb package, and then load it as shown below:

**: **

Following that, in the section below we will provide examples on how to use the crm12Comb to simulate multiple trials to obtain operating characteristics of phase I/II adaptive design for drug combinations using CRM. Additionally, we will demonstrate how to simulate a single trial for the allocation of the next patient or patient cohort based on the most recent cumulated data. Moreover, we will illustrate an example code to perform multiple simulations with different sets of inputs for trial settings.

Simulation of multiple trials

Using a $[eqn]$ dose combination for two drugs as an example, the levels of dose combinations can be obtained based on the function $[eqn]$ ,

**: **

Then, the pre-defined true probabilities of toxicity and efficacy related to the two-drug combinations can be input sequentially, corresponding to the data presented in Table 3. Other values need to be pre-specified include 1000 simulation trials, maximum sample size as 40, subset number of patients as 20, toxicity threshold $[eqn]$ , efficacy threshold is 30%, cohort size as $[eqn]$ patient, correlation between toxicity and efficacy binary endpoints as $[eqn]$ indicating independence between toxicity and efficacy, and empiric link function with normal prior distribution (mean $[eqn]$ and standard deviation $[eqn]$ ) for both toxicity and efficacy probabilities.

Table 3: Example of true (toxicity, efficacy) probabilities for a 4 ×4 combinations.

**: **

The toxicity and efficacy skeletons can be constructed by either using prior information from historical data or by employing the algorithm proposed by Lee and Cheung [35]. The $[eqn]$ function from the R package dfcrm [27], which implemented the algorithm of the empiric model and one-parameter logistic model, assuming a normal prior distribution of the parameter with a mean of zero. We adapted the $[eqn]$ function in dfcrm to build the $[eqn]$ function in our crm12Comb, making it compatible with the link functions and prior distributions in the methods section. This adaptation allows for user-defined parameters for various distributions. The toxicity and efficacy skeletons are generated using the $[eqn]$ function, and we followed the same default settings of $[eqn]$ function [27] in this example.

**: **

The steps above are the preparations for the simulations. Next, we can run simulations using the $[eqn]$ function, which outputs $[eqn]$ operating characteristics. These include:

(1)the probability of recommending safe/ineffective combinations as ODC ( $[eqn]$ in R output),(2)the probability of recommending target combinations as ODC ( $[eqn]$ in R output),(3)the probability of recommending toxic combinations as ODC ( $[eqn]$ in R output),(4)the average number of patients enrolled ( $[eqn]$ in R output),(5)the average proportion of patients allocated to target ODS(s) ( $[eqn]$ in R output),(6)the probability of simulation trials stopped early for safety ( $[eqn]$ in R output),(7)the probability of simulation trials stopped early for futility ( $[eqn]$ in R output),(8)the average observed DLT rate ( $[eqn]$ in R output),(9)the average observed response rate ( $[eqn]$ in R output).

Besides the 9 operating characteristics, we also store the data of each trial into a list through the number of simulations ( $[eqn]$ in R output). Since the DLT and efficacy data are randomly generated and the adaptive randomization phase involves randomly selecting the next allocation combination based on the randomization ratio in Equation (4), results may vary slightly across different selections of random seeds. Therefore, we recommend conducting a minimum of 1000 simulations to ensure the reliability of the results.

The $[eqn]$ function includes the following arguments: $[eqn]$ for number of simulation trials, $[eqn]$ for number of maximum enrolled patients for each trial, $[eqn]$ for pre-defined true toxicity and efficacy probabilities for each drug combinations, $[eqn]$ for drug combinations, $[eqn]$ for format of drug combinations, $[eqn]$ for number of subset patients, $[eqn]$ for pre-specified toxicity skeleton, $[eqn]$ for pre-specified efficacy skeleton, $[eqn]$ for DLT threshold to define acceptable set with default value $[eqn]$ , $[eqn]$ for efficacy threshold to define target combinations with default value $[eqn]$ , $[eqn]$ for prior probabilities of toxicity orderings (sum is 1), $[eqn]$ for prior probabilities of efficacy orderings (sum is 1), $[eqn]$ for number of patients in each cohort, $[eqn]$ for association parameter for efficacy and toxicity (0 means independent), $[eqn]$ for safety threshold for early stopping with default value $[eqn]$ , $[eqn]$ for futility threshold for early stopping with default value $[eqn]$ , $[eqn]$ for the model of link function with default model $[eqn]$ , $[eqn]$ for the prior distribution of parameters in link function with default $[eqn]$ , and other arguments are to specify the values of each parameter based on the prior distribution. For detailed description, run $[eqn]$ to check the documentation.

**: **

In addition, the dataset, including columns for dose combination, DLT outcome, and efficacy outcome for each simulation trial, is stored in a list, which can be viewed through $[eqn]$ . In this example of 1000 simulations, the $[eqn]$ is a list containing 1000 elements, each representing patient data enrolled sequentially. This command can also allow for the examination of detailed patient allocations, where efficacy outcome is denoted by ORR. For example, if we output the first 5 enrolled patients ([1:5,]) from the first simulated trial ( $[eqn]$ ), the first patient is allocated to drug combination #5 without having any DLT or ORR outcomes based on the first row in the output below.

**: **

For visualization, three types of plots can be generated. To show the sequential enrolment of patients for a single trial, the $[eqn]$ function can be utilized. For example, to visualize the patient entry of the first simulated trial from the results $[eqn]$ , we can use the syntax below to generate Fig 2, where the x-axis denotes the sequential order of patient enrolments and the y-axis represents the assigned dose combination for each patient. Patients experiencing DLT outcomes are indicated by a red-left half dot, those with efficacy outcomes by a green-right half dot, and patients without such outcomes are represented by grey dots.

Patient enrolment with toxicity and efficacy outcomes for the first trial in simulation results 𝐬𝐢𝐦𝐑𝐞𝐬.

outcomes are represented by grey dots.

**: **

The second plot, a histogram, displays the distribution of patients across each dose combination in a single simulated trial, which can be generated using the syntax $[eqn]$ funtcion. To visualize the patient allocations of the first simulated trial from the results $[eqn]$ , we can use the syntax below to generate Fig 3 with specific number of patients allocated to each dose combination indicated at the top of each bar.

Patient allocations by dose combinations for the first trial in simulation results 𝐬𝐢𝐦𝐑𝐞𝐬.

**: **

The third plot is also a histogram that represents the number of trials in which each dose combination was identified as ODC. The plot can be created using the syntax $[eqn]$ function. For the 1000 simulation trials contained in $[eqn]$ , we can use the syntax below to generate Fig 4. This figure includes the specific counts of trials where each dose combination was deemed the ODC, indicated at the top of each bar. It is worth noting that among 1000 simulation trials, 4 trials with no ODC identified has been excluded from the figure because of stopping early for futility ( $[eqn]$ ).

Dose combinations identified as ODC among all simulations in 𝐬𝐢𝐦𝐑𝐞𝐬.

**: **

The outcomes related to operating characteristics may vary slightly by choosing different random seeds for the random nature of the following factors: (1) the toxicity and efficacy outcomes are randomly generated based on a bivariate binary correlation (which was set as 0 in the example above, $[eqn]$ ), and (2) when multiple dose combinations have the same efficacy skeleton (when starting the trial and selecting dose combination for first patient allocation) and have the same maximum estimated efficacy probability (during the maximization phase), one dose combination will be chosen at random. Therefore, to ensure more reliable results while maintaining computational efficiency, conducting 1000 simulation runs is recommended.

Simulation of a single trial patient allocation

Helper functions (such as $[eqn]$ , $[eqn]$ , and $[eqn]$ ) have been created to streamline the whole simulation process that is included in the $[eqn]$ function, as described in this section. These functions can also be used to perform a single simulation to provide the allocation of the next patient based on the existing data.

An example dataset, representing six patients already enrolled using a 4x4 dose combination matrix, is provided below.

**: **

The $[eqn]$ function, as described in the section about partial ordering, can be implemented to get the possible orderings for dose combinations when there is no prior information available for the combined drugs. With prior knowledge, users can specify customized orderings by aggregating them into a list for the argument $[eqn]$ and set the argument $[eqn]$ as $[eqn]$ of the $[eqn]$ function.

**: **

From the orderings generated by the $[eqn]$ function, we can obtain the toxicity and efficacy orderings by modifying the orderings of the pre-specified skeletons, which were generated by the $[eqn]$ function described before. To condense the output, we only keep 3 digits in the example below (function with $[eqn]$ ). However, in the actual calculation, we will use the original digits of each skeleton probability.

**: **

Based on these established toxicity skeletons and orderings, the $[eqn]$ function can be used to estimate toxicity probabilities for each dose combination, given current data. The function returns a list that contains dose combinations within the acceptable set and the posterior density of toxicity orderings, which are instrumental for subsequent toxicity estimations.

**: **

Following the determination of the acceptable set of dose combinations from the toxicity estimation (as indicated by the results below $[eqn]$ ), the $[eqn]$ function can be used to estimate efficacy probabilities. This allows us to identify the most suitable dose combination for allocating the next patient or patient cohort.

**: **

The next patient or cohort of patients will be allocated to dose combination #1 and posterior density of efficacy orderings will be used for the next efficacy estimations.

Examples of multiple simulations

To process multiple lists of input values for various parameters (such as $[eqn]$ , $[eqn]$ , $[eqn]$ , toxicity/efficacy skeletons, link functions, prior distributions, etc.), a $[eqn]$ loop can be employed with the $[eqn]$ function to obtain multiple sets of results. We adopted the scenario settings, including six scenarios with pre-defined true toxicity and true efficacy probabilities [23], as shown in Table 3. Here, for illustration, we further compare the results of three different link functions: empiric model with normal prior, hyperbolic tangent model with exponential prior, and one-parameter logistic model with normal and gamma priors. Besides, we also considered various combinations of two sets of toxicity and efficacy skeletons, the maximum number of patients per simulation trial ( $[eqn]$ ), correlation between toxicity and efficacy ( $[eqn]$ as illustrated by Thall and Cook [36]), and the number of subset patients for choosing randomization or maximization phase ( $[eqn]$ ).

Using the empiric link function with a normal prior as an example, the code below demonstrates a method to derive results for a list of input values. First, we need to input the six scenarios in the form of $[eqn]$ matrices. The examples of all six scenarios are shown below, with i th row denotes the toxicity and efficacy probabilities for the i th drug combination. $[eqn]$ shows the corresponding drug combinations linked to the positions in Table 4.

Table 4: Six scenarios of true (toxicity, efficacy) probabilities [23].

**: **

Based on the above scenarios, we constructed a total of 54 conditions, derived from 2 pre-specified toxicity and efficacy skeletons, 3 maximum sample sizes, 3 numbers of patients for determination of randomization phase, and 3 correlation parameters for efficacy and toxicity.

**: **

Then, we specified two sets of toxicity and efficacy skeletons so that we could plug the skeletons into the $[eqn]$ function to store the 54 conditions into a list named $[eqn]$ .

**: **

Next, we built an output dataset with detailed names for 9 operating characteristics.

**: **

The last step iterated through the 6 scenarios and 54 conditions, each with 1000 simulations. Results for each iteration were stored into an output dataset.

**: **

The output data were saved in the crm12Comb, named as $[eqn]$ , comprising 1296 rows (6 scenarios, 4 different combinations of link functions and prior distributions, 2 sets of skeletons, 3 maximum number of patients, 3 toxicity and efficacy correlations, and 3 subset number of patients) that represent each condition by 1000 simulations, along with 15 columns that contain 9 operating characteristics, other 6 columns including $[eqn]$ , $[eqn]$ for link function and prior distribution, $[eqn]$ for maximum number of patients, $[eqn]$ , $[eqn]$ for subset number of patients, and $[eqn]$ for toxicity and efficacy correlation to separate each condition. Various plots can be generated by the $[eqn]$ function to visually explore the relationships between operating characteristics and the maximum number of patients, number of patients for determination of the randomization phase, set of skeletons, and correlation between toxicity and efficacy binary outcomes by different scenarios and link functions.

Fig 5 and Fig 6 display sample plots derived from all simulations. These plots showed the probability of recommending target dose combinations as ODC and the average proportion of patients allocated to target ODC(s) among 1000 simulation trials versus the maximum number of patients, for all six scenarios, while fixing the other three variables. An observation from these plots is that no link function consistently outperforms the other three link functions, while $[eqn]$ link function exhibited worse performance in Scenarios 2–5.

Probability of recommending target dose combinations as ODC among 1000 simulation trials.

Average proportion of patients allocated to target ODS(s) among 1000 simulation trials.

**: **

**: **

Discussion

In this paper, we introduced crm12Comb, an accessible $[eqn]$ package tailored for conducting simulation studies on phase I/II adaptive design for drug combinations within a Bayesian CRM framework. crm12Comb offers the capability to perform not only a single trial simulation of the next patient allocation given current data, but also multiple trials to obtain a comprehensive preview of operating characteristics even before a trial starts. We also demonstrated how to execute multiple simulation studies with a range of input parameters.

Building upon the algorithm using the one-parameter empiric link function with standard normal prior distribution, as suggested in [23], crm12Comb incorporates three commonly used link functions in clinical research, enhancing the package’s versatility. Users can choose from various link functions with self-determined values of parameters for different prior distributions, including the one-parameter empiric model with normal or gamma priors, one-parameter hyperbolic tangent model with exponential priors, one-parameter logistic model with normal or gamma priors, and two-parameter logistic model with normal or gamma priors. This selection allows for the customization of prior distributions based on historical data, thereby improving the package’s utility for diverse clinical scenarios.

In addition to versatile link functions with priors, crm12Comb enables users to tailor various aspects of their phase I/II drug combination trial based on CRM design, including the maximum sample size, number of dose combinations, pre-defined true toxicity and efficacy probabilities for dose combinations, cohort size, the number of patients for determination of randomization phase, pre-specified toxicity and efficacy skeletons, correlation between toxicity and efficacy binary outcomes, DLT and efficacy thresholds, and thresholds for early stopping.

Furthermore, crm12Comb facilitates the creation of visualization, including plots for single-trial patient sequential enrolment with toxicity/efficacy outcomes, single-trial patient allocation by dose combinations, and the overall dose combination identified as ODC(s) among all simulations. Sample data and plots comparing different link functions are also provided within the package to serve as practical references.

A current limitation of the package is the lack of available data from real-world phase I/II drug combination trials, particularly patient-level toxicity, efficacy, and allocation information. Our review identified only two trials aligning with our package’s focus; however, neither provided the detailed patient-level outcomes or trial design specifics needed to reproduce the trial process. To address this gap, our package crm12Comb provides flexible simulation inputs that allow users to mimic realistic trial settings and explore practical scenarios. As relevant patient-level data become available, future development will focus on integrating real trial examples, to further demonstrate and evaluate the package’s practical utility.

Supporting information

S1 FileDetailed toxicity and efficacy estimation of two-parameter logistic link function.(DOCX)