inst/doc/Simulation_Example_Binary.R

In BayesianMCPMod: Simulate, Evaluate, and Analyze Dose Finding Trials with Bayesian MCPMod

## ----setup--------------------------------------------------------------------
suppressPackageStartupMessages({
  library(BayesianMCPMod)
  library(RBesT)
  library(DoseFinding)
  library(clinDR)
  library(dplyr)
  library(tibble)
  library(reactable)
})

set.seed(7015)

## ----eval = FALSE-------------------------------------------------------------
# future::plan(future::multisession, workers = 4L)

## -----------------------------------------------------------------------------
data("metaData")
testdata    <- as.data.frame(metaData)
dataset     <- filter(testdata, bname == "XELJANZ")
histcontrol <- filter(dataset, dose == 0, primtime == 12, indication == "RHEUMATOID ARTHRITIS")

hist_data <- data.frame(
  trial = histcontrol$nctno,
  est   = histcontrol$rslt,
  se    = histcontrol$se,
  sd    = histcontrol$sd,
  r     = round(histcontrol$sampsize*histcontrol$rslt),
  n     = histcontrol$sampsize)

sd_tot <- with(hist_data, sum(sd * n) / sum(n))

## -----------------------------------------------------------------------------
dose_levels <- c(0, 2.5, 5, 10,20)

# 1) Establish MAP prior (beta mixture distribution) 
set.seed(7015) # re-sets seed only for this example; remove in your analysis script
map <- gMAP(
  cbind(hist_data$r, hist_data$n - hist_data$r) ~ 1 | histcontrol$nctno,
  family     = binomial,
  tau.dist   = "HalfNormal",
  tau.prior  = 0.5,
  beta.prior = (1 / sqrt(0.1 * 0.9)),
  warmup     = 1000,
  iter       = 10000,
  chains     = 2,
  thin       = 1
)

prior <- automixfit(map) #fits mixture distribution from MCMC samples from above

p <- summary(prior)[1]

# 2) Robustify prior
prior_rob <- RBesT::robustify(
      priormix = prior,
      mean     = 0.5,
      weight   = 0.2)

# 3) Translate prior to logit scale (to approximate via normal mixture model)
r         <- rmix(prior_rob, n = 1e4)
log_r     <- logit(r)
prior_ctr <- automixfit(log_r, type = "norm")

# Specification of reference scale (this follows the idea of [@Neuenschwander2016]). 
sigma(prior_ctr) <- sqrt(1 / (p * (1 - p)))

# Specify a prior list
prior_trt <- RBesT::mixnorm(
    comp1 = c(w = 1, m = logit(summary(prior)[1]), n = 1),
    sigma = sqrt(1/(p*(1-p))),
    param = "mn")
  
  prior_list <- c(list(prior_ctr),
                  rep(x     = list(prior_trt),
                      times = length(dose_levels[-1])))

dose_names        <- c("Ctr", paste0("DG_", seq_along(dose_levels[-1])))
names(prior_list) <- dose_names

## ----Setting Vague Prior Without Execution, eval = TRUE-----------------------
prior_list_vague <- rep(list(RBesT::mixnorm(
  comp1 = c(w = 1, m = logit(p), n = 1),
  sigma = sqrt(1 / (p * (1 - p))),
  param = "mn"
)), times = length(dose_levels))
names(prior_list_vague) <- c("Ctrl", "DG_1", "DG_2", "DG_3", "DG_4")

## -----------------------------------------------------------------------------
n_patients <- c(30, 40, 40, 40, 40)

models <- Mods(
  linear      = NULL,
  sigEmax     = c(10, 5),
  logistic    = c(11, 15),
  exponential = 10,
  emax        = 2,
  doses       = dose_levels,
  placEff     = RBesT::logit(0.18),
  maxEff      = (RBesT::logit(0.48) - RBesT::logit(0.18))
)

## -----------------------------------------------------------------------------
set.seed(7015) # re-sets seed only for this example; remove in your analysis script
success_probabilities <- assessDesign(
  n_patients        = n_patients,
  mods              = models,
  prior_list        = prior_list,
  probability_scale = TRUE,
  delta             = 0.2,
  n_sim             = 100) # speed up example run-time
success_probabilities

## -----------------------------------------------------------------------------
set.seed(7015) # re-sets seed only for this example; remove in your analysis script
success_probabilities_uneq <- assessDesign(
  n_patients        = c(40, 30, 30, 30, 50),
  mods              = models,
  prior_list        = prior_list,
  probability_scale = TRUE,
  delta             = 0.2,
  n_sim             = 100) # speed up example run-time
success_probabilities_uneq

## -----------------------------------------------------------------------------
set.seed(7015) # re-sets seed only for this example; remove in your analysis script
success_probabilities_vague <- assessDesign(
  n_patients        = n_patients,
  mods              = models,
  prior_list        = prior_list_vague,
  probability_scale = TRUE,
  delta             = 0.2,
  n_sim             = 100) # speed up example run-time
success_probabilities_vague