Bayesian Continual Reassessment Method for Phase I DoseEscalation Trials
Description
Implements a wide variety of Bayesian CRM designs, including 1parameter, 2parameter and Escalation With Overdose Control (EWOC) designs. The program can run interactively, allowing the user to enter outcomes after each cohort has been recruited, or via simulation to assess operating characteristics.
Usage
1 2 3 4 5 6 7 8  bcrm(stop = list(nmax = NULL, nmtd = NULL, precision = NULL, nmin = NULL),
data = NULL, p.tox0 = NULL, sdose = NULL, dose = NULL,
ff, prior.alpha, cohort = 3, target.tox, constrain = TRUE,
sdose.calculate = "mean", pointest = "plugin", tox.cutpoints = NULL,
loss = NULL, start = NULL, simulate = FALSE, nsims = 1, truep = NULL,
threep3 = FALSE, method = "exact", burnin.itr = 2000, production.itr = 2000,
bugs.directory = "c:/Program Files/WinBUGS14/", plot = FALSE, seed = NULL,
quietly = FALSE, file = NULL, N, tox, notox)

Arguments
stop 
A list of stopping rules for the trial. One or more of the following options should be specified

data 
A named data frame giving information about dose and toxicity from previously recruited patients. If missing, then it is assumed that no data have thus far been collected. Contains the following variables:

p.tox0 
A vector of length 
sdose 
A vector of length 
dose 
Optional vector of length 
ff 
A string indicating the functional form of the doseresponse curve. Options are

