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

This function differs from oc2 with respect to how the historical and actual control rates are determined, see the parameter controlRateDist.

Usage

oc3(
  method = c("PointMass.Bayes", "Prior.Bayes", "RCT.Bayes", "RCTvanilla.Bayes",
    "PointMass.PP", "Prior.PP", "RCT.PP", "RCTvanilla.PP"),
  nSim,
  histSize,
  trialSize,
  controlRateDist,
  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.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.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
controlRateDist: the control response rate distribution. This shall be a function without any arguments, that gives back a random draw from the control response rate distribution. It is invoked to obtain control rates for the historical and the actual control data sets.
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")

Details

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

Examples

## this looks OK now:
oc3(
  method = "Prior.PP",
  histSize = 1e5L,
  trialSize = 100,
  controlRateDist = function() {
    rbeta(n = 1, 10, 10)
  },
  nmeRate = 0.70,
  nSim = 20,
  delta = 0.15
)
#>      ExpectedN ExpectedNactive ExpectedNcontrol PrStopEarly PrEarlyEff
#> [1,]      60.4            60.4                0         0.7        0.3
#>      PrEarlyFut PrEfficacy PrFutility PrGrayZone
#> [1,]        0.4       0.45       0.55          0

## this one as well:
oc3(
  method = "Prior.Bayes",
  histSize = 1e6L,
  trialSize = 100,
  controlRateDist = function() {
    rbeta(n = 1, 10, 10)
  },
  nmeRate = 0.65,
  nSim = 20,
  delta = 0.15
)
#>      ExpectedN ExpectedNactive ExpectedNcontrol PrStopEarly PrEarlyEff
#> [1,]     91.75           91.75                0         0.2        0.1
#>      PrEarlyFut PrEfficacy PrFutility PrGrayZone
#> [1,]        0.1        0.1       0.15       0.75

## 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 come from a Beta(10, 10) distribution
  ## and the historical data will be generated randomly from the resulting
  ## binomial distribution!

  ## delta will be 0.15

  ## 200 simulations

  ## the different prior sample sizes (historical data size) e.g. 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 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,
    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, ],
        oc3(
          method = method,
          histSize = histSize,
          trialSize = trialSize,
          controlRateDist = function() {
            rbeta(n = 1, 10, 10)
          },
          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,
                        # name of the output file
                        out = "output.docx",
                        # 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")
    writeDoc(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,
    (histSize == 150) &
      (method %in% ppNames)
  )
}