Skip to main content
  • Research article
  • Open access
  • Published:

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

Abstract

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.

Peer Review reports

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 d1<d2<<dK be a set of K pre-specified dose levels for a single agent under investigation, θ be the target toxicity probability, the dose-toxicity curve F(·) be the distribution of dose levels, pk=F(dk) be the toxicity probability at dose level dk 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 dk in the first j cohorts, and \(y_{k}^{j}\) represents the number of patients experienced the dose-limiting toxicity (DLT) among \(n_{k}^{j}\) patients treated at dose level dk. Let Sj be the number of non-zero components in nj. Thus, we have SjK. 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 nj and yj.

The continual reassessment method (CRM)

Let p01<…<p0K be an initial guess of toxicity probabilities associated with dose levels d1<…<dK, where p0k is referred to as the “skeleton” of the CRM. It is assumed that there is a prior belief that dose level dk is the MTD, one may set the initial guess as p0k=θ. 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 dose-toxicity relationship F(dk) is a strictly increasing function with respect to the dose level dk, and is usually specified by a parametric model F(dk;β) in which β is an unknown parameter to be estimated, such as, logistic function, power function, and hyperbolic function.

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

$$L_{j}(\beta|D_{j})\propto\prod_{i=1}^{j}\{F(x_{i};\beta)\}^{\phantom{\dot{i}\!}y_{i}}\{1-F(x_{i};\beta)\}^{\phantom{\dot{i}\!}n_{i}-y_{i}}. $$

Let f(β) be the prior density of β. It follows from the Bayesian theorem that a Bayesian estimator of β is given by

$$\hat{\beta}_{j}=E_{\phantom{\dot{i}\!}\beta|D_{j}}(\beta)=\frac{\int\beta L_{j}(\beta|D_{j})f(\beta)d\beta}{\int L_{j}(\beta|D_{j})f(\beta)d\beta}, $$

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

$$x_{j+1}=\text{arg}\min_{\phantom{\dot{i}\!}d_{k}}|F\left(d_{k};\hat{\beta}_{j}\right)-\theta|, ~k=1,\ldots,K, $$

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(dk)DP(αF0(x|η)), where F0(x|η) is a base distribution that serves as a starting point for constructing a nonparametric distribution for F(dk), η 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(dk). Thus, the hierarchical model for the dataset {y1,…,yK} is given by

