Flexible use of copula‐type model for dose‐finding in drug combination clinical trials

Abstract Identification of the maximum tolerated dose combination (MTDC) of cancer drugs is an important objective in phase I oncology trials. Numerous dose‐finding designs for drug combination have been proposed over the years. Copula‐type models exhibit distinctive advantages in this task over other models used in existing competitive designs. For example, their application enables the consideration of dose‐limiting toxicities attributable to one of two agents. However, if a particular combination therapy demonstrates extremely synergistic toxicity, copula‐type models are liable to induce biases in toxicity probability estimators due to the associated Fréchet–Hoeffding bounds. Consequently, the dose‐finding performance may be worse than those of other competitive designs. The objective of this study is to improve the performance of dose‐finding designs based on copula‐type models while maintaining their advantageous properties. We propose an extension of the parameter space of the interaction term in copula‐type models. This releases the Fréchet–Hoeffding bounds, making the estimation of toxicity probabilities more flexible. Numerical examples in various scenarios demonstrate that the performance (e.g., the percentage of correct MTDC selection) of the proposed method is better than those exhibited by existing copula‐type models and comparable with those of other competitive designs, irrespective of the existence of extreme synergistic toxicity. The results obtained in this study could motivate the real‐world application of the proposed method in cases requiring the utilization of the properties of copula‐type models.

2019). MTDC is defined as the dose combination that has a dose-limiting toxicity (DLT) probability closest to the target toxicity probability.
Various types of dose-finding designs have been proposed to identify MTDC accurately. For two-drug combinations, Yin and Yuan (2009a) proposed a Bayesian dose-finding design using the following copula-type model based on the Clayton copula.
In the remainder of this study, ( , ) denotes the combination of the th dose ( = 1, 2, … , ) of drug and the th dose ( = 1, 2, … , ) of drug . Further, represents the joint toxicity probability corresponding to dose combination ( , ). and are the prespecified toxicity probabilities corresponding to the th dose of drug and the th dose of drug , respectively, > 0, > 0, and > 0 are the model parameters, and represents the effect of the drug-drug interaction. Yuan (2009a, 2009b) proposed two other copula-type models based on the Gumbel-Murtaugh and Gumbel-Hougaard copulas.
The copula-type model has distinctive properties that have been applied to dose-finding designs considering specific situations. For example, Jimenez et al. (2019) proposed a design based on a copula-type model for a trial in which a clinician is able to identify certain toxicities attributable to one of the two constituent agents. In addition, Wheeler et al. (2017) proposed a design applying a copula-type model to a trial in which one of the two constituent agents is administered at the beginning of the treatment cycle and the other is administered at a later instant during the cycle. This allows DLTs observed before the administration of the second drug to be attributed to the first drug. The aforementioned examples instantiate the enhanced versatility of copula-type models compared with models used in other competitive designs. However, Gasparini et al. (2010) argued that copula-type models are not suitable for the estimation of joint toxicity probabilities in cases where extreme synergistic or extreme antagonistic behavior is expected, owing to the Fréchet-Hoeffding (F-H) bounds; max( , ) ≤ ≤ min( + , 1). In response to Gasparini's comment, Yin and Yuan explained the nonexistence of lower and upper bounds in their model by arguing that the upper and lower bounds in their copula-type model is given by max( , ) ≤ ≤ min( + , 1), which allows to take any value in the range between 0 and 1 (Yin and Yuan, 2010). However, this is valid only during the estimation of the joint toxicity probability corresponding to a single dose combination. As an illustrative example, let us consider the case in which the joint toxicity probabilities, 12 , 21 , and 22 , are equal to 0.05, 0.05, and 0.30, respectively. It is evident from Equation (2) that the maxima of 12 and 21 are 1 + 2 and 2 + 1 , respectively. They can be transformed into 2 = 21 − 1 and 2 = 12 − 1 , respectively. Similarly, as the maximum of 22 is 2 + 2 , 22 = 21 − 1 + 12 − 1 holds. It can be readily verified that, when 1 = 1 = 0, 22 attains the maximum value. If 21 and 12 were estimated without the bias (i.e., 0.05), the bias of the estimator for 22 would be 0.20. On the contrary, the proximity of the estimated value of 22 to 0.30 is proportional to the bias of the estimator of 21 and 12 . Therefore, irrespective of the values of and , a copula-type model is restricted by the F-H bounds and the inequality,̂, <̂− 1, +̂, −1 , always holds in such models.
The degraded performance of dose-finding designs based on copula-type models compared with other competitive designs under scenarios involving extreme synergistic toxicity have also been verified in several studies (Wages et al., 2011;Riviere et al., 2014;Hirakawa et al., 2015;Mander and Sweeting, 2015;Yin, 2016, 2017). Certain phase I trial reports have suggested that combination therapy demonstrates extreme synergistic toxicity. For example, 62 patients were evaluated in the phase Ib doseescalation study of the combination of the oral Pan-PI3K inhibitor, buparlisib, and the oral MEK1/2 inhibitor, trametinib (Bedard et al., 2015). Let ( , ) denote the dose levels of combinations of buparlisib and trametinib. The toxicity probabilities corresponding to dose combinations (1, 3), (1, 4), (2, 3), and (2, 4) were estimated to be 0.03, 0.06, 0.15, and 0.33, respectively, using a logistic regression model (Lam et al., 2019). This suggests the existence of an extreme synergistic toxicity. In addition, extreme synergistic toxicity has been observed in the phase I study of the combination of neratinib and temsirolimus (Gandhi et al., 2014). Therefore, dose-finding design based on copula-type models may not be suitable in several practical clinical trials.
The objective of this study is to improve the performances of dose-finding designs based on copula-type models while maintaining their advantageous properties. We propose an extension of the parameter space of in copulatype models under the assumption that any function that takes values between 0 and 1 can be used to estimate . The proposed extension releases the F-H bounds, making the estimation of toxicity probabilities more flexible. We conduct numerical evaluations to verify the performance of the proposed method and compare it with those of other competitive designs.
The remainder of this paper is structured as follows. In Section 2, we describe the characteristics of a copula-type model and introduce certain copula-type models based on well-known copulas. In addition, we formulate the proposed method. Section 3 details the simulation setting along with the dose-escalation algorithm and MTDC-determination procedure of each dose-finding design. In Section 4, we explain the results of the numerical studies under various scenarios. Section 5 presents some concluding remarks and considerations.

