A nonparametric Bayesian continual reassessment method in single-agent dose-finding studies.

BACKGROUND: The main purpose of dose-finding studies in Phase I trial is to estimate maximum tolerated dose (MTD), which is the maximum test dose that can be assigned with an acceptable level of toxicity. Existing methods developed for single-agent dose-finding assume that the dose-toxicity relationship follows a specific parametric potency curve. This assumption may lead to bias and unsafe dose escalations due to the misspecification of parametric curve.
METHODS: This paper relaxes the parametric assumption of dose-toxicity relationship by imposing a Dirichlet process prior on unknown dose-toxicity curve. A hybrid algorithm combining the Gibbs sampler and adaptive rejection Metropolis sampling (ARMS) algorithm is developed to estimate the dose-toxicity curve, and a two-stage Bayesian nonparametric adaptive design is presented to estimate MTD.
RESULTS: For comparison, we consider two classical continual reassessment methods (CRMs) (i.e., logistic and power models). Numerical results show the flexibility of the proposed method for single-agent dose-finding trials, and the proposed method behaves better than two classical CRMs under our considered scenarios.
CONCLUSIONS: The proposed dose-finding procedure is model-free and robust, and behaves satisfactorily even in small sample cases.

Background

Dose-finding designs for phase I trials have been widely discussed over the past two decades. Many methods have been proposed to identify maximum tolerated dose (MTD) in single-agent dose-finding clinical trials. There are two major branches: model-based and algorithm-based methods. Among the model-based methods, the continual reassessment method (CRM) proposed by O’Quigley et al. [1] is a quite popular dose-finding approach. Since O’Quigley et al.’s pioneer work, there is considerable literature on modifying or improving CRM. For example, see Whitehead and Brunier [2] for Bayesian decision-theoretic approach, Piantadosi et al. [3] for a modified CRM, Heyd and Carlin [4] for an adaptive design improvement of the CRM, Leung and Wang [5] for a decision-theory-based extension of the CRM, Yuan et al. [6] for the quasi-likelihood approach, Yin and Yuan [7] for Bayesian model averaging CRM, Møller [8] for the up-and-down design, Fan et al. [9] for a simple Bayesian decision-theoretic design. Recently, Morita et al. [10] extended the CRM to a hierarchical Bayesian dose-toxicity model that borrows strength between subgroups under the assumption that the subgroups are exchangeable. However, in many practical applications, the true dose-toxicity curve is unknown, thus specifying an explicit dose-toxicity curve via some parametric model is artificial and may largely increase complications.

The algorithm-based approach, which can be regarded as a “nonparametric” or model-free method, has received considerable attention over the past years. There is considerable literature on model-free methods from a Bayesian nonparametric viewpoint. For example, Gelfand and Kuo [11] developed a fully Bayesian nonparametric approach to estimate the dose-toxicity curve with the Dirichlet process (DP) prior and the product-of-beta prior in the quantal bioassay problem. Mukhopadhyay [12] adopted the DP prior to make inference on the dose level with a prespecified response rate. Gasparini and Eisele [13] developed a curve-free method to improve CRM by modeling toxicity probabilities with the product-of-beta prior. The product-of-beta prior is conjugate with the binomial distribution, but it has a well-known numerical problem for exact computation and involves many hyperparameters to be determined, which may affect the performance of the considered method [14]. Ivanova and Wang [15] developed a nonparametric approach to estimate MTD for each of two strata of patients. Yan et al. [16] presented a Bayesian design for dose escalation and de-escalation. To our knowledge, no literature exists on Bayesian nonparametric dose-finding approach in single-agent clinical trials using the Dirichlet process prior.

The main purpose of this paper is to develop a flexible nonparametric statistical model to estimate unknown dose-toxicity curve for single-agent phase I clinical trials from a Bayesian viewpoint. This is implemented by treating toxicity probabilities as unknown parameters of interest which is subject only to monotonicity assumption, which is equivalent to assuming that the dose-toxicity curve is an arbitrary right continuous nondecreasing function whose range lies in [0,1], and then specifying a Dirichlet process prior for the unknown dose-toxicity curve. Also, a two-stage Bayesian nonparametric adaptive dose-finding design is developed to estimate MTD in single-agent dose-finding clinical trials. The main merits of the proposed method include that (i) the assumption of a specific dose-toxicity curve in CRM is not necessary; (ii) there are only two hyperparameters in the specified DP prior; (iii) there is more information that can be used to estimate toxicity probabilities and MTD; (iv) the escalation or deescalation of the current dose to the adjacent dose is only implemented once, thus the highest toxicity dose can be adaptively reached; (v) doses with relatively high toxicity probability may have no chance to be assigned to patients, which guarantees the safety of patients.

The rest of this paper is organized as follows. In “Methods” section, we briefly review the traditional CRM in single-agent dose-finding clinical trials, and introduce Bayesian nonparametric probability model by imposing a DP prior on the unknown dose-toxicity curve. A hybrid algorithm combining the Gibbs sampler and adaptive rejection Metropolis sampling (ARMS) algorithm is developed to estimate toxicity probabilities in “Methods” section. In “Dose-finding algorithm” section, a two-stage Bayesian nonparametric adaptive dose-finding algorithm is presented to estimate MTD in single-agent dose-finding clinical trials. Simulation studies are conducted to investigate the finite sample performance of the proposed method in “Simulation study” section. Some concluding remarks are given in “Results” section.

Methods

