inst/doc/Simulation_Comparison.R

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

## ----collapse=TRUE------------------------------------------------------------
library(BayesianMCPMod)
library(DoseFinding)
library(MCPModPack)

library(tidyr)
library(dplyr)

library(kableExtra)
library(ggplot2)

library(doFuture)

## ----echo = FALSE-------------------------------------------------------------

load_path <- file.path(getwd(), "data", "simulation_comparison")


## ----collapse=TRUE------------------------------------------------------------
doses_sim     <- c(0, 1, 2, 4, 8) 
n_sample      <- c(40, 40, 40, 40, 40)
sd_sim        <- 0.4

max_dose      <- max(doses_sim)    
plc_eff_guess <- 0    

alpha         <- 0.05

## ----collapse=TRUE------------------------------------------------------------
set.seed(7015)
n_sim <- 10000

plan(multisession)
registerDoFuture()

## ----echo = FALSE-------------------------------------------------------------

n_sim <- 10


## ----fig.height=5, collapse=TRUE----------------------------------------------
emax_guess     <- guesst(d = doses_sim[2], p = 0.6, "emax") 
exp_guess      <- guesst(d = doses_sim[2], p = 0.05, model = "exponential", Maxd = max_dose)
logit_guess    <- guesst(d = c(doses_sim[2], doses_sim[3]), p = c(0.1, 0.9), "logistic", Maxd = max_dose) 
sig_emax_guess <- guesst(d = c(doses_sim[2], doses_sim[3]), p = c(0.15, 0.75), model = "sigEmax") 

plot(Mods(linear      = NULL,
          exponential = exp_guess,
          emax        = emax_guess,
          logistic    = logit_guess,
          sigEmax     = sig_emax_guess,
          doses       = doses_sim,
          placEff     = plc_eff_guess,
          direction   = "increasing"),
     main = "Candidate Models")

## ----collapse=TRUE------------------------------------------------------------
exp_eff <- c(0.0001, 0.05, 0.1, 0.2, 0.3, 0.5)

## ----collapse=TRUE------------------------------------------------------------
# Simulation parameters
sim_parameters <- list(n            = n_sample,
                       doses        = doses_sim,
                       dropout_rate = 0.0,
                       go_threshold = 0.1,
                       nsims        = n_sim)

# Candidate dose - response models
models_MCPModPack = list(linear      = NA,
                         exponential = exp_guess,
                         emax        = emax_guess,
                         logistic    = logit_guess,
                         sigemax     = sig_emax_guess)

# Assumed dose - response models (models will be added in loop)
sim_models_part <- list(max_effect     = exp_eff,
                        sd             = rep(sd_sim, length(doses_sim)),
                        placebo_effect = plc_eff_guess)

# Parallelization across assumed dose - response models
powers_MCPModPack_eff <- foreach(
  k = seq_along(models_MCPModPack),
  .combine = cbind, 
  .options.future = list(seed = TRUE)) %dofuture% {
    
    sim_model_k <- c(models_MCPModPack[k], sim_models_part)
    
    MCPModSimulation(endpoint_type   = "Normal",
                     models          = models_MCPModPack,
                     alpha           = alpha,
                     direction       = "increasing",
                     model_selection = "aveAIC",
                     Delta           = 0.1,
                     sim_models      = sim_model_k,
                     sim_parameters  = sim_parameters)$sim_results$power
  
  }

# Post-processing result for printing
colnames(powers_MCPModPack_eff) <- names(models_MCPModPack)
results_MCPModPack_eff <- cbind(max_eff = round(exp_eff, digits = 2),
                                powers_MCPModPack_eff) %>%
  data.frame() %>%
  rename(sigEmax = sigemax) %>%
  mutate(average = rowMeans(select(., linear:sigEmax)))

## ----echo = FALSE-------------------------------------------------------------

results_MCPModPack_eff <- readRDS(file.path(load_path, "results_MCPModPack_eff.rds"))


## ----collapse=TRUE------------------------------------------------------------
# Vague prior specification
prior_list_vague <- rep(list(RBesT::mixnorm(comp1 = c(w = 1, m = 0, n = 1),
                                            sigma = sd_sim, param = "mn")),
                        times = length(doses_sim))
names(prior_list_vague) <- c("Ctrl", "DG_1", "DG_2", "DG_3", "DG_4")