PROBABILITY MODELS
2.1 Copula-type models Yuan (2009a, 2009b) introduced the categorization of the toxicity probability, , of the dose combination, ( , ), into the following four classes of joint probabilities.
(1) 11 , associated with the DLT attributed to both drugs; (2) 10 , associated with the DLT attributed to drug , but not to drug ; (3) 01 , associated with the DLT attributed to drug , but not to drug ; and (4) 00 , which is the probability of not observing a DLT (i.e., the DLT attributed to drug A or drug B or both is not observed).
Owing to the difficulty of explicit attribution of an observed DLT to drug , drug , or both, the authors focused on 00 and proposed the estimation of based on the transformation, = 1 − 00 . Let (⋅, ⋅) denote a copula connecting univariate marginal distribution functions. The joint toxicity probability, , of the dose combination, ( , ), based on a copula-type model is given by Here, ( ) = and ( ) = are the marginal cumulative distribution functions of the observed toxicity. and represent the prespecified toxicity probabilities corresponding to the th dose of drug and the th dose of drug , respectively, where 1 < 2 < … < and 1 < 2 < … < . The Clayton copula-type model given by Equation (1) is derived by the Clayton copula of ( , ) = ( − + − − 1) −1∕ . Other copulas may also be applied in a similar fashion. Yin and Lin (2015) compared the performances of dose-finding designs based on six types of copula-type models and demonstrated that there is no substantial difference between them. Table 1 illustrates the copula-type models based on wellknown one-parameter copulas (two-dimensional copulas) along with the respective parameter spaces of and the ranges of under fixed values of = = 0.15. As is evident from the table, the parameter space of and the range of are unique for each model. The Clayton and Frank model attains the upper and lower F-H bounds given by Equation (2). However, even in the case of the Clayton and Frank model, the range of is not sufficiently wide. This is because, in cases involving extreme synergistic toxicity, is expected to be greater than min( + , 1). Thus, in such cases, dose-finding designs based on copula-type models are not effective. In contrast, copula-type models exhibit the following properties: (1) if = 0 and = 0, then = 0; (2) if = 0, then = , and if = 0, then = ; (3) if = 1 or = 1, then = 1; (4) there is a unique value of for which Property (2) is notable because it asserts that is equal to the toxicity probability of one drug if the toxicity probability of the other drug within the combination is known to be zero. As discussed in Section 1, several dose-finding designs have been proposed based on this property and the concepts of copula-type models described at the beginning of this section (Wheeler et al., 2017;Jimenez et al., 2019). Property (4) implies that any copula-type model is capable of expressing with no interaction. In addition to a copula-type model, the models proposed by Thall et al. (2003) and Wang and Ivanova (2005) also satisfy the properties of (1) to (3).
In general, the parameter space of corresponding to each model is derived to ensure that the two-dimensional joint distribution is an actual probability distribution. We outline the derivation of the parameter space of the Ali-Mikhail-Haq (AMH) model (Ali et al., 1978). The two-dimensional joint distribution function (JDF), ( , ), connected by the AMH copula, which corresponds to the joint toxicity probability, 11 , is given by .  In addition, the corresponding joint toxicity probability density function is given by It is readily apparent that 0 ≤ ( ), ( ) ≤ 1, ( , 1) = ( ), (1, ) = ( ), ( , 0) = (0, ) = 0, and (1, 1) = 1. To ensure that ( , ) ≥ 0 for all and , is required to satisfy < 1. In addition, to ensure that the joint toxicity probability density function, ( , ) ≥ 0 for all and , is required to satisfy −1 ≤ < 1. Therefore, by combining the constraints introduced by the desired properties of ( , ) and ( , ), the parameter space of is identified to be −1 ≤ < 1.

