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

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


Background
Dose-finding designs for phase I trials have been widely discussed over the past two decades. Many methods have been proposed to identify maximum tolerated dose (MTD) in single-agent dose-finding clinical trials. There are two major branches: model-based and algorithmbased 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 productof-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 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 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 n j = n

The continual reassessment method (CRM)
Let p 01 < . . . < p 0K be an initial guess of toxicity probabilities associated with dose levels 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 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 β|D j is the expectation taken with respect to the posterior probability density function The (j + 1)th cohort of patients is assigned to the dose level where θ is the target toxicity probability. Once the maximum sample size N is reached, the dose with the toxicity probability closest to θ is selected as the MTD.

Nonparametric probability model
The model-based approaches such as CRM and the modified CRMs are widely used in single-agent dose-finding clinical trials. However, parametric modeling of F(·) may be problematical when there is little prior information on the shape of the dose-toxicity curve [12]. Moreover, in some experimental situations, it is recognized that parametric models do not work [18]. To address these issues, we assume that F(·) is an arbitrary right continuous nondecreasing unknown function whose range is [ 0, 1], and can be approximated using some nonparametric method. To this end, a Dirichlet process (DP) prior is adopted to approximate F(·). To wit, we assume F(d 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 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 , 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 (−∞, ∞): Then, the joint prior density π(p) of p = (p 1 , . . . , p K ) can be written as where

Conditional distributions
Under the above assumption, the likelihood function of the first j (j = 1, . . . , J) cohorts of patients (i.e., the dataset D S j = {(n 1 , y 1 ), . . . , (n S j , y S j )}) is given by It follows from Eq. (3) that the joint prior density π(p) of p = (p 1 , . . . , p S j ) can be expressed as where It follows from Eqs. (4) and (5) that the joint posterior probability density of p given the dataset D S j can be expressed as Again, it follows from Eqs. (2) and (3) that the joint conditional distribution of {α, μ, σ } given {D S j , p} is given by which indicates that the conditional distributions of α, μ and σ given {p, D 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 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 S j +1 = 1. The Gibbs sampler for drawing observations α, μ, σ , p 1 , . . . , p S j is implemented as follows. Step S j , and let p 0 = 0 and p S j +1 = 1.
Step S j + 4. Draw p (q+1) S j from the conditional distribution: Step S 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 (p −i , D 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 n 0 = {x 1 < x 2 < · · · < x 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 i be the simulated observation of p i at the qth iteration. Then, the ARMS algorithm for drawing p (q+1) i 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 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 i , and relabel the points in D n 0 +1 in ascending order and let n 0 = n 0 + 1, then go to step 2; otherwise, we set p * (q+1) i = p * i .

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 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 (j) k 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 (j+1) k . Based on the data of the first j + 1 cohorts of patients, we can simultaneously obtainp = (p 1 , . . . ,p S j+1 ) and the toxicity probability Pr(p k < θ) at the current dose level d (j+1) k via the above proposed hybrid algorithm.
(ii) If the probability Pr(p k < θ) > 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 (iii) If the probability Pr(p k < θ) < 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 (iv) Otherwise, the (j + 2)th cohort of patients continues to be treated at the dose level d k i.e., d (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.
For comparison, we consider two Bayesian model-based CRM dose-finding methods including the power model: and one-parameter 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 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 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    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 two-stage 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);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 andp 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., 20 i=1 3 m=1 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 dose-toxicity curve for three specified values of α for scenarios 1 and 3, respectively. In each figure, "− " represents the dose-toxicity curve corresponding to true dose-toxicity data, "−•" corresponds to the base curve, and "− ", "− " and "-+" correspond to the estimated dose-toxicity curves for α = 5, 10 and 20, respectively. Examination of Figs. 1 and 2 show that the estimated curves are more and more close to 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 over-MTD with a relatively lower probability than two parametric CRMs; (iii) Table 3 Selection probabilities and total numbers of toxicities observed for logistic model, power model and NCRM under six scenarios 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 ( (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.

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 singleagent dose-finding trials, and the proposed method outperforms two classical CRMs under our considered scenarios.   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 dosefinding 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.