R/get.oc.R

In BOIN: Bayesian Optimal INterval (BOIN) Design for Single-Agent and Drug- Combination Phase I Clinical Trials

#'
#' Generate operating characteristics for single agent trials
#'
#' Obtain the operating characteristics of the BOIN design for single agent trials by simulating trials.
#'
#' @usage get.oc(target, p.true, ncohort, cohortsize, n.earlystop=100,
#'               startdose=1, titration=FALSE, p.saf=0.6*target, p.tox=1.4*target,
#'               cutoff.eli=0.95,extrasafe=FALSE, offset=0.05, boundMTD=FALSE,
#'               ntrial=1000, seed=6)
#'
#' @param target the target DLT rate
#' @param p.true a vector containing the true toxicity probabilities of the
#'              investigational dose levels.
#' @param ncohort the total number of cohorts
#' @param cohortsize the cohort size
#' @param n.earlystop the early stopping parameter. If the number of patients
#'                    treated at the current dose reaches \code{n.earlystop},
#'                    stop the trial and select the MTD based on the observed data.
#'                    The default value \code{n.earlystop=100} essentially turns
#'                    off this type of early stopping.
#' @param startdose the starting dose level for the trial
#' @param titration set \code{titration=TRUE} to perform dose escalation with cohort size = 1 to accelerate dose escalation at the begining of the trial.
#' @param p.saf the highest toxicity probability that is deemed subtherapeutic
#'              (i.e. below the MTD) such that dose escalation should be undertaken.
#'              The default value is \code{p.saf=0.6*target}.
#' @param p.tox the lowest toxicity probability that is deemed overly toxic such
#'              that deescalation is required. The default value is
#'              \code{p.tox=1.4*target}).
#' @param cutoff.eli the cutoff to eliminate an overly toxic dose for safety.
#'                  We recommend the default value of (\code{cutoff.eli=0.95}) for general use.
#' @param extrasafe set \code{extrasafe=TRUE} to impose a more stringent stopping rule
#' @param offset a small positive number (between \code{0} and \code{0.5}) to control how strict the
#'               stopping rule is when \code{extrasafe=TRUE}. A larger value leads to a more
#'               strict stopping rule. The default value \code{offset=0.05} generally works well.
#' @param boundMTD set \code{boundMTD=TRUE} to impose the condition: the isotonic estimate of toxicity probability
#'                 for the selected MTD must be less than de-escalation boundary.
#' @param ntrial the total number of trials to be simulated
#' @param seed the random seed for simulation
#'
#' @details The operating characteristics of the BOIN design are generated by simulating trials
#'          under the prespecified true toxicity probabilities of the investigational doses. If
#'          \code{titration=TRUE}, we perform dose escalation with cohort size = 1 at the begining of the trial:
#'          starting from \code{startdose}, if no toxicity is observed, we escalate the dose;
#'          otherwise, the titration is completed and we switch to cohort size = \code{cohortsize}.
#'          Titration accelerates the dose escalation and is useful when low doses are believed to be safe.
#'
#'
#'          The BOIN design has two built-in stopping rules: (1) stop the trial if the lowest
#'          dose is eliminated due to toxicity, and no dose should be selected as the MTD; and
#'          (2) stop the trial and select the MTD if the number of patients treated at the current
#'          dose reaches \code{n.earlystop}. The first stopping rule is a safety rule to protect patients
#'          from the case in which all doses are overly toxic. The rationale for the second stopping
#'          rule is that when there is a large number (i.e., \code{n.earlystop}) of patients
#'          assigned to a dose, it means that the dose-finding algorithm has approximately converged.
#'          Thus, we can stop the trial early and select the MTD to save sample size and reduce the
#'          trial duration. For some applications, investigators may prefer a more strict safety
#'          stopping rule than rule (1) for extra safety when the lowest dose is overly toxic.
#'          This can be achieved by setting \code{extrasafe=TRUE}, which imposes the following more
#'          strict safety stopping rule: stop the trial if (i) the number of patients treated at the
#'          lowest dose \code{>=3}, and (ii) \eqn{Pr(toxicity\ rate\ of\ the\ lowest\ dose > \code{target} | data)
#'          > \code{cutoff.eli}-\code{offset}}. As a tradeoff, the strong stopping rule will decrease the MTD
#'          selection percentage when the lowest dose actually is the MTD.
#'
#' @return \code{get.oc()} returns the operating characteristics of the BOIN design as a list,
#'        including:
#'        (1) selection percentage at each dose level (\code{$selpercent}),
#'        (2) the number of patients treated at each dose level (\code{$npatients}),
#'        (3) the number of toxicities observed at each dose level (\code{$ntox}),
#'        (4) the average number of toxicities (\code{$totaltox}),
#'        (5) the average number of patients (\code{$totaln}),
#'        (6) the percentage of early stopping without selecting the MTD (\code{$percentstop}),
#'        (7) risk of overdosing 60\% or more of patients (\code{$overdose60}),
#'        (8) risk of overdosing 80\% or more of patients (\code{$overdose80}),
#'        (9) data.frame (\code{$simu.setup}) containing simulation parameters, such as target, p.true, etc.
#'
#' @note We should avoid setting the values of \code{p.saf} and \code{p.tox} very close to the
#'      \code{target}. This is because the small sample sizes of typical phase I trials prevent us from
#'       differentiating the target DLT rate from the rates close to it. The default values provided by
#'       \code{get.oc()} are strongly recommended, and generally yield excellent operating characteristics.
#'
#' @author Suyu Liu, Yanhong Zhou, and Ying Yuan
#'
#' @references Liu S. and Yuan, Y. (2015). Bayesian Optimal Interval Designs for Phase I
#'             Clinical Trials, \emph{Journal of the Royal Statistical Society: Series C}, 64, 507-523.
#'
#'            Yan, F., Zhang, L., Zhou, Y., Pan, H., Liu, S. and Yuan, Y. (2020).BOIN: An R Package
#'            for Designing Single-Agent and Drug-Combination Dose-Finding Trials Using Bayesian Optimal
#'            Interval Designs. \emph{Journal of Statistical Software}, 94(13),1-32.<doi:10.18637/jss.v094.i13>.
#'
#'
#'
#'             Yuan Y., Hess K.R., Hilsenbeck S.G. and Gilbert M.R. (2016) Bayesian Optimal Interval Design: A
#'        Simple and Well-performing Design for Phase I Oncology Trials, \emph{Clinical Cancer Research}, 22, 4291-4301.
#'
#'
#' @seealso Tutorial: \url{http://odin.mdacc.tmc.edu/~yyuan/Software/BOIN/BOIN2.6_tutorial.pdf}
#'
#'          Paper: \url{http://odin.mdacc.tmc.edu/~yyuan/Software/BOIN/paper.pdf}
#'
#' @examples
#'
#' ## get the operating characteristics for BOIN single agent trial
#' oc <- get.oc(target=0.3, p.true=c(0.05, 0.15, 0.3, 0.45, 0.6),
#' 			ncohort=20, cohortsize=3, ntrial=1000)
#'
#' summary(oc) # summarize design operating characteristics
#' plot(oc)   # plot flowchart of the BOIN design and design operating characteristics, including
#'            # selection percentage, number of patients, and observed toxicities at each dose
#'
#'
#' ## perform titration at the begining of the trial to accelerate dose escalation
#' oc <- get.oc(target=0.3, p.true=c(0.05, 0.15, 0.3, 0.45, 0.6),
#' 			titration=TRUE, ncohort=20, cohortsize=3, ntrial=1000)
#'
#' summary(oc)          # summarize design operating characteristics
#' plot(oc)  # plot flowchart of the BOIN design and design operating characteristics
#' @export
get.oc <- function (target, p.true, ncohort, cohortsize, n.earlystop = 100,
                    startdose = 1, titration = FALSE, p.saf = 0.6 * target, p.tox = 1.4 *
                      target, cutoff.eli = 0.95, extrasafe = FALSE, offset = 0.05,boundMTD=FALSE,
                    ntrial = 1000, seed = 6)
{
  if (target < 0.05) {
    stop("the target is too low!")

  }
  if (target > 0.6) {
    stop("the target is too high!")

  }
  if ((target - p.saf) < (0.1 * target)) {
    stop("the probability deemed safe cannot be higher than or too close to the target!")
  }
  if ((p.tox - target) < (0.1 * target)) {
    stop("the probability deemed toxic cannot be lower than or too close to the target!")
  }
  if (offset >= 0.5) {
    stop("the offset is too large!")
  }
  if (n.earlystop <= 6) {
    warning("the value of n.earlystop is too low to ensure good operating characteristics. Recommend n.earlystop = 9 to 18.")
  }
  set.seed(seed)
  if (cohortsize == 1)
    titration = FALSE
  lambda_e = log((1 - p.saf)/(1 - target))/log(target * (1 -
                                                           p.saf)/(p.saf * (1 - target)))
  lambda_d = log((1 - target)/(1 - p.tox))/log(p.tox * (1 -
                                                          target)/(target * (1 - p.tox)))
  ndose = length(p.true)
  npts = ncohort * cohortsize
  Y = matrix(rep(0, ndose * ntrial), ncol = ndose)
  N = matrix(rep(0, ndose * ntrial), ncol = ndose)
  dselect = rep(0, ntrial)

  if (cohortsize > 1) {
    temp = get.boundary(target, ncohort, cohortsize, n.earlystop=ncohort*cohortsize,
                        p.saf, p.tox, cutoff.eli, extrasafe)$full_boundary_tab
  }
  else {
    temp = get.boundary(target, ncohort, cohortsize, n.earlystop=ncohort*cohortsize,
                        p.saf, p.tox, cutoff.eli, extrasafe)$boundary_tab
  }
  b.e = temp[2, ]
  b.d = temp[3, ]
  b.elim = temp[4, ]
  for (trial in 1:ntrial) {
    y <- rep(0, ndose)
    n <- rep(0, ndose)
    earlystop = 0
    d = startdose
    elimi = rep(0, ndose)
    ft=TRUE #flag used to determine whether or not to add cohortsize-1 patients to a dose for the first time when titration is triggered.
    if (titration) {
      z <- (runif(ndose) < p.true)
      if (sum(z) == 0) {
        d = ndose
        n[1:ndose] = 1
      }
      else {
        d = which(z == 1)[1]
        n[1:d] = 1
        y[d] = 1
      }
    }
    for (i in 1:ncohort) {
      if (titration && n[d] < cohortsize && ft){
        ft=FALSE
        y[d] = y[d] + sum(runif(cohortsize - 1) < p.true[d])
        n[d] = n[d] + cohortsize - 1
      }
      else {
           newcohort = runif(cohortsize)<p.true[d];
              if((sum(n)+cohortsize) >= npts){
                nremain = npts - sum(n);
                y[d] = y[d] + sum(newcohort[1:nremain]);
                n[d] = n[d] + nremain;
                break;
            }
              else{
                  y[d] = y[d] + sum(newcohort);
                  n[d] = n[d] + cohortsize;
              }
      }


      if (!is.na(b.elim[n[d]])) {
        if (y[d] >= b.elim[n[d]]) {
          elimi[d:ndose] = 1
          if (d == 1) {
            earlystop = 1
            break
          }
        }
        if (extrasafe) {
          if (d == 1 && n[1] >= 3) {
            if (1 - pbeta(target, y[1] + 1, n[1] - y[1] +
                          1) > cutoff.eli - offset) {
              earlystop = 1
              break
            }
          }
        }
      }


      if(n[d]>=n.earlystop &&
         (
           (y[d]>b.e[n[d]] && y[d]<b.d[n[d]])||
           (d==1 && y[d]>=b.d[n[d]]) ||
           ((d==ndose||elimi[d+1]==1) && y[d]<=b.e[n[d]])
         )
      ) break;

      if (y[d] <= b.e[n[d]] && d != ndose) {
        if (elimi[d + 1] == 0)
          d = d + 1
      }
      else if (y[d] >= b.d[n[d]] && d != 1) {
        d = d - 1
      }
      else {
        d = d
      }
    }
    Y[trial, ] = y
    N[trial, ] = n
    if (earlystop == 1) {
      dselect[trial] = 99
    }
    else {
      dselect[trial] = select.mtd(target, n, y, cutoff.eli,
                                  extrasafe, offset, boundMTD = boundMTD, p.tox=p.tox)$MTD
    }
  }
  selpercent = rep(0, ndose)
  nptsdose = apply(N, 2, mean)
  ntoxdose = apply(Y, 2, mean)
  for (i in 1:ndose) {
    selpercent[i] = sum(dselect == i)/ntrial * 100
  }
  if (length(which(p.true == target)) > 0) {
    if (which(p.true == target) == ndose - 1) {
      overdosing60 = mean(N[, p.true > target] > 0.6 *
                            npts) * 100
      overdosing80 = mean(N[, p.true > target] > 0.8 *
                            npts) * 100
    }
    else {
      overdosing60 = mean(rowSums(N[, p.true > target]) >
                            0.6 * npts) * 100
      overdosing80 = mean(rowSums(N[, p.true > target]) >
                            0.8 * npts) * 100
    }
    out = list(selpercent = selpercent, npatients = nptsdose,
               ntox = ntoxdose, totaltox = sum(Y)/ntrial, totaln = sum(N)/ntrial,
               percentstop = sum(dselect == 99)/ntrial * 100, overdose60 = overdosing60,
               overdose80 = overdosing80, simu.setup = data.frame(target = target,
                                                                  p.true = p.true, ncohort = ncohort, cohortsize = cohortsize,
                                                                  startdose = startdose, p.saf = p.saf, p.tox = p.tox,
                                                                  cutoff.eli = cutoff.eli, extrasafe = extrasafe,
                                                                  offset = offset, ntrial = ntrial, dose = 1:ndose),
               flowchart = TRUE, lambda_e = lambda_e, lambda_d = lambda_d)
  }
  else {
    out = list(selpercent = selpercent, npatients = nptsdose,
               ntox = ntoxdose, totaltox = sum(Y)/ntrial, totaln = sum(N)/ntrial,
               percentstop = sum(dselect == 99)/ntrial * 100, simu.setup = data.frame(target = target,
                                                                                      p.true = p.true, ncohort = ncohort, cohortsize = cohortsize,
                                                                                      startdose = startdose, p.saf = p.saf, p.tox = p.tox,
                                                                                      cutoff.eli = cutoff.eli, extrasafe = extrasafe,
                                                                                      offset = offset, ntrial = ntrial, dose = 1:ndose),
               flowchart = TRUE, lambda_e = lambda_e, lambda_d = lambda_d)
  }
  class(out)<-"boin"
  return(out)
}