prior.alpha 
A list of length 3 containing the distributional information for the prior. The first element is a number from 14 specifying the type of distribution. Options are
The second and third elements of the list are the parameters a and b, respectively. 
cohort 
The size of each cohort of patients that are sequentially recruited to the trial. Defaults to 3 
target.tox 
The target toxicity probability. Defaults to 1/3. 
constrain 
Should a doseskipping constraint be placed on the escalation procedure, as imposed by a modified CRM? Defaults to TRUE. 
sdose.calculate 
What plugin estimate of the prior alpha should be used to calculate the standardised doses? Options are 
pointest 
Which summary estimate of the posterior distribution should be used to choose the next dose. Options are 
tox.cutpoints 
A vector of cutpoints for toxicity intervals if these are to be used to choose next dose. Defaults to NULL. For example Underdosing [0,0.2], Target dosing (0.2, 0.35], Excessive toxicity (0.35, 0.60], Unacceptable toxicity (0.60, 1.00] set 
loss 
A vector of length 
start 
Dose level used at the beginning of the trial. Required if 
simulate 
Should a simulation be conducted to assess operating characteristics? Defaults to 
nsims 
Number of simulations to perform if 
truep 
A vector of length k giving the true probabilities of the outcome (toxicity) at each dose level 
threep3 
Should operating characteristics of a standard 3+3 rulebased design be calculated, for comparison with 
method 
Optimisation method: options are 
burnin.itr 
Number of burnin iterations (default 2000). 
production.itr 
Number of production iterations (default 2000). 
bugs.directory 
Directory that contains the WinBUGS executable if 
plot 
Should the doseresponse curve be plotted after each cohort has been entered? Defaults to FALSE. 
seed 
Integer defining the state of the random number generator to allow reproducibility of results. The default is to not specify a seed. 
quietly 
Should the simulation number output be suppressed when running bcrm? Defaults to FALSE. 
file 
File name where the doseresponse plots are stored, in a pdf format. The program will ammend the current sample size to the end of the file name. 
N 
Final sample size (deprecated). To be replaced with 
tox 
(Deprecated). A vector of length 
notox 
(Deprecated). A vector of length 
Details
bcrm
implements a Bayesian continual reassessment method (CRM) (O'Quigley et al., 1990); an adaptive design in which cohorts of patients are sequentially recruited into a Phase I trial. A binary toxicity outcome is assumed (e.g. Dose Limiting Toxicity / No Dose Limiting Toxicity). The current cohort are given a dose "closest" to the specified target toxicity level, as estimated from the posterior distributions of toxicity at each dose level from the patients thus far recruited. If pointest="mean"
then the posterior mean probability of toxicity is used to choose the next dose. If pointest="plugin"
, however, the posterior mean of the model parameter(s) is pluggedinto the functional form of the dosetoxicity model. To implement an EWOC design (Babb et al., 1998), pointest
should be a quantile, q, between 0 and 0.5. The posterior distribution of the MTD (the dose in which the probability of toxicity is equal to the target toxicity) is then calculated and the next patient is given dose closest to the qth quantile of the MTD distribution.
Alternatively, escalation can be based on intervals of toxicity from the posterior distribution using a loss function, see Neuenschwander et al., 2008. To implement this approach, the user should specify the cutpoints of the toxicity intervals using tox.cutpoints
and the associated losses using loss
.
The possible choice of dosetoxicity model can be specified using ff
, and includes the 1parameter hyberbolic tangent, logistic or power "working models", and the 2parameter logistic model as follows:
 Hyperbolic Tangent
p(Toxd*)=[(tanh(d*)+1)/2]^α
 Logistic (1parameter)
p(Toxd*)=exp(3+α d*)/(1+exp(3+α d*))
 Power
p(Toxd*)=d*^α
 Logistic (2parameter)
p(Toxd*)=exp(log(α_1)+α_2 d*)/(1+exp(log(α_1)+α_2 d*))
where α>0 is the single positivevalued parameter for the 1parameter models, and log(α_1) and α_2>0 are the intercept and slope parameters of the 2parameter model.
The standardised doses, d*, are specified by the user using sdose
, or alternatively the prior probability of toxicity at each dose level is specified using p.tox0
. If the latter is used, then the standardised doses are calculated using the inverse of the functional form and a plugin estimate of the prior mean or median, as specified in sdose.calculate
, as follows
d* = f^{1}(p.tox0,α= a)
where f^{1} is the the inverse of the chosen functional form, and the parameter(s) of the model are set equal to a, either the prior mean or median of α.
Data that have already been accrued can be entered using the data
argument. A constrained CRM design can be implemented using constrain=TRUE
, in which case doseskipping is prohibited (i.e. the next cohort can only be dosed up to one dose level above the current cohort). If a constrained model is used then the starting dose must be specified using start
. Alternatively, if data have already been accrued, then the dose level of the last recruited patient determines the constraint for the next patient.
The prior is set using prior.alpha
. For example prior.alpha=list(1,1,1)
specifies a Gamma prior with shape and scale parameters both equal to one (i.e. an Exponential(1) distribution), whilst prior.alpha=list(2,0,10)
specifies a Uniform(0,10) prior.
To specify a fixed maximum sample size of size m
use
stop=list(nmax=m)
. Alternatively, the trial can stop after m2
patients have been treated at the current MTD estimate, by setting
stop=list(nmtd=m2)
. To stop the trial when the MTD estimate is within a certain level of precision, use stop=list(precision=c(l,u))
, where l
and u
are the lower and upper percentage points that the MTD 95% credible intervals for the risk of toxicity should lie within. Finally, to prevent the trial stopping too early using these rules, the argument stop=list(nmin=m3)
can be used to ensure the sample size is greater than or equal to m3
. Stopping rules can be used on their own or in combination.
The trial can be run interactively using simulate=FALSE
, where the user enters the outcomes for each new cohort, or as a simulation study when simulate=TRUE
.
The default calculations use exact methods (method="exact"
) to calculate the mean and quantiles for the posterior distributions. There are three choices for MCMC calculations:
method="rjags"
, method="BRugs"
or method="R2WinBUGS"
. The first uses the JAGS software, the second uses OpenBUGS, whilst the latter uses WinBUGS. To implement these methods, users require one or more of these packages to be installed on their system.
A simulated bcrm
design can be compared with the standard 3+3 rulebased method, see threep3
for more details.
Value
bcrm
returns an object of class "bcrm" or "bcrm.sim"; the latter occuring when a simulation has been conducted (simulate=TRUE
).
The function print
(i.e. print.bcrm
or print.bcrm.sim
) can be used to obtain summary information about the design used, the data observed, current posterior estimates of toxicity, and the next recommended dose level.
An object of class "bcrm" is a list with the following components:
dose 
Range of doses 
sdose 
Standardised doses 
tox 
A vector of length 
notox 
A vector of length 
ndose 
A list of lists containing for each cohort the components 
constrain 
Whether a constrained CRM design was used 
start 
The starting dose for the latest run of the model if 
target.tox 
The target toxicity level 
ff 
A number from 14 identifying the functional form, 1 = Hyperbolic tangent, 2 = 1parameter logistic, 3 = Power, 4 = 2parameter logistic 
method 
The calculation method used 
pointest 
The summary estimate used to choose the next dose, see 
prior.alpha 
Information about the prior used for alpha, see 
data 
A data frame with variables ‘patient’, ‘dose’ and ‘tox’ listing the doses and outcomes of all patients in the trial 
An object of class "bcrm.sim" is a list of length nsims
. Each component is itself a list with components similar to those obtained from a "bcrm" object. The print function, print.bcrm.sim
should be used to obtain operating characteristics from the simulation.
Note
Currently, the reparameterisation of the twoparameter model proposed by (Babb et al., 1998) is not implemented. Therefore, users wishing to implement an EWOC design should check whether their choice of prior for the model parameter(s) translates to a sensible prior for the MTD distribution before they implement the design. For example
1 2 3 4 5 6 
Oneparameter models are designed as working models only, and should not be used with an escalation strategy based on intervals of the posterior probabilities of toxicity.
Author(s)
Michael Sweeting mjs212@medschl.cam.ac.uk (University of Cambridge, UK), drawing on code originally developed by J. Jack Lee and Nan Chen, Department of Biostatistics, the University of Texas M. D. Anderson Cancer Center
References
Sweeting M., Mander A., Sabin T. bcrm: Bayesian Continual Reassessment Method Designs for Phase I DoseFinding Trials. Journal of Statistical Software (2013) 54: 1–26. http://www.jstatsoft.org/article/view/v054i13
O'Quigley J., Pepe M., Fisher L. Continual reassessment method: a practical design for phase I clinical trials in cancer. Biometrics (1990) 46: 33–48.
Babb J., Rogatko A., Zacks S. Cancer phase I clinical trials: efficient dose escalation with overdose control. Statistics in Medicine (1998) 17: 1103–1120.
Neuenschwander B., Branson M., Gsponer T. Critical aspects of the Bayesian approach to phase I cancer trials. Statistics in Medicine (2008) 27: 2420–2439.
See Also
print.bcrm
, print.bcrm.sim
, plot.bcrm
, plot.bcrm.sim
, threep3
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94  ## Doseescalation cancer trial example as described in Neuenschwander et al 2008.
## Predefined doses
dose<c(1,2.5,5,10,15,20,25,30,40,50,75,100,150,200,250)
## Prespecified probabilities of toxicity
## [dose levels 1115 not specified in the paper, and are for illustration only]
p.tox0<c(0.010,0.015,0.020,0.025,0.030,0.040,0.050,0.100,0.170,0.300,0.400,0.500,0.650
,0.800,0.900)
## Data from the first 5 cohorts of 18 patients
data<data.frame(patient=1:18,dose=rep(c(1:4,7),c(3,4,5,4,2)),tox=rep(0:1,c(16,2)))
## Target toxicity level
target.tox<0.30
## A 1parameter power model is used, with standardised doses calculated using
## the plugin prior median
## Prior for alpha is lognormal with mean 0 (on log scale)
## and standard deviation 1.34 (on log scale)
## The recommended dose for the next cohort if posterior mean is used
Power.LN.bcrm<bcrm(stop=list(nmax=18),data=data,p.tox0=p.tox0,dose=dose
,ff="power",prior.alpha=list(3,0,1.34^2),target.tox=target.tox,constrain=FALSE
,sdose.calculate="median",pointest="mean")
print(Power.LN.bcrm)
plot(Power.LN.bcrm)
## Simulate 10 replicate trials of size 36 (cohort size 3) using this design
## with constraint (i.e. no doseskipping) and starting at lowest dose
## True probabilities of toxicity are set to prespecified probabilities (p.tox0)
Power.LN.bcrm.sim<bcrm(stop=list(nmax=36),p.tox0=p.tox0,dose=dose,ff="power"
,prior.alpha=list(3,0,1.34^2),target.tox=target.tox,constrain=TRUE
,sdose.calculate="median",pointest="mean",start=1,simulate=TRUE,nsims=10,truep=p.tox0)
print(Power.LN.bcrm.sim)
plot(Power.LN.bcrm.sim)
## Comparing this CRM design with the standard 3+3 design
## (only considering the first 12 dose levels)
## Not run:
Power.LN.bcrm.compare.sim<bcrm(stop=list(nmax=36),p.tox0=p.tox0[1:12],dose=dose[1:12]
,ff="power",prior.alpha=list(3,0,1.34^2),target.tox=target.tox,constrain=TRUE
,sdose.calculate="median",pointest="mean",start=1,simulate=TRUE,nsims=50
,truep=p.tox0[1:12],threep3=TRUE)
print(Power.LN.bcrm.compare.sim,threep3=TRUE)
plot(Power.LN.bcrm.compare.sim,threep3=TRUE)
## End(Not run)
## A 2parameter model, using priors as specified in Neuenschwander et al 2008.
## Posterior mean used to choose the next dose
## Standardised doses using reference dose, 250mg
sdose<log(dose/250)
## Bivariate lognormal prior for two parameters
mu<c(2.15,0.52)
Sigma<rbind(c(0.84^2,0.134),c(0.134,0.80^2))
## Using rjags (requires JAGS to be installed)
TwoPLogistic.mean.bcrm<bcrm(stop=list(nmax=18),data=data,sdose=sdose
,dose=dose,ff="logit2",prior.alpha=list(4,mu,Sigma),target.tox=target.tox
,constrain=FALSE,pointest="mean",method="rjags")
print(TwoPLogistic.mean.bcrm)
plot(TwoPLogistic.mean.bcrm)
## A 2parameter model, using an EWOC design with feasibility bound (MTD quantile)
## of 0.25 to choose the next dose
## Using rjags (requires JAGS to be installed)
## Not run:
TwoPLogistic.EWOC0.25.bcrm<bcrm(stop=list(nmax=18),data=data,sdose=sdose,dose=dose
,ff="logit2",prior.alpha=list(4,mu,Sigma),target.tox=target.tox,constrain=FALSE
,pointest=0.25,method="rjags")
print(TwoPLogistic.EWOC0.25.bcrm)
plot(TwoPLogistic.EWOC0.25.bcrm)
## End(Not run)
## A 2parameter model, using a loss function based on intervals of toxicity to choose
## the next dose
## Using rjags (requires JAGS to be installed)
## Not run:
## Toxicity cutpoints
tox.cutpoints<c(0.2,0.35,0.6)
## Losses associated with toxicity intervals
## [0,0.2]=1, (0.2,0.35]=0, (0.35,0.6]=1, (0.6,1]=2
loss<c(1,0,1,2)
TwoPLogistic.tox.intervals.bcrm<bcrm(stop=list(nmax=18),data=data,sdose=sdose
,dose=dose,ff="logit2",prior.alpha=list(4,mu,Sigma),target.tox=target.tox
,constrain=FALSE,tox.cutpoints=tox.cutpoints,loss=loss,method="rjags")
print(TwoPLogistic.tox.intervals.bcrm)
plot(TwoPLogistic.tox.intervals.bcrm)
## Greater loss associated with overdosing and unacceptable toxicity
## [0,0.2]=1, (0.2,0.35]=0, (0.35,0.6]=2, (0.6,1]=4
loss2<c(1,0,2,4)
TwoPLogistic.tox.intervals.2.bcrm<bcrm(stop=list(nmax=18),data=data,sdose=sdose
,dose=dose,ff="logit2",prior.alpha=list(4,mu,Sigma),target.tox=target.tox
,constrain=FALSE,tox.cutpoints=tox.cutpoints,loss=loss2,method="rjags")
print(TwoPLogistic.tox.intervals.2.bcrm)
plot(TwoPLogistic.tox.intervals.2.bcrm)
## End(Not run)
