 Research article
 Open Access
 Published:
A nonparametric Bayesian continual reassessment method in singleagent dosefinding studies
BMC Medical Research Methodology volume 18, Article number: 172 (2018)
Abstract
Background
The main purpose of dosefinding 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 singleagent dosefinding assume that the dosetoxicity 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 dosetoxicity relationship by imposing a Dirichlet process prior on unknown dosetoxicity curve. A hybrid algorithm combining the Gibbs sampler and adaptive rejection Metropolis sampling (ARMS) algorithm is developed to estimate the dosetoxicity curve, and a twostage 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 singleagent dosefinding trials, and the proposed method behaves better than two classical CRMs under our considered scenarios.
Conclusions
The proposed dosefinding procedure is modelfree and robust, and behaves satisfactorily even in small sample cases.
Background
Dosefinding 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 singleagent dosefinding clinical trials. There are two major branches: modelbased and algorithmbased methods. Among the modelbased methods, the continual reassessment method (CRM) proposed by O’Quigley et al. [1] is a quite popular dosefinding 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 decisiontheoretic 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 decisiontheorybased extension of the CRM, Yuan et al. [6] for the quasilikelihood approach, Yin and Yuan [7] for Bayesian model averaging CRM, Møller [8] for the upanddown design, Fan et al. [9] for a simple Bayesian decisiontheoretic design. Recently, Morita et al. [10] extended the CRM to a hierarchical Bayesian dosetoxicity model that borrows strength between subgroups under the assumption that the subgroups are exchangeable. However, in many practical applications, the true dosetoxicity curve is unknown, thus specifying an explicit dosetoxicity curve via some parametric model is artificial and may largely increase complications.
The algorithmbased approach, which can be regarded as a “nonparametric” or modelfree method, has received considerable attention over the past years. There is considerable literature on modelfree methods from a Bayesian nonparametric viewpoint. For example, Gelfand and Kuo [11] developed a fully Bayesian nonparametric approach to estimate the dosetoxicity curve with the Dirichlet process (DP) prior and the productofbeta 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 curvefree method to improve CRM by modeling toxicity probabilities with the productofbeta prior. The productofbeta prior is conjugate with the binomial distribution, but it has a wellknown 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 deescalation. To our knowledge, no literature exists on Bayesian nonparametric dosefinding approach in singleagent clinical trials using the Dirichlet process prior.
The main purpose of this paper is to develop a flexible nonparametric statistical model to estimate unknown dosetoxicity curve for singleagent 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 dosetoxicity curve is an arbitrary right continuous nondecreasing function whose range lies in [0,1], and then specifying a Dirichlet process prior for the unknown dosetoxicity curve. Also, a twostage Bayesian nonparametric adaptive dosefinding design is developed to estimate MTD in singleagent dosefinding clinical trials. The main merits of the proposed method include that (i) the assumption of a specific dosetoxicity 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 singleagent dosefinding clinical trials, and introduce Bayesian nonparametric probability model by imposing a DP prior on the unknown dosetoxicity 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 “Dosefinding algorithm” section, a twostage Bayesian nonparametric adaptive dosefinding algorithm is presented to estimate MTD in singleagent dosefinding 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 d_{1}<d_{2}<⋯<d_{K} be a set of K prespecified dose levels for a single agent under investigation, θ be the target toxicity probability, the dosetoxicity curve F(·) be the distribution of dose levels, p_{k}=F(d_{k}) be the toxicity probability at dose level d_{k} for k=1,…,K. It is assumed that there are j cohorts enrolled in the current experiment. Let \(\boldmath {n}^{j}=\left (n_{1}^{j},\ldots,n_{K}^{j}\right)^{\!\top \!}\) and \(\boldmath {y}^{j}=\left (y_{1}^{j},\ldots,y_{K}^{j}\right)^{\!\top \!}\), where \(n_{k}^{j}\) is the number of patients treated at dose level d_{k} in the first j cohorts, and \(y_{k}^{j}\) represents the number of patients experienced the doselimiting toxicity (DLT) among \(n_{k}^{j}\) patients treated at dose level d_{k}. Let S_{j} be the number of nonzero components in n^{j}. Thus, we have S_{j}≤K. That is, only doses \(d_{1},\ldots,d_{\phantom {\dot {i}\!}S_{j}}\) are administered to the first j cohorts until now. Let \(D_{\phantom {\dot {i}\!}S_{j}}=\left \{\left (n_{1}^{j},y_{1}^{j}\right),\ldots, \left (n_{\phantom {\dot {i}\!}S_{j}}^{j},y_{\phantom {\dot {i}\!}S_{j}}^{j}\right)\right \}\) be the data collected from the first j cohorts of patients, \(\boldmath {p}=(p_{1},\ldots,p_{\phantom {\dot {i}\!}S_{j}})^{\!\top \!}\) be a set of toxicity probabilities corresponding to \(\{d_{1},\ldots,d_{\phantom {\dot {i}\!}S_{j}}\}\). For simplicity, in what follows, we omit superscript j in n^{j} and y^{j}.
The continual reassessment method (CRM)
Let p_{01}<…<p_{0K} be an initial guess of toxicity probabilities associated with dose levels d_{1}<…<d_{K}, where p_{0k} is referred to as the “skeleton” of the CRM. It is assumed that there is a prior belief that dose level d_{k} is the MTD, one may set the initial guess as p_{0k}=θ. In practice, clinical investigators attempt to provide their best guess on the “skeleton”. However, in simulation study, the “skeleton” can be directly generated using the function ‘ getprior(δ,θ,k,K)’ in the R package ‘dfcrm’, where δ is the desired length of the indifference interval and the optimal value of δ can be calibrated by the algorithm given in Lee and Cheung [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 dosetoxicity relationship F(d_{k}) is a strictly increasing function with respect to the dose level d_{k}, and is usually specified by a parametric model F(d_{k};β) in which β is an unknown parameter to be estimated, such as, logistic function, power function, and hyperbolic function.
Let x_{i}∈{d_{1},…,d_{K}} be the dose assigned to the ith cohort for i=1,…,j. Suppose that y_{i} patients have experienced DLT among n_{i} patients treated at dose level x_{i}. Thus, y_{i} follows the binomial distribution with the toxicity probability F(x_{i};β). Let D_{j}={y_{1},…,y_{j}} be the dataset collected from the first j cohorts of patients. The likelihood function for D_{j} 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 \(E_{\phantom {\dot {i}\!}\beta D_{j}}\) is the expectation taken with respect to the posterior probability density function \(f\left (\beta D_{j}\right)=L_{j}(\beta D_{j})f(\beta)/\int L_{j}(\beta D_{j})f(\beta)d\beta \). 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 modelbased approaches such as CRM and the modified CRMs are widely used in singleagent dosefinding clinical trials. However, parametric modeling of F(·) may be problematical when there is little prior information on the shape of the dosetoxicity 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_{k})∼DP(αF_{0}(xη)), where F_{0}(xη) is a base distribution that serves as a starting point for constructing a nonparametric distribution for F(d_{k}), η 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 F_{0}(xη) to F(d_{k}). Thus, the hierarchical model for the dataset {y_{1},…,y_{K}} 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 F_{0}(xη) is a prior guess of the unknown distribution F(·). For computational simplicity, F_{0}(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 F_{0}(xη). More discussions on the selection of F_{0}(·) see Mukhopadhyay [12]. Here, we take F_{0} as the cumulative distribution function of some standard distribution, and choose some appropriate hyperparameters such that the median of F_{0} should be consistent with the initial guess for the MTD. For example, if there is a prior belief that dose level d_{k} is the MTD, we can select appropriate values of μ and σ such that \(F_{0}(d_{k}\eta)=\Phi \left (\frac {d_{k}\mu }{\sigma }\right)=\theta \), where η={μ,σ}, μ is associated with the mean of F_{0}, and σ is related with the standard variance of F_{0}. 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 F_{0}, 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 B_{1},…,B_{K} in the support of F_{0}, the probability vector (F(B_{1}),…,F(B_{K}))^{⊤} follows the Dirichlet distribution with parameter vector (αF_{0}(B_{1}),…, αF_{0}(B_{K}))^{⊤}. For the base distribution F_{0}(xη)=N(x;μ,σ), if we consider the following finite measurable partition for the support (−∞,∞): B_{1}=(−∞,d_{1}), B_{2}=[d_{1},d_{2}),…,B_{K}=[d_{K−1},d_{K}), B_{K+1}=[d_{K},+∞), thus we have αF_{0}(B_{1})=αF_{0}(d_{1}),…,αF_{0}(B_{K+1})=α[1−F_{0}(d_{K})], and F(B_{1})=p_{1}, F(B_{2})=p_{2}−p_{1},…,F(B_{K+1})=1−p_{K}. Then, the joint prior density π(p) of p=(p_{1},…,p_{K})^{⊤} can be written as
where γ_{k}=α{F_{0}(d_{k})−F_{0}(d_{k−1})} for k=1,…,K+1, d_{0}=−∞, d_{K+1}=∞, F_{0}(−∞)=0, F_{0}(∞)=1, p_{0}=0, and p_{K+1}=1.
Conditional distributions
Under the above assumption, the likelihood function of the first j (j=1,…,J) cohorts of patients (i.e., the dataset \(D_{\phantom {\dot {i}\!}S_{j}}=\{(n_{1},y_{1}),\ldots, (n_{\phantom {\dot {i}\!}S_{j}},y_{\phantom {\dot {i}\!}S_{j}})\}\)) is given by
It follows from Eq. (3) that the joint prior density π(p) of \(\boldmath {p}=(p_{1},\ldots,p_{\phantom {\dot {i}\!}S_{j}})^{\!\top \!}\) can be expressed as
where γ_{i}=α{F_{0}(d_{i})−F_{0}(d_{i−1})} for i=1,…,S_{j}+1, d_{0}=−∞, \(d_{\phantom {\dot {i}\!}S_{j}+1}=\infty \), F_{0}(−∞)=0, F_{0}(∞)=1, p_{0}=0, and \(p_{\phantom {\dot {i}\!}S_{j}+1}=1\). It is easily seen that Eq. (5) reduces to Eq. (3) when S_{j}=K.
It follows from Eqs. (4) and (5) that the joint posterior probability density of p given the dataset \(D_{\phantom {\dot {i}\!}S_{j}}\) can be expressed as
Again, it follows from Eqs. (2) and (3) that the joint conditional distribution of {α,μ,σ} given \(\{D_{\phantom {\dot {i}\!}S_{j}},\boldmath {p}\}\) is given by
which indicates that the conditional distributions of α, μ and σ given \(\{\boldmath {p},D_{\phantom {\dot {i}\!}S_{j}}\}\) 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_{−i} be a subset vector of p with the ith element of p deleted. It follows from Eq. (6) that the conditional distribution of p_{i} given p_{−i} and \(D_{\phantom {\dot {i}\!}S_{j}}\) can be written as
where p_{i−1}<p_{i}<p_{i+1} for i=1,…,S_{j}, and p_{0}=0 and \(p_{\phantom {\dot {i}\!}S_{j}+1}=1\).
The Gibbs sampler for drawing observations \(\alpha,\mu,\sigma,p_{1},\ldots, p_{\phantom {\dot {i}\!}S_{j}}\) is implemented as follows.
Step 1. Initialize α^{(0)},μ^{(0)},σ^{(0)}, \(\boldmath {p}^{(0)}=\left (p_{1}^{(0)},p_{2}^{(0)},\ldots,\right. \left.p_{\phantom {\dot {i}\!}S_{j}}^{(0)}\right)\), and let p_{0}=0 and \(p_{\phantom {\dot {i}\!}S_{j}+1}=1\).
Step 2. At the (q+1)th iteration with current values \(\alpha ^{(q)}, \mu ^{(q)}, \sigma ^{(q)}, \boldmath {p}^{(q)}=\left (p_{1}^{(q)},p_{2}^{(q)}, \ldots,p_{\phantom {\dot {i}\!}S_{j}}^{(q)}\right)\), generate α^{∗} from the uniform distribution U(1,20) and u from the uniform distribution U(0,1), respectively, if \(u\leq \min \!\left \{\!1, \pi \!\!\left (\alpha ^{*}\boldmath {p}^{(q)},\!\mu ^{(q)}, \!\sigma ^{(q)},\!D_{\phantom {\dot {i}\!}S_{j}}\right)\!/\pi \!\left (\alpha ^{(q)}\boldmath {p}^{(q)},\!\mu ^{(q)},\sigma ^{(q)},\!D_{\phantom {\dot {i}\!}S_{j}}\right)\!\right \}\), we let α^{(q+1)}=α^{∗}; otherwise, we let α^{(q+1)}=α^{(q)}.
Step 3. Generate μ^{(q+1)} from the conditional distribution \(\pi \left (\mu \boldmath {p}^{(q)},\alpha ^{(q)}, \sigma ^{(q)}, D_{\phantom {\dot {i}\!}S_{j}}\right)\);
Step 4. Generate σ^{(q+1)} from the conditional distribution \(\pi \left (\sigma \boldmath {p}^{(q)},\alpha ^{(q+1)}, \mu ^{(q+1)},D_{\phantom {\dot {i}\!}S_{j}}\right)\);
Step 5. Simulate \(p_{1}^{(q+1)}\) from the conditional distribution: \(\pi _{\phantom {\dot {i}\!}S_{j}}\left (p_{1}\alpha ^{(q+1)},\mu ^{(q+1)}, \sigma ^{(q+1)}, p_{2}^{(q)},\ldots,p_{\phantom {\dot {i}\!}S_{j}}^{(q)},D_{\phantom {\dot {i}\!}S_{j}}\right)\).
Step 6. Generate \(p_{2}^{(q+1)}\) from the conditional distribution:
StepS_{j}+4. Draw \(p_{\phantom {\dot {i}\!}S_{j}}^{(q+1)}\) from the conditional distribution: \(\pi _{\phantom {\dot {i}\!}S_{j}}\!\left (p_{\phantom {\dot {i}\!}S_{j}}\alpha ^{(q+1)}, \!\mu ^{(q+1)},\!\sigma ^{(q+1)}, \!p_{1}^{(q+1)},\!p_{2}^{(q+1)},\ldots, \!p_{\phantom {\dot {i}\!}S_{j}1}^{(q+1)},D_{\phantom {\dot {i}\!}S_{j}}\right)\).
StepS_{j}+5. Repeat Step 2 to Step S_{j}+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_{i} given \((\boldmath {p}_{i},D_{\phantom {\dot {i}\!}S_{j}})\) 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 \(D_{\phantom {\dot {i}\!}n_{0}}=\{x_{1}< x_{2}<\cdots < x_{\phantom {\dot {i}\!}n_{0}}\}\) denote a set of dose levels in the interval [p_{i−1},p_{i+1}]. For 1≤i≤j≤n_{0}, let l_{ij}(x) be the straight line passing through two points (x_{i}, logπ(x_{i})) and (x_{j}, logπ(x_{j})); otherwise, it is assumed that l_{ij}(x) is not defined. We construct the following piecewise function h_{i}(x):
When l_{1}(x) is undefined but l_{2}(x) is defined above, we set max(l_{1},l_{2})= max(l_{2},l_{1})= min(l_{1},l_{2})= min(l_{2},l_{1})=l_{2}. When both l_{1}(x) and l_{2}(x) are undefined, we take max(l_{1},l_{2})= max(l_{2},l_{1})= min(l_{1},l_{2})= min(l_{2},l_{1})=0. Under the above assumption, the proposal distribution g_{i}(x) is defined as g_{i}(x)= exp(h_{i}(x))/ω, which is a piecewise exponential distribution, where \(\omega =\int \exp (h(x))dx\). Let \(p_{i}^{(q)}\) be the simulated observation of p_{i} at the qth iteration. Then, the ARMS algorithm for drawing \(p_{i}^{(q+1)}\) from the conditional distribution (8) at the (q+1)th iteration is as follows.
Step 1. Initialize n_{0}, which is the number of points in the interval [p_{i−1},p_{i+1}], and determine the set of dose levels \(D_{\phantom {\dot {i}\!}n_{0}}\).
Step 2. Generate \(p_{i}^{*}\) from the proposal distribution g_{i}(x).
Step 3. Generate u from the uniform distribution U(0,1).
Step 4. If \(u>\pi (p_{i}^{*})/\exp (h_{i}(p_{i}^{*}))\), we set \(D_{\phantom {\dot {i}\!}n_{0}+1}=D_{\phantom {\dot {i}\!}n_{0}}\cup \left \{p_{i}^{*}\right \}\), and relabel the points in \(D_{\phantom {\dot {i}\!}n_{0}+1}\) in ascending order and let n_{0}=n_{0}+1, then go to step 2; otherwise, we set \(p_{i}^{*(q+1)}=p_{i}^{*}\).
Step 5. Generate u^{∗} from the uniform distribution U(0,1).
Step 6. If \(u^{*}\!>\!\min \!\left \{\!\!1,\!\frac {\pi \!\left (p_{i}^{*(q+1)}\right)\min \left \{\pi \left (p_{i}^{(q)}\!\right),\, \exp \left (h\left (p_{i}^{(q)}\right)\!\right)\right \}}{\pi \!\left (p_{i}^{(q)}\right)\!\min \left \{\pi \!\left (p_{i}^{*(q+1)}\!\right), \,\exp \left (\!h\left (p_{i}^{\phantom {\dot {i}\!}*(q+1)}\right)\!\right)\right \}}\!\right \}\), we take \(p_{i}^{(q+1)}=p_{i}^{(q)}\); otherwise, we set \(p_{i}^{(q+1)}=p_{i}^{\phantom {\dot {i}\!}*(q+1)}\).
Dosefinding algorithm
In this section, we develop a twostage Bayesian nonparametric adaptive dosefinding design based on a twostage procedure and the above proposed hybrid algorithm. Let c_{e} and c_{d} be the threshold values for dose escalation and deescalation, respectively. In our numerical illustration, c_{e} and c_{d} satisfying the restriction that c_{e}+c_{d}>1 can be selected by the datadependent 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 deescalating for the next cohort of patients to only one dose level of change at a time. The twostage Bayesian nonparametric adaptive dosefinding design is described as follows.
(I) The startup stage
Patients in the first cohort are administered to the lowest dose level d_{1}. 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 d_{2}. This process is continued until at least one toxicity is observed.
(II) The second stage
For dose level d_{k} administered to the jth cohort of patients, if there is at least one toxicity outcome observed, we denote \(d_{k}^{(j)}\) 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 \(d_{k}^{(j+1)}\). Based on the data of the first j+1 cohorts of patients, we can simultaneously obtain \(\hat {\boldmath {p}}=(\hat {p}_{1},\ldots,\hat {p}_{\phantom {\dot {i}\!}S_{j+1}})\) and the toxicity probability \(\text {Pr}(\hat {p}_{k}<\theta)\) at the current dose level \(d_{k}^{(j+1)}\) via the above proposed hybrid algorithm.
(ii) If the probability \(\text {Pr}(\hat {p}_{k}<\theta)>c_{e}\), the dose level administered to the (j+2)th cohort of patients is escalated to the dose level d_{k+1}, i.e. \(d_{k+1}^{(j+2)}=d_{k+1}\). If the current dose \(d_{k}^{(j+1)}=d_{K}\), the dose administered to the (j+2)th cohort of patients is still d_{K}, that is, \(d_{K}^{(j+2)}=d_{K}\).
(iii) If the probability \(\text {Pr}(\hat {p}_{k}<\theta)< c_{d}\), the dose level administered to the (j+2)th cohort of patients is deescalated to the dose level d_{k−1}, i.e. \(d_{k1}^{(j+2)}=d_{k1}\). If the current dose \(d_{k}^{(j+1)}=d_{1}\), the dose level administered to the (j+2)th cohort of patients is still d_{1}, that is, \(d_{1}^{(j+2)}=d_{1}\).
(iv) Otherwise, the (j+2)th cohort of patients continues to be treated at the dose level d_{k} i.e., \(d_{K}^{(j+2)}=d_{K}\).
(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}=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 burnin iterations.
For comparison, we consider two Bayesian modelbased CRM dosefinding methods including the power model: \(p_{i}=d_{i}^{\exp (\beta)}\) and oneparameter logistic model: p_{i}= exp(a_{0}+βd_{i})/{1+ exp(a_{0}+βd_{i})} with a_{0}=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 datadependent 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 F_{0} 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 F_{0}. As mentioned above, we choose appropriate values of μ and σ such that median of F_{0} should be consistent with the initial guess for the MTD. To wit, if there is a prior belief that dose level d_{k} is the MTD, we can select appropriate values of μ and σ such that \(F_{0}(d_{k}\eta)=\Phi \left (\frac {d_{k}\mu }{\sigma }\right)=\theta \). Table 1 gives the values of F_{0}’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_{k}’s, and the preceding proposed twostage Bayesian nonparametric adaptive dosefinding 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_{i} (the dose administered to the current cohort of patients); y_{i1}, y_{i2}, y_{i3} (the observations for the current experiment cohort); \(\hat {p}_{k}\) (estimate of toxicity probability for k=1,…,K). Examination of Table 2 shows that (i) dose level 5 (i.e., d_{5}) is selected as the MTD and \(\hat {p}_{5}=0.2977\), which is quite close to true toxicity probability 0.3; (ii) among 11 cohorts of patients administered to dose level d_{5}, only 10 patients (i.e., \(\sum \nolimits _{i=1}^{20}\sum \nolimits _{m=1}^{3}y_{{im}}\)) are experienced toxicity; (iii) the doses with the relatively high toxicity probability (such as d_{7} and d_{8}) may have no chance to be administered to patients, which guarantees the safety of patients.
Figures 1 and 2 plot Bayesian nonparametric estimations of unknown dosetoxicity curve for three specified values of α for scenarios 1 and 3, respectively. In each figure, “ −⋆” represents the dosetoxicity curve corresponding to true dosetoxicity data, “ −∘” corresponds to the base curve, and “\(\vartriangle \)”, “\(\square \)” and “+” correspond to the estimated dosetoxicity curves for α=5, 10 and 20, respectively. Examination of Figs. 1 and 2 show that the estimated curves are more and more close to F_{0} 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 F_{0}.
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 overMTD 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.
Now we assume that α and F_{0} 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_{k}’s, and the preceding developed twostage Bayesian nonparametric adaptive dosefinding 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.
Results
According to the above presented simulation study, we have the following results. First, the estimated MTD in singleagent dosefinding clinical trials via our proposed twostage Bayesian nonparametric adaptive dosefinding 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 modelbased 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 singleagent dosefinding trials, and the proposed method outperforms two classical CRMs under our considered scenarios.
Discussion
Although this manuscript only considers a single agent dosefinding design, the proposed Bayesian NCRM can be extended to twoagent dosefinding 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 dosefinding 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 singleagent in phase I clinical trials. We relax the traditional parametric model assumption imposed on dosetoxicity relationship using a DP prior to approximate unknown distribution of dosetoxicity 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 twostage Bayesian nonparametric adaptive dosefinding design is developed to estimate the MTD. In the proposed Bayesian nonparametric adaptive dosefinding 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 dosefinding procedure is modelfree and robust, and performs better than two parametric models even in small sample sizes.
Abbreviations
 ARMS:

Adaptive rejection Metropolis sampling
 CRM:

Continual reassessment method
 DP:

Dirichlet process
 DLT:

Doselimiting toxicity
 MTD:

Maximum tolerated dose
 NCRM:

Nonparametric Continual reassessment method
References
 1
O’Quigley J, Pepe M, Fisher L. Continual reassessment method: a practical design for phase 1 clinical trials in cancer. Biometrics. 1990; 46:33–48.
 2
Whitehead J, Brunier H. Bayesian decision procedures for dose determining experiments. Stat Med. 1995; 14:885–93.
 3
Piantadosi S, Fisher JD, Grossman S. Practical implementation of a modified continual reassessment method for dosefinding trials. Cancer Chemother Pharmacol. 1998; 41:429–36.
 4
Heyd JM, Carlin BP. Adaptive design improvements in the continual reassessment method for phase I studies. Stat Med. 1999; 18:1307–21.
 5
Leung DH, Wang YG. An extension of the continual reassessment method using decision theory. Stat Med. 2002; 21:51–63.
 6
Yuan Z, Chappell R, Bailey H. The Continual Reassessment Method for Multiple Toxicity Grades: A Bayesian Quasi ikelihood Approach. Biometrics. 2007; 63:173–9.
 7
Yin G, Yuan Y. Bayesian Model Averaging Continual Reassessment Method in Phase I Clinical Trials. J Am Stat Assoc. 2009; 104:954–68.
 8
Møller S. An extension of the continual reassessment methods using a preliminary upanddown design in a dose finding study in cancer patients, in order to investigate a greater range of doses. Stat Med. 2010; 14:911–22.
 9
Fan SK, Lu Y, Wang YG. A simple Bayesian decisiontheoretic design for dosefinding trials. Stat Med. 2012; 31:3719–30.
 10
Morita S, Thall PF, Takeda K. A simulation study of methods for selecting subgroupspecific dosesin phase i trials. Pharm Stat. 2017; 16:143–56.
 11
Gelfand AE, Kuo L. Nonparametric bayesian bioassay including ordered polytomous response. Biometrika. 1991; 78:657–66.
 12
Mukhopadhyay S. Bayesian nonparametric inference on the dose level with specified response rate. Biometrics. 2000; 56:220–6.
 13
Gasparini M, Eisele J. A curvefree method for phase i clinical trials. Biometrics. 2000; 56:609–15.
 14
Cheng YK. Dose Finding by the Continual Reassessment Method. Boca Raton: Chapman and Hall/CRC; 2011, pp. 57–62.
 15
Ivanova A, Wang K. Bivariate isotonic design for dosefinding with ordered groups. Stat Med. 2006; 25:2018–26.
 16
Yan F, Mandrekar SJ, Yuan Y. Keyboard: a novel bayesian toxicity probability interval design for phase i clinical trials. Clin Cancer Res Off J Am Assoc Cancer Res. 2017; 23:3994–4003.
 17
Lee SM, Cheung YK. Modal calibration in the continual reassessment method. Clin Trials. 2009; 6:227–38.
 18
Ramsey FL. A Bayesian approach to bioassay. Biometrics. 1972; 28:841–58.
 19
Escobar MD, West M. Bayesian density estimation and inference using mixtures. Publ Am Stat Assoc. 1995; 90:577–88.
 20
West M. Hierarchial mixture models in neurological transmission analysis. J Am Stat Assoc. 1997; 92:587–608.
Acknowledgements
The authors are grateful to the Editor, an Associate Editor and two referees for their valuable suggestions and comments for largely improving the manuscript.
Funding
This work was supported by the grants from the National Natural Science Foundation of China (Grant No.: 11671349) and the Key Projects of the National Natural Science Foundation of China (Grant No.: 11731101).
Availability of data and materials
The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.
Author information
Affiliations
Contributions
NST conceived of research questions, developed methods and revised the manuscript; SJW carried out statistical analysis and drafted the manuscript; GY conducted simulation studies. All authors commented on successive drafts, and read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Tang, N., Wang, S. & Ye, G. A nonparametric Bayesian continual reassessment method in singleagent dosefinding studies. BMC Med Res Methodol 18, 172 (2018). https://doi.org/10.1186/s1287401806049
Received:
Accepted:
Published:
Keywords
 Adaptive rejection Metropolis sampling algorithm
 Continual reassessment method
 Dirichlet process prior
 Dosefinding design
 Gibbs sampler
 Maximum tolerated dose