Flexibility of application of copula-type models
In copula-type models, the likelihood is simply a product of the binomial probability mass function, which is given by where denotes the data collected until a given instant, and and denote the number of patients who received the drug combination, ( , ), and experienced DLTs corresponding to the drug combination, ( , ), respectively. We assume the DLT outcomes of all patients to be independent. In a copula-type model, is represented by a two-dimensional JDF connecting two univariate marginal distribution functions. The parameter space of in a copula-type model is determined such that it meets the two conditions mentioned in Section 2.1. By relaxing one of these conditions, 2 ( , )/ ≥ 0, we propose an extension of the parameter space of within the range satisfying the restriction that should lie within 0 and 1, thereby corresponding to the condition 0 ≤ ( , ) ≤ 1. This extension of leads to an increase in the F-H upper bound (henceforth denoted by upper ). The proposed method with this extended parameter space utilizes the regression model using the same formula (function) as the original copula-type model to estimate . The function in the original method corresponds to a JDF derived using a copula, but the function in the proposed method does not completely satisfy the condition of a JDF because the corresponding density function, ( , ) = 2 ( , )∕ , does not always provide a positive value in the extended range of . Here, ( , ) = (1 − , 1 − ), where represents a copula. The function in the proposed method is not a JDF but a mere function that takes values between 0 and 1; however, the proposed method inherits the advantageous properties of copula-type models. Web Appendix A records the properties of a copula that are preserved by the proposed method, along with the proofs of these statements, and establishes that the proposed method is equipped with the minimal properties of a copula that are required for the modeling of . Here, we would like to emphasize that any function that takes values between 0 and 1 can be used for modeling the relationship between the DLTs and dose combinations because is only a parameter of the binomial probability mass function in the likelihood of Equation (3). When 0 ≤ ≤ 1 is preserved, the resultant likelihood is proper. Table 2 lists the extended parameter spaces of and the corresponding ranges of . We selected the AMH, GUMH (Gumbel-Hougaard, AMH: Ali-Mikhail-Haq), and Joe models because their parameter spaces of can be extended easily irrespective of and . For instance, in the case of AMH, the parameter space of is (−∞, 1). However, in the case of FGM (Farlie-Gumbel-Morgenstern), it is required to be above −( ) −1 to ensure that lies between 0 and 1. It is evident from Table 2 that the extension of the parameter space of enables upper to be increased to 1. The model proposed by Wang and Ivanova (2005) was not developed as a copula-type model, but Gasparini (2013) used his concept of dose-free parameterization to prove that it is equivalent to a copula-type model if the original parameter space of 0 < is restricted to 0 < < 1. If the model is not regarded as a copula-type model, it seems natural not to restrict the parameter space, as done by Wang and Ivanova. It appears that their model with its original parameter space is one that can be constructed based on the proposed method. The stopping rule used in original copula-type models can be implemented in the proposed method as well. The trial can be terminated at an early stage for safety if even the lowest dose combination is overly toxic, as indicated by Pr( 11 > | D) > , where denotes a threshold value close to 1.

