Run simulations to obtain operating characteristics for methods that include historical data on control

Interim looks are done three times total including the last one, with approximately equal distance between the looks.

Usage

oc2(
  method = c("PointMass.Bayes", "Prior.Bayes", "RCT.Bayes", "RCT.BayesRobust",
    "RCTvanilla.Bayes", "PointMass.PP", "Prior.PP", "RCT.PP", "RCT.PProbust",
    "RCTvanilla.PP"),
  nSim,
  histSize,
  trialSize,
  drift = 0,
  controlRate,
  nmeRate,
  delta,
  tL = 0.8,
  tU = 0.8,
  tT = 0.85,
  phiL = 0.05,
  phiU = 0.95,
  parE = c(1, 1),
  weights
)

Arguments

method: 8 different methods can be selected: PointMass.Bayes: single arm trial with pointmass derived from the historical data on the control response rate, see postprob Prior.Bayes: single arm trial with prior derived from the historical data on the control response rate, see postprobDist RCT.Bayes: RCT (1:1 randomization) with prior derived from the historical data on the control response rate, see postprobDist RCT.BayesRobust: Same as RCT.Bayes, but with a robust prior, enabling dynamic borrowing of the historical information RCT.vanillaBayes: RCT, not using the historical control data PointMass.PP: single arm trial with pointmass derived from the historical data on the control response rate, and using the predictive probability method, see predprob Prior.PP: single arm trial with prior derived from the historical data on the control response rate, and using the predictive probability method, see predprobDist RCT.PP: RCT (1:1 randomization) with prior derived from the historical data on the control response rate, and using the predictive probability method, see predprobDist RCT.PProbust: Same as RCT.PP, but with a robust prior, enabling dynamic borrowing of the historical information RCT.vanillaPP: RCT and using the predictive probability method, not using the historical control data
nSim: number of trials to simulate
histSize: historical data size
trialSize: total trial sample size
drift: drift parameter: the actual control response rate is the historical control response rate plus the drift. Hence a positive drift means that the controls now have higher response rates than in the historical trial.
controlRate: the actual control response rate
nmeRate: the treatment response rate
delta: delta for stopping for efficacy to be used. Implicitly a zero delta for stopping for futility is used.
tL: for Bayes methods: probability threshold for being below control response rate (default: 0.8)
tU: for Bayes methods: probability threshold for being above control response rate + delta (default: 0.8)
tT: threshold for the probability to be above the response rate p0 at the end of the trial (default: 0.85)
phiL: lower threshold on the predictive probability (default: 0.05)
phiU: upper threshold on the predictive probability (default: 0.95)
parE: the beta parameters matrix, with K rows and 2 columns, corresponding to the beta parameters of the K components. Default is a uniform prior.
weights: the mixture weights of the beta mixture prior. Default are uniform weights across mixture components.

Value

Returned operating characteristics in a matrix include: ExpectedN: expected number of patients in the trials ExpectedNactive: expected number of patients with treatment PrStopEarly: probability to stop the trial early (before reaching the maximum sample size) PrEarlyEff: probability to decide for efficacy early PrEarlyFut: probability to decide for futility early PrEfficacy: probability to decide for efficacy PrFutility: probability to decide for futility PrGrayZone: probability of no decision at the end ("gray zone")

Examples

# this looks OK now:
oc2(
  method = "Prior.PP",
  histSize = 1e5L,
  trialSize = 100,
  drift = 0,
  controlRate = 0.5,
  nmeRate = 0.70,
  nSim = 20,
  delta = 0.15
)
#>      ExpectedN ExpectedNactive ExpectedNcontrol PrStopEarly PrEarlyEff
#> [1,]     75.25           75.25                0        0.45        0.2
#>      PrEarlyFut PrEfficacy PrFutility PrGrayZone
#> [1,]       0.25        0.5        0.5          0