$$ \begin{array}{lll} y_{k}|p_{k} \overset{\text{ind}}{\sim}{\text{Bin}(n_{k},p_{k}),\ \ k=1,\ldots,K,}\\ p_{k}=F(d_{k}),\ \ \ F|\eta,\alpha \sim \text{DP}(\alpha F_{0}(x|\eta)). \end{array} $$
(1)

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 dk 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 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 [μ0r,μ0+r] (e.g., taking ten equally spaced points in [μ0r,μ0+r]) and [σ0r,σ0+r] (e.g., taking ten equally spaced points in [σ0r,σ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

$$ \begin{array}{lll} y_{k}|p_{k} \overset{\text{ind}}{\sim} \text{Bin}(n_{k},p_{k}),\ \ k=1,\ldots,K,\\ p_{k}=F(d_{k}), \ \ \ F|\eta,\alpha \sim \text{DP}(\alpha F_{0}(x|\eta)),\\ \alpha \sim \Gamma(a,b), \ \ \ \eta \sim \pi(\mu)\pi(\sigma). \end{array} $$
(2)

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

$$ \pi(\boldmath{p})=\frac{\Gamma\left(\sum\limits_{k=1}^{K+1}\gamma_{k}\right)}{\prod\limits_{k=1}^{K+1} \Gamma(\gamma_{k})}\prod\limits_{k=1}^{K+1}(p_{k}-p_{k-1})^{\phantom{\dot{i}\!}\gamma_{k}-1}, $$
(3)

where γk=α{F0(dk)−F0(dk−1)} for k=1,…,K+1, d0=−, dK+1=, F0(−)=0, F0()=1, p0=0, and pK+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

$$ L_{\phantom{\dot{i}\!}S_{j}}(\boldmath{p}|D_{\phantom{\dot{i}\!}S_{j}})\propto\prod\limits_{i=1}^{S_{j}} p_{i}^{\phantom{\dot{i}\!}y_{i}}(1-p_{i})^{\phantom{\dot{i}\!}n_{i}-y_{i}},~1\leq S_{j}\leq K. $$
(4)

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

$$ \pi(\boldmath{p}) \propto \prod\limits_{i=1}^{S_{j}+1}(p_{i}-p_{i-1})^{\phantom{\dot{i}\!}\gamma_{i}-1},\ \ 1\leq S_{j}\leq K, $$
(5)

where γi=α{F0(di)−F0(di−1)} for i=1,…,Sj+1, d0=−, \(d_{\phantom {\dot {i}\!}S_{j}+1}=\infty \), F0(−)=0, F0()=1, p0=0, and \(p_{\phantom {\dot {i}\!}S_{j}+1}=1\). It is easily seen that Eq. (5) reduces to Eq. (3) when Sj=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

$$ \pi(\boldmath{p}|D_{\phantom{\dot{i}\!}S_{j}})\!\propto\! \prod\limits_{i=1}^{S_{j}} p_{i}^{\phantom{\dot{i}\!}y_{i}}(1-p_{i})^{\phantom{\dot{i}\!}n_{i}-y_{i}}\!\prod\limits_{i=1}^{S_{j}+1}\!\!(p_{i}-p_{i-1})^{\phantom{\dot{i}\!}\gamma_{i}-1},\ \ \!\!1\!\leq\! S_{j}\!\leq \!K. $$
(6)

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

$${}\pi(\alpha,\mu,\sigma|\boldmath{p},D_{\phantom{\dot{i}\!}S_{j}})\propto \frac{\Gamma\left(\sum\limits_{i=1}^{S_{j}+1}\gamma_{i}\right)}{\prod\limits_{i=1}^{S_{j}+1} \Gamma(\gamma_{i})}\prod_{i=1}^{S_{j}+1}(p_{i}-p_{i-1})^{\phantom{\dot{i}\!}\gamma_{i}-1}\pi(\alpha)\pi(\mu)\pi(\sigma), $$

which indicates that the conditional distributions of α, μ and σ given \(\{\boldmath {p},D_{\phantom {\dot {i}\!}S_{j}}\}\) have the following expressions:

$$ \begin{array}{lll} \pi(\alpha|\boldmath{p},\mu,\sigma,D_{\phantom{\dot{i}\!}S_{j}}) &\propto& \frac{\Gamma(\alpha)}{\prod\limits_{i=1}^{S_{j}+1}\Gamma(\gamma_{i})}\prod\limits_{i=1}^{S_{j}+1}(p_{i}-p_{i-1})^{\phantom{\dot{i}\!}\gamma_{i}-1}\alpha^{a-1}e^{-b\alpha}, \\ \pi(\mu,\sigma|\boldmath{p},\alpha,D_{\phantom{\dot{i}\!}S_{j}}) &\propto & \frac{\Gamma(\alpha)}{\prod\limits_{i=1}^{S_{j}+1}\Gamma(\gamma_{i})}\prod\limits_{i=1}^{S_{j}+1}(p_{i}-p_{i-1})^{\phantom{\dot{i}\!}\gamma_{i}-1}\pi(\mu)\pi(\sigma).\\ \end{array} $$
(7)

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 pi be a subset vector of p with the ith element of p deleted. It follows from Eq. (6) that the conditional distribution of pi given pi and \(D_{\phantom {\dot {i}\!}S_{j}}\) can be written as

$$ {}\pi(p_{i}|\boldmath{p}_{-i},D_{\phantom{\dot{i}\!}S_{j}})\!\propto p_{i}^{\phantom{\dot{i}\!}y_{i}}\!(1\!-p_{i})^{\phantom{\dot{i}\!}n_{i}-y_{i}}\!(p_{i}-p_{i-1})^{\phantom{\dot{i}\!}\gamma_{i}-\!1}\!(p_{i+1}\,-\,p_{i})^{\phantom{\dot{i}\!}\gamma_{i+1}-1}, $$
(8)

where pi−1<pi<pi+1 for i=1,…,Sj, and p0=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 p0=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:

$$ {}\pi_{\phantom{\dot{i}\!}S_{j}}\!\!\left(p_{2}|\alpha^{(q+1)}, \mu^{\!(q+1)},\!\sigma^{\!(q+1)}, \!p_{1}^{(q+1)},p_{3}^{(q)},\ldots,p_{\phantom{\dot{i}\!}S_{j}}^{(q)},\!p_{\phantom{\dot{i}\!}S_{j}+\!1},D_{\phantom{\dot{i}\!}S_{j}}\!\right)\!. $$
(9)
$$\hspace{3.5cm}\vdots$$

StepSj+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)\).

StepSj+5. Repeat Step 2 to Step Sj+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 pi 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 [pi−1,pi+1]. For 1≤ijn0, let lij(x) be the straight line passing through two points (xi, logπ(xi)) and (xj, logπ(xj)); otherwise, it is assumed that lij(x) is not defined. We construct the following piecewise function hi(x):

$$\begin{array}{*{20}l} h_{i}(x)=&\max\{l_{i,i+1}(x),\min\{l_{i-1,i}(x),l_{i+1,i+2}(x)\}\}\\ &\;\text{for}\ \ x_{i}\leq x\leq x_{i+1}. \end{array} $$
(10)

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 gi(x) is defined as gi(x)= exp(hi(x))/ω, which is a piecewise exponential distribution, where \(\omega =\int \exp (h(x))dx\). Let \(p_{i}^{(q)}\) be the simulated observation of pi 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 n0, which is the number of points in the interval [pi−1,pi+1], and determine the set of dose levels \(D_{\phantom {\dot {i}\!}n_{0}}\).