Common settings
We conducted simulation studies under the 18 scenarios listed in Table 3 to investigate the operating characteristics of a dose-finding design based on an existing copulatype model (COPULA), the proposed method based on the extension of the parameter space of (COPULA-E), and four existing dose-finding designs-two-dimensional model-based design (TDM) (Wang and Ivanova, 2005), logistic model-based design (LOGI) (Riviere et al., 2014), partial ordering CRM design (POCRM) (Wages et al., 2011), and two-dimensional Bayesian optimal interval design (BOIN) (Lin and Yin, 2017). We utilized three copula-type models, namely, AMH, GUMH, and Joe in both COPULA (the designs are denoted by AMH-O, GUMH-O, and Joe-O) and COPULA-E (the designs are denoted by AMH-E, GUMH-E, and Joe-E). The 18 scenarios comprised the 16 scenarios considered by Hirakawa et al. (2015) and 2 additional scenarios that do not involve synergistic toxicity corresponding to any of the dose combinations. Of the 18 scenarios, we considered scenarios 4, 5, 6, 7, 8, 14, 15, and 16 to exhibit extreme synergistic toxicity around true MTDCs. The scenario satisfying the condition, , was identified to be the scenario exhibiting extreme synergistic toxicity, where Pr( , ) represents the toxicity probability of the dose combination, ( , ), and the dose combination, ( , ), is an MTDC. In all of the scenarios considered, the target toxicity probability, , was set to 0.3, and dose combinations with toxicity probabilities greater than or equal to 0.35 were defined to be overdose combinations (OCs). The total sample size was taken to be 30 with a cohort size of 3. A total of 1000 trials were simulated corresponding to each scenario. We did not employ any stopping rules, that is, the MTDC is determined after attaining the maximum sample size.
We evaluated the operating characteristics in terms of six outcome metrics. The first three outcome metrics are the percentage of correct MTDC selections, percentage of OC selections, and accuracy index (AI) proposed by Cheung (2011); these represent the efficiency and accuracy of a dose-finding design. The remaining three outcome metrics are the overall percentage of observed toxicity, percentage of patients allocated to MTDC, and percentage of patients allocated to OC; these represent the safety perspectives of a trial. AI is defined as follows.
where * indicates the true toxicity probability of the dose combination, ( , ), and denotes the probability of selecting the dose combination, ( , ), as the MTDC. The maximum value of AI is 1, and larger values correspond to higher dose-finding accuracies. In this study, we multiplied all values of AI by 100.
The model parameters for COPULA, COPULA-E, LOGI, and TDM are estimated via the MCMC using PROC MCMC of SAS 9.4 (SAS Institute, Cary, NC). We record 10 000 posterior samples after 1000 burn-in iterations. In addition, the number of chains and the thinning interval are taken to be 1 and 3, respectively.