Let N be the total number of patients, J be the total number of cohorts, m be the number of patients in each of J cohorts. Let d1toxicity probability, the dose-toxicity curve F(·) be the distribution of dose levels, p=F(d) be the toxicity probability at dose level d for k=1,…,K. It is assumed that there are j cohorts enrolled in the current experiment. Let and , where is the number of patients treated at dose level d in the first j cohorts, and represents the number of patients experienced the dose-limiting toxicity (DLT) among patients treated at dose level d. Let S be the number of non-zero components in n. Thus, we have S≤K. That is, only doses are administered to the first j cohorts until now. Let be the data collected from the first j cohorts of patients, be a set of toxicity probabilities corresponding to . For simplicity, in what follows, we omit superscript j in n and y.

The continual reassessment method (CRM)

Let p01<…toxicity probabilities associated with dose levels d1<…17], θ is the target toxicity probability, k is the location of the target dose level and K is the total number of dose levels. More discussions on the guess of the “skeleton” of the CRM see [17]. In the Bayesian framework, the CRM assumes that the dose-toxicity relationship F(d) is a strictly increasing function with respect to the dose level d, and is usually specified by a parametric model F(d;β) in which β is an unknown parameter to be estimated, such as, logistic function, power function, and hyperbolic function.

Let x∈{d1,…,d} be the dose assigned to the ith cohort for i=1,…,j. Suppose that y patients have experienced DLT among n patients treated at dose level x. Thus, y follows the binomial distribution with the toxicity probability F(x;β). Let D={y1,…,y} be the dataset collected from the first j cohorts of patients. The likelihood function for D is given by

Let f(β) be the prior density of β. It follows from the Bayesian theorem that a Bayesian estimator of β is given by where is the expectation taken with respect to the posterior probability density function . The (j+1)th cohort of patients is assigned to the dose level where θ is the target toxicity probability. Once the maximum sample size N is reached, the dose with the toxicity probability closest to θ is selected as the MTD.

Nonparametric probability model

The model-based approaches such as CRM and the modified CRMs are widely used in single-agent dose-finding clinical trials. However, parametric modeling of F(·) may be problematical when there is little prior information on the shape of the dose-toxicity curve [12]. Moreover, in some experimental situations, it is recognized that parametric models do not work [18]. To address these issues, we assume that F(·) is an arbitrary right continuous nondecreasing unknown function whose range is [0,1], and can be approximated using some nonparametric method. To this end, a Dirichlet process (DP) prior is adopted to approximate F(·). To wit, we assume F(d)∼DP(αF0(x|η)), where F0(x|η) is a base distribution that serves as a starting point for constructing a nonparametric distribution for F(d), η is a parameter associated with the base distribution, α represents the weight that a researcher assigns a priori to the base distribution, and reflects the closeness of F0(x|η) to F(d). Thus, the hierarchical model for the dataset {y1,…,y} is given by

In many applications, α is usually unknown. Many methods have been proposed to select the prior of α (e.g., see [19, 20]). Here, we assume that the prior distribution of α follows a Gamma distribution with hyperparameters a and b, i.e., α∼Γ(a,b). The base distribution F0(x|η) is a prior guess of the unknown distribution F(·). For computational simplicity, F0(x|η) should be a conjugate prior, such as a uniform distribution or a mixture of two types of distributions. Also, one may consider empirical Bayesian and noninformative prior for F0(x|η). More discussions on the selection of F0(·) see Mukhopadhyay [12]. Here, we take F0 as the cumulative distribution function of some standard distribution, and choose some appropriate hyperparameters such that the median of F0 should be consistent with the initial guess for the MTD. For example, if there is a prior belief that dose level d is the MTD, we can select appropriate values of μ and σ such that , where η={μ,σ}, μ is associated with the mean of F0, and σ is related with the standard variance of F0. In this case, similar to [12], we assume that the priors of μ and σ discretely follow the uniform distribution in the intervals [μ0−r,μ0+r] (e.g., taking ten equally spaced points in [μ0−r,μ0+r]) and [σ0−r,σ0+r] (e.g., taking ten equally spaced points in [σ0−r,σ0+r]) with r∈{1,2}, respectively, where μ0 and σ0 are the mean and standard variance of F0, respectively.

According to the definitions of μ and σ, it is logical to assume that the joint prior density of μ and σ has the form of π(μ,σ)=π(μ)π(σ). Thus, the hierarchical model can be written as

By the definition of the DP prior, for any finite measurable partition B1,…,B in the support of F0, the probability vector (F(B1),…,F(B))⊤ follows the Dirichlet distribution with parameter vector (αF0(B1),…, αF0(B))⊤. For the base distribution F0(x|η)=N(x;μ,σ), if we consider the following finite measurable partition for the support (−∞,∞): B1=(−∞,d1), B2=[d1,d2),…,B=[d,d), B=[d,+∞), thus we have αF0(B1)=αF0(d1),…,αF0(B)=α[1−F0(d)], and F(B1)=p1, F(B2)=p2−p1,…,F(B)=1−p. Then, the joint prior density π(p) of p=(p1,…,p)⊤ can be written as

where γ=α{F0(d)−F0(d)} for k=1,…,K+1, d0=−∞, d=∞, F0(−∞)=0, F0(∞)=1, p0=0, and p=1.

Conditional distributions

Under the above assumption, the likelihood function of the first j (j=1,…,J) cohorts of patients (i.e., the dataset ) is given by

It follows from Eq. (3) that the joint prior density π(p) of can be expressed as

where γ=α{F0(d)−F0(d)} for i=1,…,S+1, d0=−∞, , F0(−∞)=0, F0(∞)=1, p0=0, and . It is easily seen that Eq. (5) reduces to Eq. (3) when S=K.