Step 2. Generate \(p_{i}^{*}\) from the proposal distribution gi(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 n0=n0+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)}\).

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 ce and cd be the threshold values for dose escalation and de-escalation, respectively. In our numerical illustration, ce and cd satisfying the restriction that ce+cd>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 dk 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 dk+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 dK, 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 dk−1, i.e. \(d_{k-1}^{(j+2)}=d_{k-1}\). If the current dose \(d_{k}^{(j+1)}=d_{1}\), the dose level administered to the (j+2)th cohort of patients is still d1, 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 dk 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, dk=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

For comparison, we consider two Bayesian model-based CRM dose-finding methods including the power model: \(p_{i}=d_{i}^{\exp (\beta)}\) and one-parameter logistic model: pi= exp(a0+βdi)/{1+ exp(a0+βdi)} 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 dk 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 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 pk’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); xi (the dose administered to the current cohort of patients); yi1, yi2, yi3 (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., d5) 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 d5, 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 d7 and d8) may have no chance to be administered to patients, which guarantees the safety of patients.

Table 2 Bayesian estimates of pk’s via NCRM under scenario 1

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 “\(-\vartriangle \)”, “\(-\square \)” 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
figure 1

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

Fig. 2
figure 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
Table 4 Average numbers of patients treated at each of eights doses for logistic model, power model and NCRM

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 pk’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
Table 6 Average numbers of patients treated at each of eight doses for NCRM under scenario 1 and 3 when α and (μ,σ) are unknown

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.

Abbreviations

ARMS:

Adaptive rejection Metropolis sampling

CRM:

Continual reassessment method

DP:

Dirichlet process

DLT:

Dose-limiting 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.

    Article  Google Scholar 

  2. Whitehead J, Brunier H. Bayesian decision procedures for dose determining experiments. Stat Med. 1995; 14:885–93.

    Article  CAS  Google Scholar 

  3. Piantadosi S, Fisher JD, Grossman S. Practical implementation of a modified continual reassessment method for dose-finding trials. Cancer Chemother Pharmacol. 1998; 41:429–36.

    Article  CAS  Google Scholar 

  4. Heyd JM, Carlin BP. Adaptive design improvements in the continual reassessment method for phase I studies. Stat Med. 1999; 18:1307–21.

    Article  CAS  Google Scholar 

  5. Leung DH, Wang YG. An extension of the continual reassessment method using decision theory. Stat Med. 2002; 21:51–63.

    Article  Google Scholar 

  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.

    Article  CAS  Google Scholar 

  7. Yin G, Yuan Y. Bayesian Model Averaging Continual Reassessment Method in Phase I Clinical Trials. J Am Stat Assoc. 2009; 104:954–68.

    Article  CAS  Google Scholar 

  8. Møller S. An extension of the continual reassessment methods using a preliminary up-and-down design in a dose finding study in cancer patients, in order to investigate a greater range of doses. Stat Med. 2010; 14:911–22.

    Article  Google Scholar 

  9. Fan SK, Lu Y, Wang YG. A simple Bayesian decision-theoretic design for dose-finding trials. Stat Med. 2012; 31:3719–30.

    Article  Google Scholar 

  10. Morita S, Thall PF, Takeda K. A simulation study of methods for selecting subgroup-specific dosesin phase i trials. Pharm Stat. 2017; 16:143–56.

    Article  Google Scholar 

  11. Gelfand AE, Kuo L. Nonparametric bayesian bioassay including ordered polytomous response. Biometrika. 1991; 78:657–66.

    Article  Google Scholar 

  12. Mukhopadhyay S. Bayesian nonparametric inference on the dose level with specified response rate. Biometrics. 2000; 56:220–6.

    Article  CAS  Google Scholar 

  13. Gasparini M, Eisele J. A curve-free method for phase i clinical trials. Biometrics. 2000; 56:609–15.

    Article  CAS  Google Scholar 

  14. Cheng YK. Dose Finding by the Continual Reassessment Method. Boca Raton: Chapman and Hall/CRC; 2011, pp. 57–62.

    Book  Google Scholar 

  15. Ivanova A, Wang K. Bivariate isotonic design for dose-finding with ordered groups. Stat Med. 2006; 25:2018–26.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  17. Lee SM, Cheung YK. Modal calibration in the continual reassessment method. Clin Trials. 2009; 6:227–38.

    Article  Google Scholar 

  18. Ramsey FL. A Bayesian approach to bioassay. Biometrics. 1972; 28:841–58.

    Article  CAS  Google Scholar 

  19. Escobar MD, West M. Bayesian density estimation and inference using mixtures. Publ Am Stat Assoc. 1995; 90:577–88.

    Article  Google Scholar 

  20. West M. Hierarchial mixture models in neurological transmission analysis. J Am Stat Assoc. 1997; 92:587–608.

    Article  Google Scholar 

Download references

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

Authors and Affiliations

Authors

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

Correspondence to Niansheng Tang.

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.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Tang, N., Wang, S. & Ye, G. A nonparametric Bayesian continual reassessment method in single-agent dose-finding studies. BMC Med Res Methodol 18, 172 (2018). https://doi.org/10.1186/s12874-018-0604-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12874-018-0604-9

Keywords