Dose-escalation algorithm and MTDC determination
In all of the dose-finding designs considered, the first cohort is treated at the lowest dose combination (1, 1). Following is the description of the dose combination for the next cohort. In the case of BOIN, the dose-escalation algorithm proposed in the original article is applied. The doseescalation algorithms of other designs comprise two stages. During the first stage, the dose is escalated based on a prespecified rule without using statistical models. This is continued until one DLT is observed (alternatively, in the case of POCRM, until the maximum likelihood estimates are obtained). During the second stage, the dose is escalated or de-escalated based on the values predicted by the models. Although the dose-finding designs usually utilize distinct first stages, the first stage proposed by Wages et al. (2011) is employed for all designs in this study, to ensure homogeneity of conditions for a fair comparison. During this stage, the dose is sequentially escalated from the lowest dose group to the highest dose group. A dose group is defined to be a set of dose combinations lying on the same diagonal in the drug combination matrix. Two drug combinations lie on the same diagonal if the sums of their corresponding dose levels are equal. When multiple dose combinations exist in a dose group, one of them is sampled without replacement for the next cohort. This sampling is continued until all dose combinations in a dose group are selected.
During the second stage, the dose-escalation algorithm proposed in each original article is applied to the corresponding dose-finding design. With respect to the current dose combination, ( , ), the second stage of COPULA and COPULA-E can be described as follows.
• If Pr( < | ) ≤ and Pr( > | ) ≤ , we retain the current dose combination for the next cohort, where and represent the threshold for dose escalation and dose de-escalation, respectively, satisfying + > 1. We take = 0.75 and = 0.55 for COPULA and COPULA-E. LOGI employs the same algorithm as COP-ULA, and we use their default settings-= 0.85 and = 0.55-in its case.
The MTDC-determination procedure proposed in each original article is applied to the corresponding dosefinding design. COPULA and COPULA-E are taken to employ the determination method used in LOGI. In LOGI, the MTDC is selected from the dose combinations used to treat patients and the dose combination with the highest posterior probability, Pr( ∈ [ − , + ]| ), is taken to be the MTDC, where we set = 0.1.

