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.

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 d₁<d₂<⋯<d_K 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, 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 dose-limiting toxicity (DLT) among $n_{k}^{j}$ patients treated at dose level d_k. Let S_j be the number of non-zero 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₀₁<…<p_0K be an initial guess of toxicity probabilities associated with dose levels d₁<…<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 dose-toxicity 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₁,…,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₁,…,y_j} be the dataset collected from the first j cohorts of patients. The likelihood function for D_j 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(d_k)∼DP(αF₀(x|η)), where F₀(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₀(x|η) to F(d_k). Thus, the hierarchical model for the dataset {y₁,…,y_K} 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 F₀(x|η) is a prior guess of the unknown distribution F(·). For computational simplicity, F₀(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₀(x|η). More discussions on the selection of F₀(·) see Mukhopadhyay [12]. Here, we take F₀ as the cumulative distribution function of some standard distribution, and choose some appropriate hyperparameters such that the median of F₀ 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₀, and σ is related with the standard variance of F₀. In this case, similar to [12], we assume that the priors of μ and σ discretely follow the uniform distribution in the intervals [μ₀−r,μ₀+r] (e.g., taking ten equally spaced points in [μ₀−r,μ₀+r]) and [σ₀−r,σ₀+r] (e.g., taking ten equally spaced points in [σ₀−r,σ₀+r]) with r∈{1,2}, respectively, where μ₀ and σ₀ are the mean and standard variance of F₀, 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 B₁,…,B_K in the support of F₀, the probability vector (F(B₁),…,F(B_K))^⊤ follows the Dirichlet distribution with parameter vector (αF₀(B₁),…, αF₀(B_K))^⊤. For the base distribution F₀(x|η)=N(x;μ,σ), if we consider the following finite measurable partition for the support (−∞,∞): B₁=(−∞,d₁), B₂=[d₁,d₂),…,B_K=[d_K−1,d_K), B_K+1=[d_K,+∞), thus we have αF₀(B₁)=αF₀(d₁),…,αF₀(B_K+1)=α[1−F₀(d_K)], and F(B₁)=p₁, F(B₂)=p₂−p₁,…,F(B_K+1)=1−p_K. Then, the joint prior density π(p) of p=(p₁,…,p_K)^⊤ 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=α{F₀(d_k)−F₀(d_k−1)} for k=1,…,K+1, d₀=−∞, d_K+1=∞, F₀(−∞)=0, F₀(∞)=1, p₀=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

$$ 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=α{F₀(d_i)−F₀(d_i−1)} for i=1,…,S_j+1, d₀=−∞, $d_{\phantom {\dot {i}\!}S_{j}+1}=\infty$, F₀(−∞)=0, F₀(∞)=1, p₀=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

$$ \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 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

$$ {}\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 p_i−1<p_i<p_i+1 for i=1,…,S_j, and p₀=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 α⁽⁰⁾,μ⁽⁰⁾,σ⁽⁰⁾, $\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 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$$

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₀, 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):

$$\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 l₁(x) is undefined but l₂(x) is defined above, we set max(l₁,l₂)= max(l₂,l₁)= min(l₁,l₂)= min(l₂,l₁)=l₂. When both l₁(x) and l₂(x) are undefined, we take max(l₁,l₂)= max(l₂,l₁)= min(l₁,l₂)= min(l₂,l₁)=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₀, 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₀=n₀+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 c_e and c_d be the threshold values for dose escalation and de-escalation, respectively. In our numerical illustration, c_e and c_d satisfying the restriction that c_e+c_d>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 d₁. 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₂. 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_{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 d₁, 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 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: p_i= exp(a₀+βd_i)/{1+ exp(a₀+βd_i)} with a₀=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⁽⁰⁾ as its corresponding skeleton in implementing Gibbs sampler. When α and F₀ 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₀. As mentioned above, we choose appropriate values of μ and σ such that median of F₀ 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₀’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 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_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₅) 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₅, 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₇ and d₈) may have no chance to be administered to patients, which guarantees the safety of patients.

Table 2 Bayesian estimates of p_k’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 F₀ 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₀.

Fig. 1

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

Fig. 2

Nonparametric Bayesian estimation of dose-toxicity curve under scenarios 3 when α and F₀ 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 F₀ 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 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.