## this one as well:
oc2(
  method = "Prior.Bayes",
  histSize = 1e6L,
  trialSize = 100,
  drift = 0,
  controlRate = 0.5,
  nmeRate = 0.65,
  nSim = 20,
  delta = 0.15
)
#>      ExpectedN ExpectedNactive ExpectedNcontrol PrStopEarly PrEarlyEff
#> [1,]       100             100                0           0          0
#>      PrEarlyFut PrEfficacy PrFutility PrGrayZone
#> [1,]          0          0          0          1

## this example takes longer to run, therefore we don't execute it here:
if (FALSE) {
  ## the different methods we are looking at:
  methods <- c(
    "PointMass.Bayes", "Prior.Bayes", "RCT.Bayes", "RCTvanilla.Bayes",
    "PointMass.PP", "Prior.PP", "RCT.PP", "RCTvanilla.PP"
  )

  ## the control rate will always be 50%
  ## but the historical data will be generated randomly from that
  ## binomial distribution!

  ## delta will be 0.15

  ## 200 simulations

  ## the different prior sample sizes (historical data size, histSizes <- c(10, 20, 50, 150))
  ##
  ## later: add very large size, but with heterogeneity =>
  ## mixture of beta priors
  histSizes <- c(20, 150)

  ## the trial sample size (total)
  trialSizes <- c(40)

  ## the drift
  drift <- c(0, 0.1, 0.2)

  ## the different true rates for the NME: nmeRates <- seq(from=0.05, to=0.95, by=0.05)
  ##
  nmeRates <- c(0.4, 0.5, 0.6, 0.7, 0.8)

  ## so the whole grid is
  wholeGrid <- expand.grid(
    method = methods,
    drift = drift,
    histSize = histSizes,
    trialSize = trialSizes,
    nmeRate = nmeRates, stringsAsFactors = FALSE
  )

  ## summarize the methods

  ## compute the operating characteristics for all combinations
  if (file.exists(savefile <- "allSims4.RData")) {
    load(savefile)
  } else {
    ## setup the result list
    allOcs <- vector(
      mode = "list",
      length = nrow(wholeGrid)
    )

    for (i in seq_len(nrow(wholeGrid))) {
      set.seed(i)
      allOcs[[i]] <- with(
        wholeGrid[i, ],
        oc2(
          method = method,
          histSize = histSize,
          trialSize = trialSize,
          drift = drift,
          controlRate = 0.5,
          nmeRate = nmeRate,
          nSim = 200,
          delta = 0.15
        )
      )
    }

    save(allOcs,
      file = savefile
    )
  }

  length(allOcs)
  allOcs

  allOcs[[1]]

  allPossibleNames <- colnames(allOcs[[1]])
  allPossibleNames

  allOcsMatrix <-
    t(sapply(
      allOcs,
      function(x) {
        res <- rep(NA, length(allPossibleNames))
        names(res) <- allPossibleNames
        res[colnames(x)] <- x
        res
      }
    ))


  allOcsMatrix
  wholeMatrix <- cbind(wholeGrid, allOcsMatrix)

  library(reshape)


  ## write a table to a Word file in the current directory
  outputTab <- function(tab,
                        out = "output.docx", # name of the output file
                        # how many digits for numbers?
                        digits = 2) {
    library(officer)
    doc <- read_docx()
    dat <- as.data.frame(tab)
    whichNum <- which(sapply(dat, is.numeric))
    whichInt <- which(sapply(dat, function(x) {
      is.numeric(x) && all(x == round(x))
    }))
    for (i in whichNum) {
      if (i %in% whichInt) {
        dat[[i]] <- format(dat[[i]], nsmall = 0)
      } else {
        dat[[i]] <- format(dat[[i]], digits = digits, nsmall = digits)
      }
    }

    doc <- body_add_table(doc, dat, style = "Normal Table")
    print(doc, target = paste0(getwd(), "/myreport.docx"))
  }


  ## --------------------------------------------------
  ## derive the tables
  ## (here change drift and histSize values manually to get the
  ## different tables for the slides)
  tab1 <- subset(
    wholeMatrix,
    (drift == 0.2) &
      (histSize == 150) &
      (method %in% ppNames)
  )
}