COPULA versus COPULA-E (proposed method)
In this subsection, we compare the performances of the original method (COPULA: AMH-O, GUMH-O, and Joe-O) and the proposed method (COPULA-E: AMH-E, GUMH-E and Joe-E). Table 4 presents the measured values of the first three outcome metrics for the efficiency of a dose-finding design. The three outcome metrics for the safety aspects are depicted in Web Table 2 in Web Appendix B.
(1) Average performance: Across all the 18 scenarios, the MTDC selection percentages were observed to be, respectively, 34% and 44% for AMH-O and AMH-E (10% improvement); 32% and 41% for GUMH-O and GUMH-E (9% improvement); and 33% and 38% for Joe-O and Joe-E (5% improvement). Moreover, AMH-E, GUMH-E, and Joe-E exhibited OC selection percentages that were 16%, 14%, and 7% smaller than those exhibited by AMH-O, GUMH-O, and Joe-O, respectively. AMH-E, GUMH-E, and Joe-E exhibited improvements by 17, 18, and 10 in AI compared to AMH-O, GUMH-O, and Joe-O, respectively. Overall, AMH-E performed the best with respect to the majority of outcome metrics, including the percentage of correct MTDC selection.
(2) Performance of each scenario: AMH-E improved the MTDC selection percentage by more than 10% compared to 4,5,6,8,12,14,15,16,and 17. There was no scenario in which AMH-E was inferior to AMH-O by greater than 4% in terms of this metric. The OC selection percentage achieved by AMH-E was lower than that of AMH-O by more than 10% in all scenarios except 2, 3, and 11. GUMH-E also improved the MTDC selection percentage by more than 10% compared to GUMH-O in scenarios 2, 4, 5, 6, 8, 10, 12, and 15. The OC selection percentage achieved by GUMH-E was less than that of GUMH-O by at least 10% in all scenarios except 2, 3, 11, and 12. Joe-E improved the MTDC selection percentage by more than 10% compared to Joe-O in scenarios 5, 6, 8, 10, and 15. The OC selection percentage achieved by Joe-E was less than that of Joe-O by at least 10% in scenarios 9, 10, 13, 15, 16, and 18. More importantly, in all scenarios suggesting the presence of extreme synergistic toxicity, AMH-E and GUMH-E exhibited substantial improvements over AMH-O and GUMH-O, respectively. Furthermore, COPULA-E exhibited better MTDC selection percentages corresponding to all three copula-type models even in scenarios 17 and 18, which involved no synergistic toxicity.
(1) Average performance: Across all 18 scenarios, all dosefinding designs except TDM exhibited similar performances, with MTDC selection percentages ranging between 43% and 45%. TDM performed the worst with a score of 34%. In addition, AMH-E and LOGI exhibited the lowest percentage of OC selection at 28%, followed by POCRM (29%) and BOIN (30%) and, finally, TDM (42%). Similarly, TDM exhibited the lowest AI value of 36 and the other designs exhibited similar performance, with AIs ranging between 53 and 56. Although LOGI performed the best in terms of the first three outcome metrics on average across 18 scenarios, the differences between the performances of most of the methods were small. (2) Performance of each scenario: The MTDC selection percentages of AMH-E, LOGI, and BOIN were observed to be comparable in each scenario. However, BOIN exhibited an OC selection percentage that was 10% higher than those of AMH-E and LOGI in scenarios 1, 12, 15, and 17. The MTDC selection percentage of POCRM was 10% lower than those of AMH-E, LOGI, and BOIN in scenarios 2 and 4, even though it exhibited the highest MTDC selection percentages in scenarios 1, 9, 10, 13, and 16. TDM was observed to perform poorly in most scenarios in terms of both MTDC and OC selection percentages. The values of AI in each scenario were in agreement with both the MTDC and OC selection percentages. AMH-E was not observed to be significantly inferior to any of the other dose-finding designs in terms of any metric in any scenario.