It follows from Eqs. (4) and (5) that the joint posterior probability density of p given the dataset can be expressed as

Again, it follows from Eqs. (2) and (3) that the joint conditional distribution of {α,μ,σ} given is given by which indicates that the conditional distributions of α, μ and σ given have the following expressions:

It is easily seen from Eqs. (6) and (7) that the conditional distributions of p and α are not standard distributions. Thus, it is quite difficult to draw observations from these conditional distributions. To address the issue, the Gibbs sampler is employed to simulate observations from these conditional distributions. Thus, Bayesian estimate of p can be obtained by the simulated observations via the Gibbs sampler.

Implementation of Gibbs sampler

Let p− be a subset vector of p with the ith element of p deleted. It follows from Eq. (6) that the conditional distribution of p given p− and can be written as

where p

The Gibbs sampler for drawing observations is implemented as follows.

Step 1. Initialize α(0),μ(0),σ(0), , and let p0=0 and .

Step 2. At the (q+1)th iteration with current values , generate α∗ from the uniform distribution U(1,20) and u from the uniform distribution U(0,1), respectively, if , we let α(=α∗; otherwise, we let α(=α(.

Step 3. Generate μ( from the conditional distribution ;

Step 4. Generate σ( from the conditional distribution ;

Step 5. Simulate from the conditional distribution: .

Step 6. Generate from the conditional distribution:

Step+4. Draw from the conditional distribution: .

Step+5. Repeat Step 2 to Step S+4 until the convergence of the algorithm.

It is easily seen from Eq. (8) that it is impossible to directly draw observations from the conditional distribution of p given because of nonstandard distribution involved. To address the issue, the Adaptive Rejection Metropolis Sampling (ARMS) is employed to sample observations from the conditional distribution (8).

To implement the ARMS algorithm, we require constructing a proposal distribution from which one can easily sample observations. Let denote a set of dose levels in the interval [p,p]. For 1≤i≤j≤n0, let l(x) be the straight line passing through two points (x, logπ(x)) and (x, logπ(x)); otherwise, it is assumed that l(x) is not defined. We construct the following piecewise function h(x):

When l1(x) is undefined but l2(x) is defined above, we set max(l1,l2)= max(l2,l1)= min(l1,l2)= min(l2,l1)=l2. When both l1(x) and l2(x) are undefined, we take max(l1,l2)= max(l2,l1)= min(l1,l2)= min(l2,l1)=0. Under the above assumption, the proposal distribution g(x) is defined as g(x)= exp(h(x))/ω, which is a piecewise exponential distribution, where . Let be the simulated observation of p at the qth iteration. Then, the ARMS algorithm for drawing from the conditional distribution (8) at the (q+1)th iteration is as follows.

Step 1. Initialize n0, which is the number of points in the interval [p,p], and determine the set of dose levels .

Step 2. Generate from the proposal distribution g(x).

Step 3. Generate u from the uniform distribution U(0,1).

Step 4. If , we set , and relabel the points in in ascending order and let n0=n0+1, then go to step 2; otherwise, we set .

Step 5. Generate u∗ from the uniform distribution U(0,1).

Step 6. If , we take ; otherwise, we set .

Dose-finding algorithm

In this section, we develop a two-stage Bayesian nonparametric adaptive dose-finding design based on a two-stage procedure and the above proposed hybrid algorithm. Let c and c be the threshold values for dose escalation and de-escalation, respectively. In our numerical illustration, c and c satisfying the restriction that c+c>1 can be selected by the data-dependent approach via simulation studies rather than some fixed values such that the trail has some desirable operating characteristics, for example, a relatively high accuracy index, which is defined in Equation (6.1) of Cheng (2011). For safety, we restrict the dose escalating or de-escalating for the next cohort of patients to only one dose level of change at a time. The two-stage Bayesian nonparametric adaptive dose-finding design is described as follows.

(I) The start-up stage

Patients in the first cohort are administered to the lowest dose level d1. If at least one toxicity is observed, the first stage is stopped and the second stage is conducted. If no toxicity is observed, patients in the second cohort are administered to the dose level d2. This process is continued until at least one toxicity is observed.

(II) The second stage

For dose level d administered to the jth cohort of patients, if there is at least one toxicity outcome observed, we denote as the dose level administered to the jth cohort of patients, i.e., the superscript j denotes the numerical order of cohort and the subscript k represents the dose level.

(i) Patients in the (j+1)th cohort are administered to dose level . Based on the data of the first j+1 cohorts of patients, we can simultaneously obtain and the toxicity probability at the current dose level via the above proposed hybrid algorithm.

(ii) If the probability , the dose level administered to the (j+2)th cohort of patients is escalated to the dose level d, i.e. . If the current dose , the dose administered to the (j+2)th cohort of patients is still d, that is, .

(iii) If the probability , the dose level administered to the (j+2)th cohort of patients is deescalated to the dose level d, i.e. . If the current dose , the dose level administered to the (j+2)th cohort of patients is still d1, that is, .

(iv) Otherwise, the (j+2)th cohort of patients continues to be treated at the dose level d i.e., .

(v) Once the maximum sample size N is attained, the dose level with the probability of toxicity being closest to θ is selected as the MTD.

Simulation study

To investigate the finite sample performance of the nonparametric continual reassessment method (NCRM), several simulation studies are conducted for six toxicity scenarios together with eight dose levels (i.e., K=8), which are given in Table 1. Here, for simplicity, the dose levels are identified by a number from 1 to 8, that is, d=k for k=1,…,8; and we take the target toxicity probability as θ=0.3. In Table 1, Scenario 1 indicates that the toxicity probability is steadily increasing and the target dose is level 5; Scenario 2 shows that the first five dose levels are the same as those given in Scenario 1, but the toxicity probability for dose level 6 suddenly jumps to an unacceptably high level 0.6; Scenario 3 has a flat relationship with the toxicity probability never attaining unacceptable levels, and the target dose is level 8; Scenario 4 has a flat relationship with the toxicity probability never reaching unacceptable levels, and there is no target dose but level 8 is quite close to the target dose given in Scenario 3; Scenario 5 implies that all the dose levels have an unacceptable toxicity probability, and the target dose is level 1; Scenario 6 is used to examine the situation that all the doses are over toxic, and there is no target dose but level 1 is quite close to the target dose. In practice, the trial would be stopped and doses will be reformulated once too many toxicities are occurred. We take sample size of each trial as N=60, the number of cohorts as J=20, and m=3. In implementing the Gibbs sampler, we collect 1000 observations after 700 burn-in iterations.

Table 1

Six toxicity scenarios for a single-agent trial with θ=0.3

		Dose level
Scenario	Method	1	2	3	4	5	6	7	8
1	true	0.05	0.08	0.12	0.20	0.30 a	0.45	0.60	0.70
	skeleton	0.03	0.06	0.12	0.20	0.30	0.40	0.50	0.59
	F 0	0.006	0.02	0.07	0.16	0.31	0.50	0.69	0.84
2	true	0.05	0.08	0.12	0.20	0.30 a	0.60	0.80	0.90
	skeleton	0.002	0.01	0.06	0.16	0.30	0.45	0.59	0.71
	F 0	0.006	0.02	0.07	0.16	0.31	0.50	0.69	0.84
3	true	0.01	0.05	0.10	0.14	0.18	0.22	0.25	0.30 a
	skeleton	0.02	0.04	0.06	0.10	0.14	0.18	0.24	0.30
	F 0	0.04	0.05	0.08	0.12	0.16	0.21	0.27	0.34
4	true	0.01	0.05	0.08	0.12	0.16	0.2	0.24	0.26 a
	skeleton	0.003	0.01	0.03	0.05	0.10	0.15	0.22	0.30
	F 0	0.04	0.05	0.08	0.12	0.16	0.21	0.27	0.34
5	true	0.30 a	0.40	0.50	0.60	0.70	0.80	0.90	0.95
	skeleton	0.30	0.44	0.58	0.69	0.78	0.84	0.89	0.92
	F 0	0.31	0.40	0.50	0.60	0.69	0.77	0.84	0.89
6	true	0.40	0.45	0.50	0.55	0.60	0.65	0.70	0.80
	skeleton	0.30	0.40	0.50	0.59	0.67	0.74	0.80	0.84
	F 0	0.31	0.40	0.50	0.60	0.69	0.77	0.84	0.89

aNumbers in boldface are the target MTDs

For comparison, we consider two Bayesian model-based CRM dose-finding methods including the power model: and one-parameter logistic model: p= exp(a0+βd)/{1+ exp(a0+βd)} with a0=3, where the prior of β is assumed to follow the normal distribution with mean zero and variance 1.34, i.e., β∼N(0,1.34). For each of the above considered six scenarios, given θ, k and K, a good “skeleton” can be directly generated using the function ‘getprior’ in the R package with the optimal value of δ, which can be obtained by the algorithm of Lee and Cheung [17]. But, its computational burden is too expensive. To address this issue, we here use a fixed value, which is evaluated using the data-dependent approach via simulation studies from a prespecified indifference interval so that the trial has some desirable operating characteristics, for example, the prior expectation of η is close to its true values, to replace the optimal value of δ. For our considered six scenarios, simulation studies evidence that we can take δ as 0.05, 0.075, 0.03, 0.04, 0.07 and 0.05, respectively.

For each of the above considered six scenarios, we set the initial value p(0) as its corresponding skeleton in implementing Gibbs sampler. When α and F0 are assumed to be known, to investigate the effect of the selection of α, we consider three cases of α: α=5, 10 and 20, corresponding to small, moderate and large values of α, respectively, and a standard normal distribution assumption for F0. As mentioned above, we choose appropriate values of μ and σ such that median of F0 should be consistent with the initial guess for the MTD. To wit, if there is a prior belief that dose level d is the MTD, we can select appropriate values of μ and σ such that . Table 1 gives the values of F0’s corresponding to six scenarios together with eight dose levels, where μ=6 and σ=2 for scenario 1 and 2, μ=10 and σ=5 for scenario 3 and 4, and μ=3 and σ=4 for scenario 5 and 6, respectively.

The preceding proposed hybrid algorithm is used to calculate Bayesian estimates of p’s, and the preceding proposed two-stage Bayesian nonparametric adaptive dose-finding algorithm is employed to determine MTD. To illustrate how the NCRM works, we present results of one simulation trial for scenario 1 together with α=5 in Table 2. Included in this table are: i (the serial number of the current experiment cohort); x (the dose administered to the current cohort of patients); y, y, y (the observations for the current experiment cohort); (estimate of toxicity probability for k=1,…,K). Examination of Table 2 shows that (i) dose level 5 (i.e., d5) is selected as the MTD and , which is quite close to true toxicity probability 0.3; (ii) among 11 cohorts of patients administered to dose level d5, only 10 patients (i.e., ) are experienced toxicity; (iii) the doses with the relatively high toxicity probability (such as d7 and d8) may have no chance to be administered to patients, which guarantees the safety of patients.

Table 2

Bayesian estimates of p’s via NCRM under scenario 1

True					0.05	0.08	0.12	0.20	0.30	0.45	0.60	0.70
i	x i	y i1	y i2	y i3	p^_1	p^_2	p^_3	p^_4	p^_5	p^_6	p^_7	p^_8
1	d 1	0	0	1
2	d 1	0	0	0	0.0156	-	-	-	-	-	-	-
3	d 2	0	0	0	0.0278	0.0380	-	-	-	-	-	-
4	d 3	1	0	0	0.0439	0.0573	0.0789	-	-	-	-	-
5	d 4	0	0	0	0.0501	0.0703	0.1014	0.1336	-	-	-	-
6	d 5	0	0	0	0.0585	0.0809	0.1169	0.1534	0.2109	-	-	-
7	d 6	0	0	1	0.0589	0.0852	0.1248	0.1682	0.2377	0.3426	-	-
8	d 5	0	1	0	0.0595	0.0734	0.1196	0.1672	0.2430	0.3457	-	-
9	d 6	1	1	1	0.0711	0.1018	0.1496	0.2006	0.2841	0.4060	-	-
10	d 5	0	1	0	0.0763	0.0954	0.1484	0.2023	0.2864	0.4063	-	-
11	d 5	0	1	1	0.0681	0.0995	0.1528	0.2112	0.3116	0.4156	-	-
12	d 5	0	0	1	0.0814	0.1144	0.1668	0.2227	0.3190	0.4196	-	-
13	d 4	0	0	0	0.0702	0.0979	0.1474	0.1965	0.3072	0.4144	-	-
14	d 5	0	0	1	0.0670	0.0979	0.1431	0.1958	0.3104	0.4151	-	-
15	d 5	0	0	0	0.0672	0.0960	0.1421	0.1915	0.2904	0.4062	-	-
16	d 5	0	0	1	0.0767	0.1073	0.1549	0.2026	0.2989	0.4094	-	-
17	d 5	0	1	0	0.0704	0.1004	0.1450	0.1941	0.2958	0.4070	-	-
18	d 5	1	1	0	0.0723	0.1039	0.1546	0.2066	0.3196	0.4205	-	-
19	d 4	0	0	0	0.0632	0.0916	0.1369	0.1842	0.3128	0.4170	-	-
20	d 5	0	0	0	0.0730	0.1046	0.1428	0.1868	0.2977	0.4103	-	-

Figures 1 and 2 plot Bayesian nonparametric estimations of unknown dose-toxicity curve for three specified values of α for scenarios 1 and 3, respectively. In each figure, “ −⋆” represents the dose-toxicity curve corresponding to true dose-toxicity data, “ −∘” corresponds to the base curve, and “”, “” and “-+” correspond to the estimated dose-toxicity curves for α=5, 10 and 20, respectively. Examination of Figs. 1 and 2 show that the estimated curves are more and more close to F0 with the increase of α when the doses administered to patients are less than but close to the target dose, which is consistent with the conclusion that the large value of α reflects a prior belief that F is tight around F0.

Fig. 1

Nonparametric Bayesian estimation of dose-toxicity curve under scenarios 1 when α and F0 are known

Fig. 2

Nonparametric Bayesian estimation of dose-toxicity curve under scenarios 3 when α and F0 are known

For the aforementioned NCRM, we calculate the selection probabilities that a dose level is selected as the MTD, the total numbers of toxicities observed, and average numbers of patients that are administered to each of eight dose levels for the above specified six scenarios together with α=5 and 20. For comparison, we also calculate the results corresponding to logistic model and power model. Results for 1000 simulated trials are given in Tables 3 and 4. Examination of Tables 3 and 4 show that the above developed Bayesian NCRM performs better than two parametric CRMs (i.e., logistic and power models) in terms of the following four aspects: (i) Bayesian NCRM generally selects the MTD with a relatively higher probability than two parametric CRMs; (ii) Bayesian NCRM consistently selects over-MTD with a relatively lower probability than two parametric CRMs; (iii) the total number of toxicities observed are almost identical for all three methods under our considered cases; (iv) Bayesian NCRM has a higher percentage of patients treated at MTD than two parametric CRMs, except for scenario 3 with α=5, but Bayesian NCRM treats more patients at dose levels below MTD and less patients at dose levels above MTD than two parametric CRMs; (v) the selection probabilities for scenarios 3 and 4 are smaller than those for other four scenarios because the locations of target doses for scenarios 3 and 4 are different from those for other scenarios, which indicates that the highest dose should be carefully administered to patients for safety, and patients should be administered to the lower dose level.

Table 3

Selection probabilities and total numbers of toxicities observed for logistic model, power model and NCRM under six scenarios

		Dose level								# of	# of
Scenario	Method	1	2	3	4	5	6	7	8	Tox.	Pat.
1	Logistic	0.000	0.000	0.008	0.227	0.637	0.128	0.000	0.000	15	60
	Power	0.000	0.000	0.004	0.212	0.664	0.120	0.000	0.000	15	60
	NCRM5a	0.000	0.002	0.001	0.161	0.773	0.063	0.000	0.000	15	60
	NCRM2a	0.000	0.000	0.000	0.007	0.980	0.013	0.000	0.000	16	60
2	Logistic	0.000	0.001	0.042	0.362	0.571	0.024	0.000	0.000	13	60
	Power	0.000	0.001	0.028	0.369	0.581	0.021	0.000	0.000	14	60
	NCRM5	0.002	0.000	0.005	0.292	0.698	0.003	0.000	0.000	15	60
	NCRM2	0.000	0.000	0.000	0.186	0.814	0.000	0.000	0.000	14	60
3	Logistic	0.000	0.000	0.000	0.002	0.043	0.108	0.306	0.541	13	60
	Power	0.000	0.000	0.000	0.001	0.031	0.104	0.311	0.553	13	60
	NCRM5	0.000	0.001	0.007	0.019	0.028	0.117	0.373	0.455	13	60
	NCRM2	0.000	0.000	0.000	0.000	0.000	0.008	0.379	0.613	13	60
4	Logistic	0.000	0.000	0.000	0.002	0.036	0.143	0.286	0.533	11	60
	Power	0.000	0.000	0.000	0.000	0.028	0.128	0.321	0.523	11	60
	NCRM5	0.000	0.001	0.003	0.007	0.011	0.086	0.289	0.603	12	60
	NCRM2	0.000	0.000	0.000	0.000	0.000	0.004	0.263	0.733	12	60
5	Logistic	0.804	0.189	0.007	0.000	0.000	0.000	0.000	0.000	20	60
	Power	0.811	0.183	0.006	0.000	0.000	0.000	0.000	0.000	20	60
	NCRM5	0.898	0.100	0.002	0.000	0.000	0.000	0.000	0.000	19	60
	NCRM2	0.971	0.029	0.000	0.000	0.000	0.000	0.000	0.000	19	60
6	Logistic	0.970	0.027	0.003	0.000	0.000	0.000	0.000	0.000	24	60
	Power	0.972	0.028	0.000	0.000	0.000	0.000	0.000	0.000	24	60
	NCRM5	0.987	0.013	0.000	0.000	0.000	0.000	0.000	0.000	24	60
	NCRM2	1.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	24	60

aNCRM5 and NCRM2 denote NCRM method with α=5 and 20, respectively

Table 4

Average numbers of patients treated at each of eights doses for logistic model, power model and NCRM

		Dose level
Scenario	Method	1	2	3	4	5	6	7	8
1	Logistic	3.579	4.008	5.916	14.463	23.304	8.124	0.576	0.030
	Power	3.489	3.954	5.622	14.739	23.889	7.755	0.543	0.009
	NCRM5a	3.474	3.615	3.912	15.462	27.387	6.009	0.141	0.000
	NCRM2a	3.423	3.558	3.708	8.298	35.454	5.493	0.066	0.000
2	Logistic	3.582	4.632	8.868	19.206	21.147	2.535	0.030	0.000
	Power	3.486	4.158	8.004	20.496	21.663	2.178	0.015	0.000
	NCRM5	3.447	3.648	4.569	18.543	25.698	4.065	0.300	0.000
	NCRM2	3.447	3.513	3.693	19.977	27.825	1.497	0.048	0.000
3	Logistic	3.108	3.486	4.137	5.274	6.882	9.096	11.724	16.293
	Power	3.066	3.480	4.089	5.157	6.552	8.970	12.342	16.344
	NCRM5	3.090	3.435	3.684	3.897	4.908	8.901	17.457	14.628
	NCRM2	3.057	3.420	3.660	3.636	3.531	3.930	21.168	17.598
4	Logistic	3.105	3.483	3.942	4.977	7.266	10.125	11.610	15.492
	Power	3.066	3.468	3.843	4.737	6.573	10.266	12.861	15.186
	NCRM5	3.135	3.366	3.453	3.777	4.395	7.809	16.281	17.784
	NCRM2	3.084	3.387	3.573	3.612	3.522	3.630	18.879	20.313
5	Logistic	44.637	12.972	2.229	0.156	0.006	0.000	0.000	0.000
	Power	45.189	12.822	1.848	0.135	0.006	0.000	0.000	0.000
	NCRM5	49.740	8.958	1.176	0.114	0.012	0.000	0.000	0.000
	NCRM2	53.178	6.366	0.408	0.042	0.006	0.000	0.000	0.000
6	Logistic	54.468	4.467	0.894	0.141	0.030	0.000	0.000	0.000
	Power	54.573	4.500	0.810	0.102	0.015	0.000	0.000	0.000
	NCRM5	57.676	2.865	0.414	0.039	0.006	0.000	0.000	0.000
	NCRM2	58.149	1.596	0.237	0.018	0.000	0.000	0.000	0.000

aNCRM5 and NCRM2 denote NCRM method with α=5 and 20, respectively

Now we assume that α and F0 are unknown. As an illustration of the above presented Bayesian NCRM procedure, here we only consider scenarios 1 and 3. For scenario 1, we consider four different priors on α: Γ(2,5), Γ(2,2), Γ(5,2) and Γ(10,2), which correspond to small and moderate expectations of α, and four discrete uniform priors for (μ,σ) on the rectangles (5,7)×(1,3), (4,8)×(1,3), (5,7)×(0,4) and (4,8)×(0,4), which indicate that prior expectations of μ and σ are 6 and 2, respectively. For scenario 3, we consider the same priors on α as scenario 1, but the following four different discrete uniform priors of (μ,σ) on the rectangles (9,11)×(4,6), (8,12)×(4,6), (9,11)×(3,7) and (8,12)×(3,7), which imply that prior expectations of μ and σ are 10 and 5, respectively.

Based on the above considered settings, the preceding introduced hybrid algorithm is adopted to evaluate Bayesian estimation of p’s, and the preceding developed two-stage Bayesian nonparametric adaptive dose-finding algorithm is employed to determine the MTD. Similarly, we also calculate the selection probabilities, total numbers of toxicities observed and average numbers of patients treated at each of eight dose levels for scenarios 1 and 3. Results for 1000 simulated trials are given in Tables 5 and 6. Examination of Table 5 and 6 shows that (i) the selection probabilities and the number of patients treated at the MTD increase with the increase of prior expectation of α; (ii) the selection probabilities and the number of patients treated at the dose level closet to the MTD decrease with the increase of prior expectation of α for scenario 1, but increase with the increase of prior expectation of α for scenario 3; (iii) the total number of toxicities observed are almost equal regardless of the priors of α and (μ,σ), which shows that there is little effect of the selection of the priors of α, μ and σ on the total number of toxicities observed.

Table 5

Selection probabilities and total numbers of toxicities observed for NCRM under scenarios 1 and 3 when α and (μ,σ) are unknown

Prior			Dose Level								# of	# of
Case	(μ,σ)a	α	1	2	3	4	5	6	7	8	Tox.	Pat.
1	A	Γ(2,5)	0.012	0.016	0.016	0.200	0.646	0.110	0.000	0.000	15	60
		Γ(2,2)	0.002	0.001	0.001	0.167	0.771	0.057	0.001	0.000	15	60
		Γ(5,2)	0.000	0.001	0.001	0.092	0.853	0.053	0.000	0.000	15	60
		Γ(10,2)	0.000	0.000	0.000	0.064	0.896	0.040	0.000	0.000	16	60
	B	Γ(2,5)	0.009	0.008	0.015	0.236	0.647	0.081	0.004	0.000	14	60
		Γ(2,2)	0.000	0.000	0.000	0.139	0.817	0.044	0.000	0.000	15	60
		Γ(5,2)	0.000	0.000	0.000	0.125	0.826	0.049	0.000	0.000	15	60
		Γ(10,2)	0.000	0.000	0.000	0.092	0.861	0.047	0.000	0.000	15	60
	C	Γ(2,5)	0.000	0.000	0.002	0.172	0.776	0.050	0.000	0.000	15	60
		Γ(2,2)	0.000	0.000	0.000	0.162	0.783	0.055	0.000	0.000	15	60
		Γ(5,2)	0.000	0.000	0.001	0.107	0.831	0.061	0.000	0.000	15	60
		Γ(10,2)	0.000	0.000	0.000	0.064	0.881	0.055	0.000	0.000	16	60
	D	Γ(2,5)	0.000	0.000	0.004	0.219	0.726	0.050	0.001	0.000	15	60
		Γ(2,2)	0.000	0.000	0.001	0.177	0.783	0.039	0.000	0.000	15	60
		Γ(5,2)	0.000	0.000	0.001	0.152	0.796	0.051	0.000	0.000	15	60
		Γ(10,2)	0.000	0.000	0.000	0.129	0.834	0.037	0.000	0.000	15	60
3	E	Γ(2,5)	0.001	0.005	0.015	0.048	0.094	0.157	0.242	0.438	13	60
		Γ(2,2)	0.001	0.005	0.008	0.032	0.059	0.138	0.294	0.463	12	60
		Γ(5,2)	0.001	0.004	0.009	0.013	0.026	0.095	0.377	0.475	13	60
		Γ(10,2)	0.000	0.000	0.000	0.000	0.000	0.020	0.400	0.580	13	60
	F	Γ(2,5)	0.001	0.007	0.018	0.041	0.077	0.162	0.239	0.455	13	60
		Γ(2,2)	0.001	0.004	0.013	0.017	0.044	0.118	0.347	0.456	13	60
		Γ(5,2)	0.000	0.001	0.001	0.000	0.008	0.059	0.410	0.521	13	60
		Γ(10,2)	0.000	0.000	0.000	0.000	0.000	0.046	0.425	0.529	13	60
	G	Γ(2,5)	0.003	0.009	0.015	0.039	0.083	0.145	0.258	0.448	13	60
		Γ(2,2)	0.000	0.004	0.017	0.031	0.048	0.127	0.318	0.455	13	60
		Γ(5,2)	0.001	0.002	0.004	0.004	0.017	0.071	0.376	0.525	13	60
		Γ(10,2)	0.000	0.001	0.000	0.000	0.000	0.014	0.357	0.628	13	60
	H	Γ(2,5)	0.001	0.001	0.017	0.032	0.073	0.144	0.264	0.468	12	60
		Γ(2,2)	0.001	0.005	0.005	0.014	0.019	0.089	0.350	0.517	13	60
		Γ(5,2)	0.000	0.001	0.001	0.002	0.006	0.064	0.364	0.562	13	60
		Γ(10,2)	0.000	0.000	0.000	0.000	0.000	0.024	0.406	0.570	13	60

aNote: A=(5,7)×(1,3), B=(4,8)×(1,3), C=(5,7)×(0,4), D=(4,8)×(0,4), E=(9,11)×(4,6), F=(8,12)×(4,6), G=(9,11)×(3,7), H=(8,12)×(3,7)

Table 6

Average numbers of patients treated at each of eight doses for NCRM under scenario 1 and 3 when α and (μ,σ) are unknown

Prior			Dose Level
Case	(μ,σ)a	α	1	2	3	4	5	6	7	8
1	A	Γ(2,5)	3.474	3.915	6.273	15.879	21.567	7.983	0.879	0.030
		Γ(2,2)	3.342	3.606	4.263	15.705	26.751	6.030	0.303	0.000
		Γ(5,2)	3.390	3.672	3.705	12.576	30.048	6.426	0.177	0.006
		Γ(10,2)	3.450	3.603	3.675	10.635	31.983	6.582	0.066	0.006
	B	Γ(2,5)	3.537	3.957	6.381	17.718	21.378	6.450	0.555	0.024
		Γ(2,2)	3.342	3.618	3.810	15.129	28.473	5.367	0.165	0.006
		Γ(5,2)	3.393	3.510	3.702	14.523	29.253	5.460	0.159	0.000
		Γ(10,2)	3.363	3.594	3.717	13.014	31.074	5.157	0.075	0.006
	C	Γ(2,5)	3.417	3.570	4.155	14.847	26.706	7.101	0.204	0.00
		Γ(2,2)	3.411	3.573	3.948	14.589	28.089	6.339	0.051	0.000
		Γ(5,2)	3.372	3.570	3.783	12.495	29.781	6.819	0.180	0.000
		Γ(10,2)	3.402	3.567	3.663	10.887	31.743	6.621	0.105	0.012
	D	Γ(2,5)	3.435	3.525	4.404	18.249	25.404	4.761	0.210	0.012
		Γ(2,2)	3.444	3.558	3.927	16.455	28.017	4.455	0.144	0.000
		Γ(5,2)	3.429	3.501	3.798	16.017	28.530	4.566	0.153	0.006
		Γ(10,2)	3.408	3.612	3.684	15.390	29.757	4.053	0.090	0.006
3	E	Γ(2,5)	3.075	3.399	3.687	4.335	6.141	9.936	15.465	13.962
		Γ(2,2)	3.108	3.486	3.714	4.305	5.871	9.753	16.014	13.749
		Γ(5,2)	3.078	3.363	3.810	3.870	4.554	8.478	18.222	14.625
		Γ(10,2)	3.090	3.396	3.633	3.684	3.576	4.617	21.420	16.584
	F	Γ(2,5)	3.090	3.471	3.741	4.326	6.333	10.383	15.444	13.212
		Γ(2,2)	3.111	3.375	3.735	4.113	5.418	9.525	17.049	13.674
		Γ(5,2)	3.081	3.450	3.633	3.645	3.801	6.630	19.911	15.849
		Γ(10,2)	3.105	3.387	3.657	3.708	3.588	5.532	21.4189	15.834
	G	Γ(2,5)	3.096	3.477	3.789	4.275	6.015	9.840	15.312	14.196
		Γ(2,2)	3.096	3.450	3.651	4.245	5.706	9.441	16.191	14.220
		Γ(5,2)	3.372	3.375	3.792	3.735	4.005	6.873	19.104	16.056
		Γ(10,2)	3.066	3.474	3.681	3.597	3.570	4.347	20.088	18.177
	H	Γ(2,5)	3.060	3.411	3.753	4.272	5.952	9.834	15.462	14.256
		Γ(2,2)	3.081	3.426	3.672	3.849	4.434	7.320	18.411	15.807
		Γ(5,2)	3.051	3.477	3.648	3.672	3.747	6.306	19.278	16.821
		Γ(10,2)	3.075	3.396	3.666	3.648	3.633	4.785	20.613	17.184

aNote: A=(5,7)×(1,3), B=(4,8)×(1,3), C=(5,7)×(0,4), D=(4,8)×(0,4), E=(9,11)×(4,6), F=(8,12)×(4,6), G=(9,11)×(3,7), H=(8,12)×(3,7)

Results

According to the above presented simulation study, we have the following results. First, the estimated MTD in single-agent dose-finding clinical trials via our proposed two-stage Bayesian nonparametric adaptive dose-finding algorithm is quite close to true toxicity probability under our considered settings. Second, the doses with the relatively high toxicity probability may have no chance to be assigned to patients, which guarantees the safety of patients. Third, the value of the weight α in the DP prior is a measure of a prior belief on the base distribution. Fourth, when the parameters in the DP prior are fixed, the proposed Bayesian NCRM behaves better than traditional model-based CRM; but when the parameters in the DP prior are unknown, the selection of the weight α has a positive effect on the selection probabilities and the number of patients treated at the MTD, while the selection of parameters in the DP prior has little effect on the total number of toxicities observed. In a word, numerical results show the flexibility of the proposed method for single-agent dose-finding trials, and the proposed method outperforms two classical CRMs under our considered scenarios.

Discussion

Although this manuscript only considers a single agent dose-finding design, the proposed Bayesian NCRM can be extended to two-agent dose-finding studies, which is our further research topic. On the other hand, this paper only considers the evaluation of the toxicity of novel drug treatment, i.e., phase I clinical trial, but the developed Bayesian NCRM procedure can be extended to Bayesian nonparametric phase I/II dose-finding trial design that simultaneously allows for toxicity and efficiency of novel drug treatment for precision medicine.

Conclusions

This paper proposes a Bayesian nonparametric continual reassessment method to estimate the MTD for a single-agent in phase I clinical trials. We relax the traditional parametric model assumption imposed on dose-toxicity relationship using a DP prior to approximate unknown distribution of dose-toxicity curve. A Bayesian method is developed to estimate toxicity probabilities of dose levels considered. A hybrid algorithm combining the Gibbs sampler and adaptive rejection Metropolis sampling algorithm is developed to generate observations from joint conditional distributions required in evaluating Bayesian estimates of toxicity probabilities of dose levels. A two-stage Bayesian nonparametric adaptive dose-finding design is developed to estimate the MTD. In the proposed Bayesian nonparametric adaptive dose-finding design, the dose administered to next cohort of patients can be escalated or deescalated to adjacent dose level more safety; generally, doses with a relatively higher toxicity probability may have no chance to be administered to patients, which guarantees the safety of patients. Simulation studies evidence that the proposed dose-finding procedure is model-free and robust, and performs better than two parametric models even in small sample sizes.