# Parallelization across the expected effects for maximum dose
success_rates_BayesianMCPMod_eff <- foreach(
  k               = seq_along(exp_eff),
  .combine        = rbind, 
  .options.future = list(seed = TRUE)) %dofuture% {
    
    exp_eff_k <- exp_eff[k]
    
    models_BayesianMCPMod <- Mods(linear      = NULL,
                                  exponential = exp_guess,
                                  emax        = emax_guess,
                                  logistic    = logit_guess,
                                  sigEmax     = sig_emax_guess,
                                  doses       = doses_sim,
                                  placEff     = plc_eff_guess,
                                  maxEff      = exp_eff_k,
                                  direction   = "increasing")
       
    # Optimal contrasts
    contr <- getContr(mods         = models_BayesianMCPMod,
                      dose_levels  = doses_sim,
                      prior_list   = prior_list_vague,
                      dose_weights = rep(1, length(doses_sim)))
  
    # Perform Simulations
    sim_result <- assessDesign(n_patients     = n_sample,
                               mods           = models_BayesianMCPMod,
                               prior_list     = prior_list_vague,
                               sd             = sd_sim,
                               n_sim          = n_sim,
                               alpha_crit_val = alpha,
                               contr          = contr)
    
    c(sapply(sim_result, attr, "successRate"),
      average = attr(sim_result, "avgSuccessRate"))
  
  }

# Post-processing result for printing
rownames(success_rates_BayesianMCPMod_eff) <- NULL
results_BayesianMCPMod_eff <- data.frame(cbind(max_eff = round(exp_eff, digits = 2),
                                           success_rates_BayesianMCPMod_eff))

## ----echo = FALSE-------------------------------------------------------------

results_BayesianMCPMod_eff <- readRDS(file.path(load_path, "results_BayesianMCPMod_eff.rds"))


## ----fig.height=5-------------------------------------------------------------
## pre-processing the data
df_plot_eff <- rbind(results_MCPModPack_eff %>% 
                     mutate(package_name = "MCPModPack"),
                     results_BayesianMCPMod_eff %>% 
                     mutate(package_name = "BayesianMCPMod")) %>%
  pivot_longer(cols      = names(results_BayesianMCPMod_eff)[-1],
               names_to  = "model_shape",
               values_to = "success_rate") %>%
  filter(model_shape != "average") %>%
  spread(key = package_name, value = success_rate)  %>%
  mutate(model_shape = as.factor(model_shape),
         max_eff     = as.factor(max_eff)) %>%
  group_by(max_eff) %>%
  mutate(offset = 0.1 * (seq_along(model_shape) - ceiling(length(model_shape) / 2)))

# Plot with short horizontal dashes for each model_shape
ggplot(df_plot_eff, aes(x = as.numeric(max_eff) + offset, y = BayesianMCPMod, color = model_shape)) +
  geom_segment(aes(x = as.numeric(max_eff) + offset - 0.05, xend = as.numeric(max_eff) + offset + 0.05,
                   y = BayesianMCPMod, yend = BayesianMCPMod)) +
  geom_segment(aes(xend = as.numeric(max_eff) + offset, yend = MCPModPack)) +
  scale_x_continuous(breaks = unique(as.numeric(df_plot_eff$max_eff)), labels = levels(df_plot_eff$max_eff)) +
  labs(title    = "Comparing Power and Success Rate",
       subtitle = "For Different Expected Effects for Maximum Dose",
       x        = "Expected Effect for Maximum Dose",
       y        = "Success Rate / Power",
       color    = "Assumed True Model Shapes",
       linetype = "Type",
       caption  = "Horizontal lines represent the BayesianMCPMod success rates and vertical lines mark the distance to the MCPModPack power.") +
  theme_minimal() + 
  theme(legend.position = "bottom")

## -----------------------------------------------------------------------------
kable(results_MCPModPack_eff) %>%
  kable_classic() %>%
    add_header_above(c("Power Values Across Different Expected Effects" = 7),
                     font_size = 15, bold = TRUE) %>%
    add_header_above(c("MCPModPack" = 7), font_size = 15, bold = TRUE)

## -----------------------------------------------------------------------------
kable(results_BayesianMCPMod_eff) %>%
  kable_classic(full_width = TRUE) %>%
    add_header_above(c("Success Rates Across Different Expected Effects" = 7),
                     font_size = 15, bold = TRUE) %>%
    add_header_above(c("BayesianMCPMod" = 7), font_size = 15, bold = TRUE) 

## -----------------------------------------------------------------------------
exp_eff_fix <- 0.2

## -----------------------------------------------------------------------------
n_sim_vec <- seq(2, 10, 1)^4

## ----echo = FALSE-------------------------------------------------------------
n_sim_vec <- c(5, 10)