DISCUSSION AND CONCLUSION
The performance of the dose-finding designs based on existing copula-type models have been reported to be somewhat inferior to those of other dose-finding designs (Wages et al., 2011;Riviere et al., 2014;Mander and Sweeting, 2015;Yin, 2016, 2017). Hirakawa et al. (2015) reported the MTDC selection percentage of the dose-finding design using the Clayton copulatype model to be 34% on average across 16 scenarios, which was approximately 10% lower than that of other dose-finding designs. This is primarily attributed to the presence of extreme synergistic toxicity in the majority of scenarios considered by the authors. In particular, copula-type models are incapable of accurately estimating joint toxicity probabilities in such scenarios due to F-H bounds, which reduce the MTDC selection percentages.
In this study, we proposed a relatively more flexible method to estimate toxicity probabilities by extending the parameter space of in copula-type models. For dose-finding designs based on the copula-type models, AMH, GUMH, and Joe, we conducted numerical studies corresponding to 18 scenarios, including the 16 scenarios considered by Hirakawa et al. (2015). The simulation results revealed that the MTDC selection percentages were improved by approximately 10% by using the proposed method (AMH-E and GUMH-E) compared with the utilization of COPULA (AMH-O and GUMH-O) on average across 18 scenarios. Further, we verified that the performance of AMH-E is comparable to those of three other existing dose-finding designs (LOGI, POCRM, and BOIN) and superior to that of TDM. In addition, we evaluated the absolute performance of dose-finding designs featured in the paper using the nonparametric optimal benchmark proposed by Mozgunov et al. (2021). Based on the benchmark assessment, the best methods are LOGI and AMH-E, followed by POCRM, BOIN, GUMH-E, and Joe-E. The original methods, AMH-O, GUMH-O, and Joe-O, and TDM, perform poorly compared to the other methods (see Web Appendix C).
In order to evaluate the applicability of the proposed method, further numerical studies were conducted. We verified that the stopping rule of the proposed method worked well under the 18 scenarios listed in Table 3 and four additional scenarios where all dose combinations are considered to be overly toxic (see Web Appendix D). Further, in the four scenarios where the toxicity probabilities of all dose combinations are extremely low, we verified that the proposed design escalates the dose to the highest dose combination as smoothly as the original design (see Web Appendix E).
For AMH-E, we conducted a further numerical study by applying different standard deviation (SD) of 50, 100, 500, 1000, 1500, and 3000 to a prior distribution of . The corresponding results are presented in Web Table 10 in Web Appendix F. In ascending order of SD values, the average MTDC selection percentages and OC selection percentages over all 18 scenarios were 40%, 40%, 42%, 43%, 44%, 44%, and 34%, 33%, 30%, 29%, 28%, and 28%, respectively (SD = 1500 was applied in Section 3). It appears that AMH-E with a prior distribution of assuming a larger drug-drug interaction exhibited better performance. The averaged posterior means of over the 18 scenarios were −38, −76, −369, −729, −1084, and −2144, respectively. Although = 0 indicates no synergistic toxicity in the AMH copula, the posterior means were quite far from 0 even in scenarios 17 and 18, which involved no synergistic toxicity. The posterior mean was strongly influenced by the prior settings. We evaluated the bias of a toxicity probability estimator. The results indicate that the AMH-E can estimate toxicity probabilities with less bias on average than the AMH-O with or without the existence of extreme synergistic toxicity (see Appendix G). In addition, the larger the SD, the smaller the bias. However, if a prior distribution of with a larger SD is considered, the estimated marginal toxicity probabilities tend to be smaller. Given these insights, we recommend using a prior distribution of with a larger SD in AMH-E (e.g., SD between 500 and 1500). However, if the estimation of the marginal toxicity probabilities is important, the following two solutions should be considered. The first is the use of smaller SDs for prior distributions of , although the design performance is expected to be somewhat reduced in this case. The second is to seek more appropriate sets of , and that do not produce unexpectedly large DLT probabilities, while ensuring the production of predictive values similar to the DLT probabilities predicted under a prior distribution of with a larger SD (e.g., SD = 1500). This approach is somewhat complicated, but we confirmed that it achieves similar design performance to that obtained with a prior distribution of with a larger SD (see Web Appendix H). Similarly, GUMH-E and Joe-E with a prior distribution of presuming a large drug-drug interaction exhibited better performance (see Web Tables 11 and 12 in Web Appendix F).
Numerous dose-finding designs for drug combinations have been proposed over the years, however, none of them are decidedly better than the other in terms of crucial outcome measures (e.g., the percentage of correct MTDC selection). Considering the growing popularity of combining anticancer drugs with different action mechanisms in cancer treatment, more complicated dose-finding designs may be required in the future. In the absence of exceptional dose-finding designs, development of dosefinding designs with additional features, such as BOIN and Keyboard design (Pan et al., 2020), which are easily implementable and exhibit comparable performance to the other designs, may be beneficial. In particular, as with copula-type models, the proposed method exhibits distinctive additional advantages as dose-finding designs as they are applicable to situations in which certain toxicities can be attributed to one of the two agents and they can also naturally handle binary DLT responses when either of the two drugs is administered as a single agent. Future studies are warranted to apply the proposed method to such situations. However, considering future expandability, the results obtained in this study could motivate the real-world application of the proposed method in cases requiring the utilization of the properties of copula-type models.

A C K N O W L E D G M E N T S
We would like to thank the associate editor and the two referees for their many insightful and constructive comments that substantially improved the article.

D ATA AVA I L A B I L I T Y S TAT E M E N T
Data sharing is not applicable to this article as no new data were created or analyzed in this paper.

O P E N R E S E A R C H B A D G E S
This article has earned an Open Materials badge for making publicly available the components of the research methodology needed to reproduce the reported procedure and analysis.