## -----------------------------------------------------------------------------
# Updating the assumed dose - response models for the fixed expected effect (models will be added in loop)
sim_models_part$max_effect <- exp_eff_fix

# Parallelization across assumed dose - response models & number of simulations
## Reversing the numbers speeds up parallelization in case there are not enough workers
sim_grid <- expand.grid(names(models_MCPModPack), rev(n_sim_vec))
powers_MCPModPack_conv <- foreach(
  k = seq_len(nrow(sim_grid)),
  .combine = c,
  .options.future = list(seed = TRUE)) %dofuture% {
    
    sim_model_k <- c(models_MCPModPack[sim_grid[k, 1]], sim_models_part)
    
    sim_params_k <- sim_parameters
    sim_params_k$nsims <- sim_grid[k, 2]
    
    MCPModSimulation(endpoint_type   = "Normal",
                     models          = models_MCPModPack,
                     alpha           = alpha,
                     direction       = "increasing",
                     model_selection = "aveAIC",
                     Delta           = 0.1,
                     sim_models      = sim_model_k,
                     sim_parameters  = sim_params_k)$sim_results$power
  
  }

#Post-processing
results_MCPModPack_conv <- cbind(sim_grid, powers_MCPModPack_conv) %>%
  arrange(desc(row_number())) %>%
  rename(model_name   = Var1,
         n_sim        = Var2,
         success_rate = powers_MCPModPack_conv) %>%
  mutate(model_name = if_else(model_name == "sigemax",
                             "sigEmax", model_name)) %>%
  mutate(package_name = "MCPModPack")

## ----echo = FALSE-------------------------------------------------------------

results_MCPModPack_conv <- readRDS(file.path(load_path, "results_MCPModPack_conv.rds"))


## -----------------------------------------------------------------------------
# Model specifications with fixed expected effect
models_BayesianMCPMod <- Mods(linear      = NULL,
                              exponential = exp_guess,
                              emax        = emax_guess,
                              logistic    = logit_guess,
                              sigEmax     = sig_emax_guess,
                              doses       = doses_sim,
                              placEff     = plc_eff_guess,
                              maxEff      = exp_eff_fix,
                              direction   = "increasing")
       
# Optimal contrasts
contr <- getContr(mods         = models_BayesianMCPMod,
                  dose_levels  = doses_sim,
                  prior_list   = prior_list_vague,
                  dose_weights = rep(1, length(doses_sim)))

# Perform Simulations
sim_result <- assessDesign(n_patients     = n_sample,
                           mods           = models_BayesianMCPMod,
                           prior_list     = prior_list_vague,
                           sd             = sd_sim,
                           n_sim          = max(n_sim_vec),
                           alpha_crit_val = alpha,
                           contr          = contr)

# Getting success rates at different numbers of simulations
results_BayesianMCPMod_conv <- sapply(sim_result, function (model_result) {
  
  sapply(n_sim_vec, function (n_sim_x) {
    
    success_rate_x <- mean(model_result[seq_len(n_sim_x), 1])
    
  })
  
}) %>%
  # Post-processing
  as.data.frame %>%
  mutate(n_sim = n_sim_vec) %>%
  pivot_longer(cols = -n_sim,
               names_to = "model_name",
               values_to = "success_rate") %>%
  mutate(package_name = "BayesianMCPMod")

## ----echo = FALSE-------------------------------------------------------------

results_BayesianMCPMod_conv <- readRDS(file.path(load_path, "results_BayesianMCPMod_conv.rds"))


## ----fig.height=5-------------------------------------------------------------
df_plot_conv <- inner_join(results_BayesianMCPMod_conv,
                           results_MCPModPack_conv,
                           by = c("model_name", "n_sim")) %>%
  mutate(success_rate_diff = success_rate.x - success_rate.y) %>%
  select(model_name, n_sim, success_rate_diff)

ggplot(df_plot_conv, aes(x = n_sim, y = success_rate_diff, color = model_name)) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey")+
  scale_x_continuous(breaks = unique(df_plot_conv$n_sim)[-c(2, 3)]) +
  scale_y_continuous(breaks = c(-0.05, 0, 0.05, 0.1, 0.15)) +
  labs(
    title   = "Success Rate Difference vs Number of Simulations",
    x       = "Number of Simulations",
    y       = "Success Rate Difference",
    color   = "True Model Shape",
    caption = "The success rate difference is the difference of the success rates calculated 
    with BayesianMCPMod and the power values calculated with MCPModPack.") +
  theme_minimal() + 
  theme(panel.grid.minor.x = element_blank())