Skip to contents

Simulation Example of Bayesian MCPMod and MCPMod

Source: vignettes/Comparison_vignette.qmd
setup
#' Display Parameters Table
#'
#' This function generates a markdown table displaying the names and values of parameters
#' from a named list.
#'
#' @param named_list A named list where each name represents a parameter name and the list
#'   element represents the parameter value. Date values in the list are automatically
#'   converted to character strings for display purposes.
#'
#' @return Prints a markdown table with two columns: "Parameter Name" and "Parameter Values".
#'   The function does not return a value but displays the table directly to the output.
#'
#' @importFrom knitr kable
#' @examples
#' params <- list("Start Date" = as.Date("2020-01-01"),
#'                "End Date" = as.Date("2020-12-31"),
#'                "Threshold" = 10)
#' display_params_table(params)
#'
#' @export
display_params_table <- function(named_list) {
  display_table <- data.frame()
  value_names <- data.frame()
  for (i in 1:length(named_list)) {
    # dates will display as numeric by default, so convert to char first
    if (class(named_list[[i]]) == "Date") {
      named_list[[i]] = as.character(named_list[[i]])
    }
    if (!is.null(names(named_list[[i]]))) {
      value_names <- rbind(value_names, paste(names(named_list[[i]]), collapse = ', '))
    }
    values <- data.frame(I(list(named_list[[i]])))
    display_table <- rbind(display_table, values)
  }
  
  round_numeric <- function(x, digits = 3) {
    if (is.numeric(x)) {
      return(round(x, digits))
    } else {
      return(x)
    }
  }
  
  display_table[1] <- lapply(display_table[1], function(sublist) {
    lapply(sublist, round_numeric)
  })
  
  class(display_table[[1]]) <- "list"
  
  if (nrow(value_names) == 0) {
    knitr::kable(
      cbind(names(named_list), display_table),
      col.names = c("Name", "Value")
    )
  } else {
    knitr::kable(
      cbind(names(named_list), value_names, display_table),
      col.names = c("Name", "Value Labels", "Value")
    )
  }
}
# function to solve power results for tables (different max eff) BayesianMCPMod
# return(successrates models, average)
extract_success_rates <- function(results_list, models) {
  success_rates <- list()

  for (i in seq_along(results_list)) {
    success_rate <- c()
    for (model in models) {
      success_rate <- c(success_rate, attr(results_list[[i]][[model]]$BayesianMCP,"successRate"))
    }
    success_rates[[paste0("Bay_",  attr(results_list[[i]],"maxEff"))]] <- c(success_rate, attr(results_list[[i]], "avgSuccessRate"))
  }

  return(success_rates)
  }

# function to solve power results for tables (different nsample) BayesianMCPMod
#models <c-(min_scenario)
extract_success_rates_nsample <- function(results_list, models) {
  success_rates <- list()

  for (i in seq_along(results_list)) {
    success_rate <- c()
    for (model in models) {
      success_rate <- c(success_rate, attr(results_list[[i]][[model]]$BayesianMCP,"successRate"))
    }
    success_rates[[paste0("Bay_", i )]] <- c(success_rate,  attr(results_list[[i]],"avgSuccessRate"))
  }

  return(success_rates)
}

#input: result_list, models =
extract_success_rates_nsim <- function(results_list, models, n_sim) {
  success_rates <- list()
  for (model in models) {
    success_rate <- c()
    for (i in seq_along(n_sim)) {

      success_rate <- c(success_rate, attr(results_list[[i]][[model]]$BayesianMCP,"successRate"))
    }
    success_rates[[paste0("Bay_", model)]] <- c(success_rate)
  }

  return(success_rates)
}


print_result_Bay_max_eff <- function(results, scenario, variable) {

  result_table <- t(data.table(
    Bay_0.05 = results$Bay_0.05,
    Bay_0.1 = results$Bay_0.1,
    Bay_0.2 = results$Bay_0.2,
    Bay_0.3 = results$Bay_0.3,
    Bay_0.5 = results$Bay_0.5))

  result_table <- as.data.table(result_table)
  names(result_table) <- scenario
  #return(result_table)

  kable_result <- kable(cbind(variable, result_table))%>%
    kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different expected effects " = length(scenario)+1), font_size = 15, bold = TRUE)%>%
    add_header_above(c("BayesianMCPMod " = length(scenario)+1), font_size = 15, bold = TRUE)

  list(result_table = result_table, kable_result = kable_result)
}



print_result_Bay_nsample <- function(results, scenario, variable) {

  result_table <- t(data.table(
    Bay_1 = results$Bay_1,
    Bay_2 = results$Bay_2,
    Bay_3 = results$Bay_3,
    Bay_4 = results$Bay_4))

  result_table <- as.data.table(result_table)
  names(result_table) <- scenario
  #return(result_table)

  kable_result <- kable(cbind(variable, result_table))%>%
    kable_classic(full_width = TRUE)%>%
    add_header_above(c("Success probability results different sample sizes " = length(scenario)+1), font_size = 15, bold = TRUE)%>%
    add_header_above(c("BayesianMCPMod " = length(scenario)+1), font_size = 15, bold = TRUE)

  list(result_table = result_table, kable_result = kable_result)
}

plot_power_deviation <- function(data, x, xlab){
  plot <- ggplot2::ggplot(data, aes(x = x, y = value, color = variable, group = variable))+
    geom_point()+
    geom_line() +
    scale_x_continuous(breaks = c(0, 100, 500, 1000, 2500, 5000, 10000),
    labels = c("", "100", "", "1000", "2500", "5000", "10000")) +
    geom_hline(aes(yintercept = 0), linetype = 2)+
    geom_hline(aes(yintercept = -0.05), linetype = 2, color = "darkgrey")+
    geom_hline(aes(yintercept = 0.05), linetype = 2, color = "darkgrey")+
    geom_hline(aes(yintercept = -0.1), linetype = 2, color = "darkgrey")+
    geom_hline(aes(yintercept = 0.1), linetype = 2, color = "darkgrey")+
    scale_color_manual(name = "assumed true model", values = c("linear" = "red", "exponential" = "blue", "emax" = "darkgreen", "logistic" = "orange", "sigemax" = "purple", "beta" = "deepskyblue", "quadratic" = "deeppink"))+
    labs(x = xlab, y = "power deviation")+
    ylim(-0.15, 0.15)+        
    theme_classic()
  return(plot)
}




# parallel----------------------------------------------
chunkVector <- function (x, n_chunks) {
  if (n_chunks <= 1) {
    chunk_list <- list(x)
  } else {
    chunk_list <- unname(split(x, cut(seq_along(x), n_chunks, labels = FALSE)))
  }
  return(chunk_list)
}

library(BayesianMCPMod)
library(RBesT)
library(clinDR)
library(dplyr)
library(tibble)
library(reactable)
library(DoseFinding)
library(MCPModPack)
library(kableExtra)
library(data.table)
library(doFuture)
library(doRNG)

registerDoFuture()
plan(multisession)



set.seed(7015)

1 Introduction

This vignette demonstrates the application of the {BayesianMCPMod} package for sample size calculations and the comparison with the {MCPModPack} package. As for other bayesian approaches BMCPMod is set up in a way that it mimics the results (and operating characteristics) of the frequentist MCPMod for non-informative priors. This characteristic is illustrated in the following sections focusing on the trial planning.

In order to compare BayesianMCPMod and MCPModPack success probabilities four different dose-finding scenarios are considered. For the first three scenarios the general trial design is the same, the only difference is the considered set of DR models. These design choices are:

  • four dose levels plus placebo (0 mg, 1 mg, 2 mg, 4 mg, 8 mg)

  • total sample size of N = 200

  • equal allocation ratio for each dose group

  • standard deviation of 0.4 for every dose group

The fourth scenario, the variability scenario, has some different design choices for the number of dose levels and the total sample size: - three dose levels plus placebo (0 mg, 1 mg, 4 mg, 8 mg) - total sample size of N = 100

Building on these scenarios the following varying additional parameters are used:

  • expected effect for maximum dose of (0.05, 0.1, 0.2, 0.3, 0.5)

  • different total sample size with different allocation:

  • (40,20,20,20,40), (30,30,30,30,30), (48,24,24,24,48), (36,36,36,36,36) (for the first three scenarios)

  • (40,20,20,40), (30,30,30,30), (48,24,24,48), (36,36,36,36) (for the variability scenario)

In a subsequent step, we extend our analysis to test the convergence of power values as the number of simulations increases. For this part of the study, we assume an expected effect for the maximum dose of 0.2. The sample size is set to 200 for the linear, monotonic and non-monotonic scenario and 100 for the variability scenario using equal allocation.

For each of these simulation cases, 10000 simulation runs are performed. Given the number of simulations, and using basic mathematical approximations (like the law of large numbers) the difference in success probabilities should be in the range of 1-3%.

setup
#######  General assumptions  #########
# simulations
n_sim <- 100 #to be upscaled to 10000

# Define the list of sample sizes
nsample_list = list(20*c(2,1,1,1,2),30*c(1,1,1,1,1), 24*c(2,1,1,1,2),36*c(1,1,1,1,1))
nsample_vector <- c("(40,20,20,20,40)", "(30,30,30,30,30)", "(48,24,24,24,48)",  "(36,36,36,36,36)")

# define input parameters
doses.sim <- c(0,1,2,4,8)   # dose levels
max.dose <- max(doses.sim) # specify max dose possibly available
plc.guess <-  0 # expected placebo effect
Nsample <-c(40,40,40,40,40) #,36*c(1,1,1,1)) #list(20*c(2,1,1,2),30*c(1,1,1,1), 24*c(2,1,1,2),36*c(1,1,1,1))
expectedEffect_fix <- 0.2  #c(0.05,0.1,0.2,0.3,0.5)
expectedEffect <- c(0.05,0.1,0.2,0.3,0.5)
sd.sim <- c(0.4)

#input paramters variability scenario
doses_var <- c(0, 1, 4, 8)
sd.sim_var <- 0.4
Nsample_var <- c(25, 25, 25, 25)
nsample_list_var <- list(20*c(2,1,1,2),30*c(1,1,1,1), 24*c(2,1,1,2),36*c(1,1,1,1))


# define model functions  
emax.g <- guesst(d = doses.sim[2], p = 0.6, "emax") 
exp.g <- guesst(d=doses.sim[2],p=0.05,model="exponential", Maxd=max.dose)
logit.g <- guesst(d=c(doses.sim[2],doses.sim[3]),p=c(0.1,0.9),"logistic",Maxd=max.dose) 
sigEmax.g <- guesst(d=c(doses.sim[2],doses.sim[3]),p=c(0.15,0.75), model = "sigEmax") 
quad.g <- guesst(d=doses.sim[2],p=0.35,model="quadratic")
beta.g <- guesst(d=doses.sim[2],p=0.2,model="betaMod", dMax=5.5, scal=9.6, Maxd=max.dose)

### define models to be tested in MCPMod
models <- Mods(linear=NULL, 
               emax = emax.g,
               exponential = exp.g,
               logistic=logit.g,
               sigEmax = sigEmax.g,
               quadratic = quad.g,
               betaMod = beta.g,
               doses = doses.sim,
               placEff = plc.guess, 
               maxEff = expectedEffect_fix,
               direction = "increasing",
               addArgs = list(scal=9.6) 
)

alpha <- 0.05

# display_params_table(models)


# kable(t(as.data.table(cbind(doses.sim, Nsample))))%>%
#     kable_classic(full_width = TRUE)%>%
#     add_header_above(c( "Doses" = 6), font_size = 15, bold = TRUE)
# 
# 
# kable(t(data.table(placebo_guess = plc.guess,
#                  sd = sd.sim[1],
#                  alpha = alpha )))%>%
#     kable_classic(full_width = TRUE)

1.1 Scenarios

In the following the candidate models for the different scenarios are specified and plotted.

setup
#Specification of considered models for different scenarios for BayesianMCPMod

min_scenario <- c("linear", "exponential", "emax")
min_models <- Mods(linear=NULL, 
                   exponential = exp.g,
                  emax = emax.g,
                  doses = doses.sim,
                  placEff = plc.guess, 
                  maxEff = expectedEffect_fix,
                  direction = "increasing"
)
plot(min_models, main = "minimal Scenario")

setup
monotonic_scenario <- c("linear", "exponential", "emax", "logistic", "sigEmax")

monotonic_models <- Mods(linear=NULL, 
               exponential = exp.g,
               emax = emax.g,
               logistic=logit.g,
               sigEmax = sigEmax.g,
               doses = doses.sim,
               placEff = plc.guess, 
               maxEff = expectedEffect_fix,
               direction = "increasing"
)
plot(monotonic_models, main = "monotonic Scenario")

setup
non_monotonic_scenario <- c("linear", "emax", "sigEmax", "quadratic", "betaMod")

non_monotonic_models <- Mods(linear=NULL, 
               emax = emax.g,
               sigEmax = sigEmax.g,
               quadratic = quad.g,
               betaMod = beta.g,
               doses = doses.sim,
               placEff = plc.guess, 
               maxEff = expectedEffect_fix,
               direction = "increasing",
               addArgs = list(scal=9.6)
)
plot(non_monotonic_models,main = "non monotonic Scenario")

setup
variability_scenario <- c("linear","exponential", "emax")

kable(t(as.data.table(cbind(doses_var, Nsample_var))))%>%
    kable_classic(full_width = TRUE)%>%
    add_header_above(c( "Doses variability scenario" = 5), font_size = 15)
Doses variability scenario
doses_var 0 1 4 8
Nsample_var 25 25 25 25
setup
var_models <- Mods(linear=NULL, 
               exponential = exp.g,
               emax = emax.g,
               doses = doses_var,
               placEff = plc.guess, 
               maxEff = expectedEffect_fix,
               direction = "increasing"
)
plot(var_models, main = "variability Scenario")

2 MCPModPack

setup
#Specification of considered models for different scenarios for MCPModPack

min_modelsPack = list(linear = NA,
              exponential = 4.447149,
              emax = 0.6666667)

monotonic_modelsPack = list(linear = NA,
              exponential = 4.447149,
              emax = 0.6666667,
              logistic = c(1.5, 0.2275598),
              sigEmax =  c(1.528629, 4.087463))

non_monotonic_modelsPack = list(linear = NA,
              emax = 0.6666667 ,
              sigEmax =  c(1.528629, 4.087463),
              quadratic = -0.09688711)
#beta

var_modelsPack = list(linear = NA,
              exponential = 4.447149,
              emax = 0.6666667,
              logistic = c(1.5, 0.2275598),
              sigEmax =  c(1.528629, 4.087463))

The following simulations will be conducted utilizing the MCPModPack package (with varying expected effect for maximum dose and sample sizes).

2.1 Minimal scenario

2.1.1 varying expected effect for maximum dose

setup
# assumed dose - response model
sim_models_min_linear = list(linear = NA,
                             exponential = 4.447149,
                             max_effect = c(0.05,0.1,0.2,0.3,0.5),
                             sd = rep(sd.sim, length(doses.sim)),
                             placebo_effect = plc.guess)

sim_models_min_exp = list(exponential = 4.447149,
                  max_effect = c(0.05,0.1,0.2,0.3,0.5),
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)

sim_models_min_emax = list( emax = 0.6666667,
                  max_effect = c(0.05,0.1,0.2,0.3,0.5),
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)

# Simulation parameters
sim_parameters = list(n = Nsample, 
                      doses = doses.sim,
                      dropout_rate = 0.0, 
                      go_threshold = 0.1, 
                      nsims = n_sim)

# Initialize list to store results
results_list_min <- list()
  
# Run simulation
func_sim <- function(models ,sim_models, sim_parameters){

sim_result = MCPModSimulation(endpoint_type = "Normal", 
                                 models = models, 
                                 alpha = alpha, 
                                 direction = "increasing", 
                                 model_selection = "aveAIC", 
                                 Delta = 0.1,
                                 sim_models = sim_models,
                                 sim_parameters = sim_parameters)
}


list_models <- list(sim_models_min_linear, sim_models_min_exp, sim_models_min_emax)
chunks <- chunkVector(seq_along(list_models), getDoParWorkers())

results_list_min <- foreach(k = chunks, .combine = c)  %dorng% {
  
  lapply(k, function (i) {
    
    func_sim(min_modelsPack, list_models[[i]], sim_parameters) 
    
  })
  
}




#store results
results_min_MCP <- data.table(
  max_eff = c(0.05, 0.1, 0.2, 0.3, 0.5),
  linear = c(results_list_min[[1]]$sim_results$power),
  exponential = c(results_list_min[[2]]$sim_results$power),
  emax = c(results_list_min[[3]]$sim_results$power),
  Average = c(sum(results_list_min[[1]]$sim_results$power[1], results_list_min[[2]]$sim_results$power[1], results_list_min[[3]]$sim_results$power[1])/3,
               sum(results_list_min[[1]]$sim_results$power[2], results_list_min[[2]]$sim_results$power[2], results_list_min[[3]]$sim_results$power[2])/3,
               sum(results_list_min[[1]]$sim_results$power[3], results_list_min[[2]]$sim_results$power[3], results_list_min[[3]]$sim_results$power[3])/3,
               sum(results_list_min[[1]]$sim_results$power[4], results_list_min[[2]]$sim_results$power[4], results_list_min[[3]]$sim_results$power[4])/3,
               sum(results_list_min[[1]]$sim_results$power[5], results_list_min[[2]]$sim_results$power[5], results_list_min[[3]]$sim_results$power[5])/3))


kable(results_min_MCP)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different expected effects " = 5), font_size = 15, bold = TRUE)%>%
    add_header_above(c("Minimal scenario " = 5), font_size = 15, bold = TRUE)
Minimal scenario
Power results different expected effects
max_eff linear exponential emax Average
0.05 0.13 0.16 0.17 0.1533333
0.10 0.29 0.43 0.31 0.3433333
0.20 0.80 0.77 0.82 0.7966667
0.30 0.99 0.97 0.98 0.9800000
0.50 1.00 1.00 1.00 1.0000000

2.1.2 varying sample size

setup
#list of different samole sizes
nsample_list = list(20*c(2,1,1,1,2),30*c(1,1,1,1,1), 24*c(2,1,1,1,2),36*c(1,1,1,1,1))

# assumed dose - response model
sim_models_min_linear = list(linear = NA,
                  max_effect = expectedEffect_fix,
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)

sim_models_min_exp = list(exponential = 4.447149,
                  max_effect = expectedEffect_fix,
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)

sim_models_min_emax = list( emax = 0.6666667,
                  max_effect = expectedEffect_fix,
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)

# Simulation parameters
sim_parameters = list(n = Nsample, 
                      doses = doses.sim,
                      dropout_rate = 0.0, 
                      go_threshold = 0.1, 
                      nsims = n_sim)


# store results
power_min_nsample_linear <- vector()
power_min_nsample_exp <- vector()
power_min_nsample_emax <- vector()





# # Loop over nsample values
# for (i in seq_along(nsample_list)) {
#   
#   # Update nsample in simulation parameters
#   sim_parameters$n <- nsample_list[[i]]
#   
#   # Set max_effect to 0.2
#   min_modelsPack$max_effect <- expectedEffect_fix
#   
#   # power values for different assumed true models
#   power_nsample_linear <- func_sim(min_modelsPack, sim_models_min_linear, sim_parameters)
#   power_min_nsample_linear[i] <- power_nsample_linear$sim_results$power
#   
#   power_nsample_exp <- func_sim(min_modelsPack, sim_models_min_exp, sim_parameters)
#   power_min_nsample_exp[i] <-power_nsample_exp$sim_results$power
#   
#   power_nsample_emax <- func_sim(min_modelsPack, sim_models_min_emax, sim_parameters)
#   power_min_nsample_emax[i] <- power_nsample_exp$sim_results$power
#   
# }
# Set max_effect to 0.2
min_modelsPack$max_effect <- expectedEffect_fix
#list of different samole sizes
nsample_list = list(20*c(2,1,1,1,2),30*c(1,1,1,1,1), 24*c(2,1,1,1,2),36*c(1,1,1,1,1),
                    20*c(2,1,1,1,2),30*c(1,1,1,1,1), 24*c(2,1,1,1,2),36*c(1,1,1,1,1),
                    20*c(2,1,1,1,2),30*c(1,1,1,1,1), 24*c(2,1,1,1,2),36*c(1,1,1,1,1))

list_models <- list(sim_models_min_linear,sim_models_min_linear,sim_models_min_linear,sim_models_min_linear,
                    sim_models_min_exp, sim_models_min_exp,sim_models_min_exp,sim_models_min_exp,
                    sim_models_min_emax, sim_models_min_emax, sim_models_min_emax, sim_models_min_emax)

#ind_tasks <- 1:12

chunks <- chunkVector(seq_along(list_models), getDoParWorkers())

results_list_min_nsample <- foreach(k = chunks, .combine = c, .export = c(as.character(min_modelsPack), as.character(nsample_list)))  %dorng% {
  lapply(k, function (i) {
    sim_parameters = list(n = nsample_list[[i]], 
                      doses = doses.sim,
                      dropout_rate = 0.0, 
                      go_threshold = 0.1, 
                      nsims = n_sim)
    
    func_sim(min_modelsPack, list_models[[i]], sim_parameters) 
    
  })
  
}

# store results
results_min_MCP_nsample <- data.table(
  N_sample = c("(40,20,20,20,40)","(30,30,30,30,30)", "(48,24,24,24,48)","(36,36,36,36,36)"),
  Linear = c(results_list_min_nsample[[1]]$sim_results$power, results_list_min_nsample[[2]]$sim_results$power, results_list_min_nsample[[3]]$sim_results$power, results_list_min_nsample[[4]]$sim_results$power),
  Exponential = c(results_list_min_nsample[[5]]$sim_results$power, results_list_min_nsample[[6]]$sim_results$power, results_list_min_nsample[[7]]$sim_results$power, results_list_min_nsample[[8]]$sim_results$power),
  Emax = c(results_list_min_nsample[[9]]$sim_results$power, results_list_min_nsample[[10]]$sim_results$power, results_list_min_nsample[[11]]$sim_results$power, results_list_min_nsample[[12]]$sim_results$power),
  Average = c(sum(results_list_min_nsample[[1]]$sim_results$power, results_list_min_nsample[[5]]$sim_results$power, results_list_min_nsample[[9]]$sim_results$power)/3,
              sum(results_list_min_nsample[[2]]$sim_results$power, results_list_min_nsample[[6]]$sim_results$power, results_list_min_nsample[[10]]$sim_results$power)/3,
              sum(results_list_min_nsample[[3]]$sim_results$power, results_list_min_nsample[[7]]$sim_results$power, results_list_min_nsample[[11]]$sim_results$power)/3,
              sum(results_list_min_nsample[[4]]$sim_results$power, results_list_min_nsample[[8]]$sim_results$power, results_list_min_nsample[[12]]$sim_results$power)/3)
)

kable(results_min_MCP_nsample)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different sample sizes " = 5), font_size = 15, bold = TRUE)%>%
    add_header_above(c("Minimal scenario " = 5), font_size = 15, bold = TRUE)
Minimal scenario
Power results different sample sizes
N_sample Linear Exponential Emax Average
(40,20,20,20,40) 0.67 0.71 0.70 0.6933333
(30,30,30,30,30) 0.68 0.72 0.65 0.6833333
(48,24,24,24,48) 0.77 0.89 0.80 0.8200000
(36,36,36,36,36) 0.69 0.78 0.72 0.7300000

2.2 Monotonic scenario

2.2.1 varying expected effect for maximum dose

setup
# assumed dose - response models
sim_models_monotonic_linear = list(linear = NA,
                  max_effect = c(0.05,0.1,0.2,0.3,0.5),
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)

sim_models_monotonic_exp = list( exponential = 4.447149,
                  max_effect = c(0.05,0.1,0.2,0.3,0.5),
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)

sim_models_monotonic_logistic = list(logistic = c(1.5, 0.2275598),
                  max_effect = c(0.05,0.1,0.2,0.3,0.5),
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)

sim_models_monotonic_sigemax= list(sigemax =  c(1.528629, 4.087463),
                  max_effect = c(0.05,0.1,0.2,0.3,0.5),
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)

sim_models_monotonic_emax= list(emax = 0.6666667,
                  max_effect = c(0.05,0.1,0.2,0.3,0.5),
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)

# Simulation parameters
sim_parameters = list(n = Nsample, 
                      doses = doses.sim,
                      dropout_rate = 0.0, 
                      go_threshold = 0.1, 
                      nsims = n_sim)

# Initialize list to store results
results_list_monotonic <- list()

list_models <- list(sim_models_monotonic_linear, sim_models_monotonic_exp, sim_models_monotonic_emax, sim_models_monotonic_logistic, sim_models_monotonic_sigemax)
chunks <- chunkVector(seq_along(list_models), getDoParWorkers())

results_list_monotonic <- foreach(k = chunks, .combine = c) %dorng% {
  
  lapply(k, function (i) {
    
    func_sim(monotonic_modelsPack, list_models[[i]], sim_parameters) 
    
  })
  
}
  

results_monotonic_MCP <- data.table(
  max_eff = c(0.05, 0.1, 0.2, 0.3, 0.5),
  linear = c(results_list_monotonic[[1]]$sim_results$power),
  exp = c(results_list_monotonic[[2]]$sim_results$power),
  emax = c(results_list_monotonic[[3]]$sim_results$power),
  logistic = c(results_list_monotonic[[4]]$sim_results$power),
  sigEmax = c(results_list_monotonic[[5]]$sim_results$power),
  Average = c((sum(results_list_monotonic[[1]]$sim_results$power[1], results_list_monotonic[[2]]$sim_results$power[1], results_list_monotonic[[3]]$sim_results$power[1], results_list_monotonic[[4]]$sim_results$power[1], results_list_monotonic[[5]]$sim_results$power[1])/5),
  (sum(results_list_monotonic[[1]]$sim_results$power[2], results_list_monotonic[[2]]$sim_results$power[2], results_list_monotonic[[3]]$sim_results$power[2], results_list_monotonic[[4]]$sim_results$power[2], results_list_monotonic[[5]]$sim_results$power[2])/5),
              (sum(results_list_monotonic[[1]]$sim_results$power[3], results_list_monotonic[[2]]$sim_results$power[3], results_list_monotonic[[3]]$sim_results$power[3], results_list_monotonic[[4]]$sim_results$power[3], results_list_monotonic[[5]]$sim_results$power[3])/5),
(sum(results_list_monotonic[[1]]$sim_results$power[4], results_list_monotonic[[2]]$sim_results$power[4], results_list_monotonic[[3]]$sim_results$power[4], results_list_monotonic[[4]]$sim_results$power[4], results_list_monotonic[[5]]$sim_results$power[4])/5),
(sum(results_list_monotonic[[1]]$sim_results$power[5], results_list_monotonic[[2]]$sim_results$power[5], results_list_monotonic[[3]]$sim_results$power[5], results_list_monotonic[[4]]$sim_results$power[5], results_list_monotonic[[5]]$sim_results$power[5])/5))
)

kable(results_monotonic_MCP)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different expected effects " = 7), font_size = 15, bold = TRUE)%>%
    add_header_above(c("Monotonic scenario " = 7), font_size = 15, bold = TRUE)
Monotonic scenario
Power results different expected effects
max_eff linear exp emax logistic sigEmax Average
0.05 0.23 0.08 0.11 0.22 0.13 0.154
0.10 0.35 0.28 0.32 0.34 0.31 0.320
0.20 0.74 0.78 0.76 0.92 0.92 0.824
0.30 0.98 0.98 0.97 1.00 1.00 0.986
0.50 1.00 1.00 1.00 1.00 1.00 1.000

2.2.2 varying sample size

setup
# assumed dose- response model with assumed maximum effect = 0.2
sim_models_monotonic_linear$max_effect <- expectedEffect_fix
sim_models_monotonic_exp$max_effect <- expectedEffect_fix
sim_models_monotonic_emax$max_effect <- expectedEffect_fix
sim_models_monotonic_logistic$max_effect <- expectedEffect_fix
sim_models_monotonic_sigemax$max_effect <- expectedEffect_fix

# store results
power_monotonic_nsample_linear <- vector()
power_monotonic_nsample_exp <- vector()
power_monotonic_nsample_emax <- vector()
power_monotonic_nsample_logistic <- vector()
power_monotonic_nsample_sigemax <- vector()



# Set max_effect to 0.2
monotonic_modelsPack$max_effect <- expectedEffect_fix
#list of different samole sizes
nsample_list = list(20*c(2,1,1,1,2),30*c(1,1,1,1,1), 24*c(2,1,1,1,2),36*c(1,1,1,1,1),
                    20*c(2,1,1,1,2),30*c(1,1,1,1,1), 24*c(2,1,1,1,2),36*c(1,1,1,1,1),
                    20*c(2,1,1,1,2),30*c(1,1,1,1,1), 24*c(2,1,1,1,2),36*c(1,1,1,1,1),
                    20*c(2,1,1,1,2),30*c(1,1,1,1,1), 24*c(2,1,1,1,2),36*c(1,1,1,1,1),
                    20*c(2,1,1,1,2),30*c(1,1,1,1,1), 24*c(2,1,1,1,2),36*c(1,1,1,1,1))

list_models <- list(sim_models_monotonic_linear, sim_models_monotonic_linear, sim_models_monotonic_linear, sim_models_monotonic_linear,
                    sim_models_monotonic_exp, sim_models_monotonic_exp, sim_models_monotonic_exp, sim_models_monotonic_exp,
                    sim_models_monotonic_emax, sim_models_monotonic_emax, sim_models_monotonic_emax, sim_models_monotonic_emax,
                    sim_models_monotonic_logistic, sim_models_monotonic_logistic, sim_models_monotonic_logistic, sim_models_monotonic_logistic,
                    sim_models_monotonic_sigemax, sim_models_monotonic_sigemax, sim_models_monotonic_sigemax, sim_models_monotonic_sigemax)


chunks <- chunkVector(seq_along(list_models), getDoParWorkers())

results_list_monotonic_nsample <- foreach(k = chunks, .combine = c, .export = c(as.character(monotonic_modelsPack), as.character(nsample_list)))  %dorng% {
  lapply(k, function (i) {
    sim_parameters = list(n = nsample_list[[i]], 
                      doses = doses.sim,
                      dropout_rate = 0.0, 
                      go_threshold = 0.1, 
                      nsims = n_sim)
    
    func_sim(monotonic_modelsPack, list_models[[i]], sim_parameters) 
    
  })
  
}


# store results
results_monotonic_MCP_nsample <- data.table(
  N_sample = c("(40,20,20,20,40)","(30,30,30,30,30)", "(48,24,24,24,48)","(36,36,36,36,36)"),
  Linear = c(results_list_monotonic_nsample[[1]]$sim_results$power, results_list_monotonic_nsample[[2]]$sim_results$power, results_list_monotonic_nsample[[3]]$sim_results$power, results_list_monotonic_nsample[[4]]$sim_results$power),
  Exponential = c(results_list_monotonic_nsample[[5]]$sim_results$power, results_list_monotonic_nsample[[6]]$sim_results$power, results_list_monotonic_nsample[[7]]$sim_results$power, results_list_monotonic_nsample[[8]]$sim_results$power),
  Emax = c(results_list_monotonic_nsample[[9]]$sim_results$power, results_list_monotonic_nsample[[10]]$sim_results$power, results_list_monotonic_nsample[[11]]$sim_results$power, results_list_monotonic_nsample[[12]]$sim_results$power),
  Logistic = c(results_list_monotonic_nsample[[13]]$sim_results$power, results_list_monotonic_nsample[[14]]$sim_results$power, results_list_monotonic_nsample[[15]]$sim_results$power, results_list_monotonic_nsample[[16]]$sim_results$power),
  sigEmax = c(results_list_monotonic_nsample[[17]]$sim_results$power, results_list_monotonic_nsample[[18]]$sim_results$power, results_list_monotonic_nsample[[19]]$sim_results$power, results_list_monotonic_nsample[[20]]$sim_results$power),
Average = c((sum(results_list_monotonic_nsample[[1]]$sim_results$power, results_list_monotonic_nsample[[5]]$sim_results$power, results_list_monotonic_nsample[[9]]$sim_results$power, ... =   results_list_monotonic_nsample[[13]]$sim_results$power, results_list_monotonic_nsample[[17]]$sim_results$power)/5),
              (sum(results_list_monotonic_nsample[[2]]$sim_results$power, results_list_monotonic_nsample[[6]]$sim_results$power, results_list_monotonic_nsample[[10]]$sim_results$power, results_list_monotonic_nsample[[14]]$sim_results$power, results_list_monotonic_nsample[[18]]$sim_results$power)/5),
            (sum(results_list_monotonic_nsample[[3]]$sim_results$power, results_list_monotonic_nsample[[7]]$sim_results$power, results_list_monotonic_nsample[[11]]$sim_results$power, results_list_monotonic_nsample[[15]]$sim_results$power, results_list_monotonic_nsample[[19]]$sim_results$power)/5),
            (sum(results_list_monotonic_nsample[[4]]$sim_results$power, results_list_monotonic_nsample[[8]]$sim_results$power, results_list_monotonic_nsample[[12]]$sim_results$power, results_list_monotonic_nsample[[16]]$sim_results$power, results_list_monotonic_nsample[[20]]$sim_results$power)/5)))


kable(results_monotonic_MCP_nsample)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different sample sizes" = 7), font_size = 15, bold = TRUE)%>%
    add_header_above(c("Monotonic scenario" = 7), font_size = 15, bold = TRUE)
Monotonic scenario
Power results different sample sizes
N_sample Linear Exponential Emax Logistic sigEmax Average
(40,20,20,20,40) 0.73 0.74 0.80 0.81 0.76 0.768
(30,30,30,30,30) 0.57 0.66 0.57 0.84 0.71 0.670
(48,24,24,24,48) 0.76 0.79 0.83 0.91 0.89 0.836
(36,36,36,36,36) 0.58 0.71 0.74 0.90 0.80 0.746

2.3 Non-monotonic scenario

2.3.1 varying expected effect for maximum dose

For the simulations of the non-monotonic scenario, the R package ‘DoseFinding’ was used instead of ‘MCPModPack’. In particular the powMCT function was utilized to calculate success probabilities for the various scenarios.

setup
# linear with DoseFinding package
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(linear = NULL)
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.05,0.1,0.2,0.3,0.5), sd = c(0.4), n_total = c(200))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_linear <- result$linear
setup
# emax with DoseFinding package
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(emax = 0.6666667)
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.05,0.1,0.2,0.3,0.5), sd = c(0.4), n_total = c(200))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_emax <- result$`emax (ED50=0.6666667)`
setup
# sigemax with DoseFinding package
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(sigEmax = rbind(c(1.528629, 4.087463)))
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.05,0.1,0.2,0.3,0.5), sd = c(0.4), n_total = c(200))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_sigemax <- result$`sigEmax (ED50=1.528629,h=4.087463)`
setup
#quadratic with DoseFinding package
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(quadratic = -0.09688711)
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.05,0.1,0.2,0.3,0.5), sd = c(0.4), n_total = c(200))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_quadratic <- result$`quadratic (delta=-0.09688711)`
setup
# beta with DoseFinding package
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(betaMod = rbind(c(1.396434, 1.040978)))
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.05,0.1,0.2,0.3,0.5), sd = c(0.4), n_total = c(200))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_beta <- result$`betaMod (delta1=1.396434,delta2=1.040978,scal=9.6)`
setup
#linear, emax, sigemax, quadratic
results_non_monotonic_MCP <- data.table(
  Max_eff = c(0.05, 0.1, 0.2, 0.3, 0.5),
  linear = power_non_monotonic_linear,
  emax = power_non_monotonic_emax,
  sigEmax = power_non_monotonic_sigemax,
  quadratic = power_non_monotonic_quadratic,
  beta = power_non_monotonic_beta,
  average = c((sum(power_non_monotonic_linear[1],power_non_monotonic_emax[1], power_non_monotonic_sigemax[1], power_non_monotonic_quadratic[1], power_non_monotonic_beta[1])/5),
    (sum(power_non_monotonic_linear[2],power_non_monotonic_emax[2], power_non_monotonic_sigemax[2], power_non_monotonic_quadratic[2], power_non_monotonic_beta[2])/5),
    (sum(power_non_monotonic_linear[3],power_non_monotonic_emax[3], power_non_monotonic_sigemax[3], power_non_monotonic_quadratic[3], power_non_monotonic_beta[3])/5),
    (sum(power_non_monotonic_linear[4],power_non_monotonic_emax[4], power_non_monotonic_sigemax[4], power_non_monotonic_quadratic[4], power_non_monotonic_beta[4])/5),
    (sum(power_non_monotonic_linear[5],power_non_monotonic_emax[5], power_non_monotonic_sigemax[5], power_non_monotonic_quadratic[5], power_non_monotonic_beta[5])/5))
  )


kable(results_non_monotonic_MCP)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different expected effects" = 7), font_size = 15, bold = TRUE)%>%
    add_header_above(c("Non-monotonic scenario" = 7), font_size = 15, bold = TRUE)
Non-monotonic scenario
Power results different expected effects
Max_eff linear emax sigEmax quadratic beta average
0.05 0.1361826 0.1438824 0.1720734 0.1306460 0.1272818 0.1420132
0.10 0.2985730 0.3174089 0.4040673 0.2768840 0.2649841 0.3123835
0.20 0.7408952 0.7645582 0.8819700 0.6924798 0.6605961 0.7480998
0.30 0.9677230 0.9740876 0.9957499 0.9464217 0.9301074 0.9628180
0.50 0.9999913 0.9999952 1.0000000 0.9999419 0.9998595 0.9999576

2.3.2 varying sample size

setup
# store results
power_non_monotonic_nsample_linear <- vector()
power_non_monotonic_nsample_emax <- vector()
power_non_monotonic_nsample_quadratic <- vector()
power_non_monotonic_nsample_sigemax <- vector()
power_non_monotonic_nsample_beta <- vector()
setup
# linear with DoseFinding
# allocation = c(1,1,1,1,1)
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(linear = NULL)
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(150,180))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_nsample_linear[2] <- result$linear[1]
power_non_monotonic_nsample_linear[4] <- result$linear[2]

# allocation = c(2,1,1,1,2)
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(2,1,1,1,2)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(linear = NULL)
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(140,168))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_nsample_linear[1] <- result$linear[1]
power_non_monotonic_nsample_linear[3] <- result$linear[2]
setup
# emax with DoseFinding
#allocation = c(1,1,1,1,1)
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(emax = 0.6666667)
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(150,180))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_nsample_emax[2] <- result$`emax (ED50=0.6666667)`[1]
power_non_monotonic_nsample_emax[4] <- result$`emax (ED50=0.6666667)`[2]

# allocation = c(2,1,1,1,2)
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(2,1,1,1,2)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(emax = 0.6666667)
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(140,168))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_nsample_emax[1] <- result$`emax (ED50=0.6666667)`[1]
power_non_monotonic_nsample_emax[3] <- result$`emax (ED50=0.6666667)`[2]
setup
# sigemax with DoseFinding
#allocation = c(1,1,1,1,1)
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(sigEmax = rbind(c(1.528629, 4.087463)))
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(150,180))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_nsample_sigemax[2] <- result$`sigEmax (ED50=1.528629,h=4.087463)`[1]
power_non_monotonic_nsample_sigemax[4] <- result$`sigEmax (ED50=1.528629,h=4.087463)`[2]

# allocation = c(2,1,1,1,2)
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(2,1,1,1,2)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(sigEmax = rbind(c(1.528629, 4.087463)))
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(140,168))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_nsample_sigemax[1] <- result$`sigEmax (ED50=1.528629,h=4.087463)`[1]
power_non_monotonic_nsample_sigemax[3] <- result$`sigEmax (ED50=1.528629,h=4.087463)`[2]
setup
# quadratic with DoseFinding
#allocation = c(1,1,1,1,1)
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(quadratic = -0.09688711)
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(150,180))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_nsample_quadratic[2] <- result$`quadratic (delta=-0.09688711)`[1]
power_non_monotonic_nsample_quadratic[4] <- result$`quadratic (delta=-0.09688711)`[2]

# allocation = c(2,1,1,1,2)
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(2,1,1,1,2)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(quadratic = -0.09688711)
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(140,168))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_nsample_quadratic[1] <- result$`quadratic (delta=-0.09688711)`[1]
power_non_monotonic_nsample_quadratic[3] <- result$`quadratic (delta=-0.09688711)`[2]
setup
# beta with DoseFinding
#allocation = c(1,1,1,1,1)
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(betaMod = rbind(c(1.396434, 1.040978)))
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(150,180))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_nsample_beta[2] <- result$`betaMod (delta1=1.396434,delta2=1.040978,scal=9.6)`[1]
power_non_monotonic_nsample_beta[4] <- result$`betaMod (delta1=1.396434,delta2=1.040978,scal=9.6)`[2]

# allocation = c(2,1,1,1,2)
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(2,1,1,1,2)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(betaMod = rbind(c(1.396434, 1.040978)))
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(140,168))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
power_non_monotonic_nsample_beta[1] <- result$`betaMod (delta1=1.396434,delta2=1.040978,scal=9.6)`[1]
power_non_monotonic_nsample_beta[3] <- result$`betaMod (delta1=1.396434,delta2=1.040978,scal=9.6)`[2]
setup
#store results
results_non_monotonic_MCP_nsample <- data.table(
  N_sample = c("(40,20,20,20,40)","(30,30,30,30,30)", "(48,24,24,24,48)","(36,36,36,36,36)"),
  Linear = power_non_monotonic_nsample_linear,
  Emax = power_non_monotonic_nsample_emax,
  sigEmax = power_non_monotonic_nsample_sigemax,
  quadratic = power_non_monotonic_nsample_quadratic,
  beta = power_non_monotonic_nsample_beta
)


kable(results_non_monotonic_MCP_nsample)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different sample sizes" = 6), font_size = 15, bold = TRUE)%>%
    add_header_above(c("Non-monotonic scenario" = 6), font_size = 15, bold = TRUE)
Non-monotonic scenario
Power results different sample sizes
N_sample Linear Emax sigEmax quadratic beta
(40,20,20,20,40) 0.7182897 0.7452513 0.8118555 0.5991157 0.5492169
(30,30,30,30,30) 0.6265819 0.6539679 0.7872977 0.5801699 0.5521925
(48,24,24,24,48) 0.7869752 0.8122048 0.8716953 0.6689755 0.6177505
(36,36,36,36,36) 0.6994504 0.7248996 0.8504362 0.6509720 0.6195153

2.4 Variability scenario

2.4.1 varying expected effect for maximum dose

setup
# assumed dose - response model
sim_models_var_linear = list(linear = NA,
                  max_effect = c(0.05,0.1,0.2,0.3,0.5),
                  sd = rep(sd.sim_var, length(doses_var)),
                  placebo_effect = plc.guess)

sim_models_var_exp = list( exponential = 4.447149,
                  max_effect = c(0.05,0.1,0.2,0.3,0.5),
                  sd = rep(sd.sim_var, length(doses_var)),
                  placebo_effect = plc.guess)

sim_models_var_emax = list(emax = 0.6666667,
                  max_effect = c(0.05,0.1,0.2,0.3,0.5),
                  sd = rep(sd.sim_var, length(doses_var)),
                  placebo_effect = plc.guess)



# Simulation parameters
sim_parameters_var = list(n = Nsample_var, 
                      doses = doses_var,
                      dropout_rate = 0.0, 
                      go_threshold = 0.1, 
                      nsims = n_sim)

# perform Simulations

list_models <- list(sim_models_var_linear, sim_models_var_exp, sim_models_var_emax)
chunks <- chunkVector(seq_along(list_models), getDoParWorkers())

results_list_var <- foreach(k = chunks, .combine = c) %dorng% {
  
  lapply(k, function (i) {
    
    func_sim(var_modelsPack, list_models[[i]], sim_parameters_var) 
    
  })
  
}
  

# store results
results_var_MCP <- data.table(
  Max_eff = c(0.05, 0.1, 0.2, 0.3, 0.5),
  linear = results_list_var[[1]]$sim_results$power,
  exponential = results_list_var[[2]]$sim_results$power,
  emax = results_list_var[[3]]$sim_results$power,
  average = c((sum(results_list_var[[1]]$sim_results$power[1], results_list_var[[2]]$sim_results$power[1],results_list_var[[3]]$sim_results$power[1])/3),
              (sum(results_list_var[[1]]$sim_results$power[2], results_list_var[[2]]$sim_results$power[2],results_list_var[[3]]$sim_results$power[2])/3),
              (sum(results_list_var[[1]]$sim_results$power[3], results_list_var[[2]]$sim_results$power[3],results_list_var[[3]]$sim_results$power[3])/3),
              (sum(results_list_var[[1]]$sim_results$power[4], results_list_var[[2]]$sim_results$power[4],results_list_var[[3]]$sim_results$power[4])/3),
              (sum(results_list_var[[1]]$sim_results$power[5], results_list_var[[2]]$sim_results$power[5],results_list_var[[3]]$sim_results$power[5])/3))
)

kable(results_var_MCP)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different expected effects" = 5), font_size = 15, bold = TRUE)%>%
    add_header_above(c("Variability scenario" = 5), font_size = 15, bold = TRUE)
Variability scenario
Power results different expected effects
Max_eff linear exponential emax average
0.05 0.12 0.12 0.12 0.1200000
0.10 0.19 0.28 0.20 0.2233333
0.20 0.60 0.65 0.52 0.5900000
0.30 0.82 0.87 0.84 0.8433333
0.50 1.00 1.00 1.00 1.0000000

2.4.2 varying sample size

setup
# assume maximum effect = 0.2
sim_models_var_linear$max_effect <- expectedEffect_fix
sim_models_var_exp$max_effect <- expectedEffect_fix
sim_models_var_emax$max_effect <- expectedEffect_fix


# store results
power_var_nsample_linear <- vector()
power_var_nsample_exp <- vector()
power_var_nsample_emax <- vector()

# # Loop over nsample values
# for (i in seq_along(nsample_list_var)) {
#   
#   # Update nsample in simulation parameters
#   sim_parameters_var$n <- nsample_list_var[[i]]
# 
#   
#   # power values for different assumed true models
#   power_nsample_linear <- func_sim(var_modelsPack, sim_models_var_linear, sim_parameters_var)
#   power_var_nsample_linear[i] <- power_nsample_linear$sim_results$power
#   
#   power_nsample_exp <- func_sim(var_modelsPack, sim_models_var_exp, sim_parameters_var)
#   power_var_nsample_exp[i] <- power_nsample_exp$sim_results$power
#   
#   power_nsample_emax <- func_sim(var_modelsPack, sim_models_var_emax, sim_parameters_var)
#   power_var_nsample_emax[i] <-power_nsample_emax$sim_results$power
#   
#   
# }

var_modelsPack$max_effect <- expectedEffect_fix
#list of different samole sizes
nsample_list_var = list(20*c(2,1,1,2),30*c(1,1,1,1), 24*c(2,1,1,2),36*c(1,1,1,1),
                    20*c(2,1,1,2),30*c(1,1,1,1), 24*c(2,1,1,2),36*c(1,1,1,1),
                    20*c(2,1,1,2),30*c(1,1,1,1), 24*c(2,1,1,2),36*c(1,1,1,1))

list_models <- list(sim_models_var_linear, sim_models_var_linear, sim_models_var_linear, sim_models_var_linear, 
                    sim_models_var_exp, sim_models_var_exp, sim_models_var_exp, sim_models_var_exp, 
                    sim_models_var_emax, sim_models_var_emax, sim_models_var_emax, sim_models_var_emax)


chunks <- chunkVector(seq_along(list_models), getDoParWorkers())

results_list_var_nsample <- foreach(k = chunks, .combine = c, .export = c(as.character(var_modelsPack), as.character(nsample_list)))  %dorng% {
  lapply(k, function (i) {
    sim_parameters = list(n = nsample_list_var[[i]], 
                      doses = doses_var,
                      dropout_rate = 0.0, 
                      go_threshold = 0.1, 
                      nsims = n_sim)
    
    
    func_sim(var_modelsPack, list_models[[i]], sim_parameters) 
    
  })
  
}


results_var_MCP_nsample <- data.table(
  N_sample = c("(40,20,20,20,40)","(30,30,30,30,30)", "(48,24,24,24,48)","(36,36,36,36,36)"),
  Linear = c(results_list_var_nsample[[1]]$sim_results$power, results_list_var_nsample[[2]]$sim_results$power, results_list_var_nsample[[3]]$sim_results$power, results_list_var_nsample[[4]]$sim_results$power),
  Exponential =c(results_list_var_nsample[[5]]$sim_results$power, results_list_var_nsample[[6]]$sim_results$power, results_list_var_nsample[[7]]$sim_results$power, results_list_var_nsample[[8]]$sim_results$power),
  Emax = c(results_list_var_nsample[[9]]$sim_results$power, results_list_var_nsample[[10]]$sim_results$power, results_list_var_nsample[[11]]$sim_results$power, results_list_var_nsample[[12]]$sim_results$power),
  Average = c((sum(results_list_var_nsample[[1]]$sim_results$power, results_list_var_nsample[[5]]$sim_results$power, results_list_var_nsample[[9]]$sim_results$power)/3),
              (sum(results_list_var_nsample[[2]]$sim_results$power, results_list_var_nsample[[6]]$sim_results$power, results_list_var_nsample[[10]]$sim_results$power)/3),
              (sum(results_list_var_nsample[[3]]$sim_results$power, results_list_var_nsample[[7]]$sim_results$power, results_list_var_nsample[[11]]$sim_results$power)/3),
              (sum(results_list_var_nsample[[4]]$sim_results$power, results_list_var_nsample[[8]]$sim_results$power, results_list_var_nsample[[12]]$sim_results$power)/3))
)

kable(results_var_MCP_nsample)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different sample sizes" = 5), font_size = 15, bold = TRUE)%>%
    add_header_above(c("Variability scenario" = 5), font_size = 15, bold = TRUE)
Variability scenario
Power results different sample sizes
N_sample Linear Exponential Emax Average
(40,20,20,20,40) 0.83 0.71 0.69 0.7433333
(30,30,30,30,30) 0.64 0.62 0.68 0.6466667
(48,24,24,24,48) 0.81 0.76 0.85 0.8066667
(36,36,36,36,36) 0.76 0.60 0.62 0.6600000

3 BayesianMCPMod

Following simulations will be conducted utilizing the ‘BayesianMCPMod’ package using the same scenarios as for the frequentist evaluations.

4 Prior Specification

In a first step an uninformative prior is specified.

setup
uninf_prior_list <- list(
                   Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, n = 1), sigma = sd.sim, param = "mn"),
                   DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 0, n = 1), sigma = sd.sim, param = "mn"),
                   DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 0, n = 1), sigma = sd.sim, param = "mn"),
                   DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 0, n = 1), sigma = sd.sim, param = "mn"),
                   DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 0, n = 1), sigma = sd.sim, param = "mn")
                   
)

To calculate success probabilities for the different assumed dose-response models and the specified trial design we will apply the assessDesign function.

4.1 Minimal scenario

4.1.1 varying expected effect for maximum dose

setup
# list to store results
results_list_Bay_min <- list()


list_max_eff <- as.list(c(0.05,0.1,0.2,0.3,0.5))
chunks <- chunkVector(seq_along(list_max_eff), getDoParWorkers())


  
results_list_Bay_min <- foreach(k = chunks, .combine = c) %dorng% {
  
  lapply(k, function (i) {
    
     min_models <- Mods(linear=NULL, 
                     exponential = exp.g,
                     emax = emax.g,
                     doses = doses.sim,
                     placEff = plc.guess, 
                     maxEff = list_max_eff[[i]],
                     direction = "increasing")
  #optimal contrasts
contM <- getContr(mods = min_models, 
                    dose_levels = doses.sim, 
                    prior_list = uninf_prior_list, 
                    dose_weights = c(1,1,1,1,1))

 # Simulation step
  success_probabilities_min <- assessDesign(
    n_patients  = Nsample,
    mods        = min_models,
    prior_list  = uninf_prior_list,
    sd          = sd.sim,
    n_sim       = n_sim,
    alpha_crit_val = alpha, 
    contr = contM) 
    
  })
  
}
  
results_min_Bay <- extract_success_rates(results_list_Bay_min, min_scenario)

minimal_Bay <- print_result_Bay_max_eff(results_min_Bay, c(min_scenario, "average"), expectedEffect)

minimal_Bay$kable_result
BayesianMCPMod
Power results different expected effects
variable linear exponential emax average
0.05 0.09 0.09 0.05 0.0766667
0.10 0.28 0.26 0.30 0.2800000
0.20 0.80 0.81 0.78 0.7966667
0.30 0.97 0.97 0.99 0.9766667
0.50 1.00 1.00 1.00 1.0000000

4.1.2 varying sample size

setup
# Define the list of sample sizes
nsample_list = list(20*c(2,1,1,1,2),30*c(1,1,1,1,1), 24*c(2,1,1,1,2),36*c(1,1,1,1,1))

#optimal contrasts
#allocation c(1,1,1,1,1)
contM <- getContr(mods = min_models, 
                  dose_levels = doses.sim, 
                  prior_list = uninf_prior_list, 
                  dose_weights = c(1,1,1,1,1))

#allocation c(2,1,1,1,2)
contM_2 <- getContr(mods = min_models, 
                  dose_levels = doses.sim, 
                  prior_list = uninf_prior_list, 
                  dose_weights = c(2,1,1,1,2))
contrasts_list = list(contM_2, contM, contM_2, contM)

# Initialize list to store results
results_list_nsample_min_Bay <- list()

#Models with maxEff = 0.2
min_models <- Mods(linear=NULL, 
                     exponential = exp.g,
                     emax = emax.g,
                     doses = doses.sim,
                     placEff = plc.guess, 
                     maxEff = expectedEffect_fix,
                     direction = "increasing")
  


chunks <- chunkVector(seq_along(nsample_list), getDoParWorkers())
  
results_list_nsample_min_Bay <- foreach(k = chunks, .combine = c, .export = c(as.character(min_models), as.character(contrasts_list)))  %dorng% {
  
  lapply(k, function (i) {
    
      success_probabilities_min <- assessDesign(
    n_patients  = nsample_list[[i]],
    mods        = min_models,
    prior_list  = uninf_prior_list,
    sd          = sd.sim,
    n_sim       = n_sim,
    alpha_crit_val = alpha, 
    contr = contrasts_list[[i]])
    
  })
}
  

results_min_Bay_nsample <- extract_success_rates_nsample(results_list_nsample_min_Bay, min_scenario)
  
minimal_nsample_Bay <- print_result_Bay_nsample(results_min_Bay_nsample, c(min_scenario, "average"), nsample_vector)
minimal_nsample_Bay$kable_result
BayesianMCPMod
Success probability results different sample sizes
variable linear exponential emax average
(40,20,20,20,40) 0.67 0.73 0.67 0.6900000
(30,30,30,30,30) 0.66 0.64 0.66 0.6533333
(48,24,24,24,48) 0.81 0.81 0.81 0.8100000
(36,36,36,36,36) 0.67 0.70 0.73 0.7000000

4.2 Monotonic scenario

4.2.1 varying expected effect for maximum dose

setup
results_list_Bay_monotonic <- list()




list_max_eff <- as.list(c(0.05,0.1,0.2,0.3,0.5))
chunks <- chunkVector(seq_along(list_max_eff), getDoParWorkers())


  
results_list_Bay_monotonic <- foreach(k = chunks, .combine = c) %dorng% {
  
  lapply(k, function (i) {
    
    monotonic_models <- Mods(linear=NULL, 
                         exponential = exp.g,
                         emax = emax.g,
                         logistic=logit.g,
                         sigEmax = sigEmax.g,
                         doses = doses.sim,
                         placEff = plc.guess, 
                         maxEff = list_max_eff[[i]],
                         direction = "increasing")
       
  # optimal contrasts
  contM <- getContr(mods = monotonic_models, 
                    dose_levels = doses.sim, 
                    prior_list = uninf_prior_list, 
                    dose_weights = c(1,1,1,1,1))
  
  # perform Simulations
  success_probabilities_monotonic <- assessDesign(
    n_patients  = Nsample,
    mods        = monotonic_models,
    prior_list  = uninf_prior_list,
    sd          = sd.sim,
    n_sim       = n_sim,
    alpha_crit_val = alpha, 
    contr = contM)
    

    
  })
  
}
  

##power results in vector for tables
results_monotonic_Bay <- extract_success_rates(results_list_Bay_monotonic, monotonic_scenario)


monotonic_Bay <- print_result_Bay_max_eff(results_monotonic_Bay, c(monotonic_scenario, "average"), expectedEffect)
  

monotonic_Bay$kable_result
BayesianMCPMod
Power results different expected effects
variable linear exponential emax logistic sigEmax average
0.05 0.15 0.15 0.15 0.16 0.15 0.152
0.10 0.30 0.31 0.32 0.37 0.35 0.330
0.20 0.78 0.79 0.79 0.91 0.87 0.828
0.30 0.99 0.97 1.00 1.00 1.00 0.992
0.50 1.00 1.00 1.00 1.00 1.00 1.000

4.2.2 varying sample size

setup
# Initialize list to store results
results_list_nsample_monotonic_Bay <- list()

# dose - response models
monotonic_models <- Mods(linear=NULL, 
                         exponential = exp.g,
                         emax = emax.g,
                         logistic=logit.g,
                         sigEmax = sigEmax.g,
                         doses = doses.sim,
                         placEff = plc.guess, 
                         maxEff = expectedEffect_fix,
                         direction = "increasing")
 
#optimal contrasts
#allocation c(1,1,1,1,1)
contM <- getContr(mods = monotonic_models, 
                  dose_levels = doses.sim, 
                  prior_list = uninf_prior_list, 
                  dose_weights = c(1,1,1,1,1))

#allocation c(2,1,1,1,2)
contM_2 <- getContr(mods = monotonic_models, 
                  dose_levels = doses.sim, 
                  prior_list = uninf_prior_list, 
                  dose_weights = c(2,1,1,1,2))
contrasts_list = list(contM_2, contM, contM_2, contM)


chunks <- chunkVector(seq_along(nsample_list), getDoParWorkers())
  
results_list_nsample_monotonic_Bay <- foreach(k = chunks, .combine = c, .export = c(as.character(monotonic_models), as.character(contrasts_list)))  %dorng% {
  
  lapply(k, function (i) {
    
    success_probabilities_monotonic <- assessDesign(
    n_patients  = nsample_list[[i]],
    mods        = monotonic_models,
    prior_list  = uninf_prior_list,
    sd          = sd.sim,
    n_sim       = n_sim,
    alpha_crit_val = alpha,
    contr = contrasts_list[[i]])
  })
}
  

results_monotonic_Bay_nsample <- extract_success_rates_nsample(results_list_nsample_monotonic_Bay, monotonic_scenario)


monotonic_nsample_Bay <- print_result_Bay_nsample(results_monotonic_Bay_nsample, c(monotonic_scenario, "average"), nsample_vector)
monotonic_nsample_Bay$kable_result
BayesianMCPMod
Success probability results different sample sizes
variable linear exponential emax logistic sigEmax average
(40,20,20,20,40) 0.80 0.80 0.77 0.88 0.84 0.818
(30,30,30,30,30) 0.65 0.70 0.63 0.79 0.75 0.704
(48,24,24,24,48) 0.84 0.85 0.78 0.87 0.86 0.840
(36,36,36,36,36) 0.77 0.76 0.74 0.91 0.85 0.806

4.3 Non-monotonic scenario

4.3.1 varying expected effect for maximum dose

setup
results_list_Bay_non_monotonic <- list()


list_max_eff <- as.list(c(0.05,0.1,0.2,0.3,0.5))
chunks <- chunkVector(seq_along(list_max_eff), getDoParWorkers())


  
results_list_Bay_non_monotonic <- foreach(k = chunks, .combine = c) %dorng% {
  
  lapply(k, function (i) {
    
    
     non_monotonic_models <- Mods(linear=NULL, 
                               emax = emax.g,
                               sigEmax = sigEmax.g,
                               quadratic = quad.g,
                               betaMod = beta.g,
                               doses = doses.sim,
                               placEff = plc.guess,
                               maxEff = list_max_eff[[i]],
                               direction = "increasing",
                               addArgs = list(scal=9.6))

  #optimal contrasts
  contM <- getContr(mods = non_monotonic_models,
                    dose_levels = doses.sim,
                    prior_list = uninf_prior_list,
                    dose_weights = c(1,1,1,1,1))

  success_probabilities_non_monotonic <- assessDesign(
    n_patients  = Nsample,
    mods        = non_monotonic_models,
    prior_list  = uninf_prior_list,
    sd          = sd.sim,
    n_sim       = n_sim,
    alpha_crit_val = alpha,
    contr = contM)


  })
  
}
  
results_non_monotonic_Bay <- extract_success_rates(results_list_Bay_non_monotonic, non_monotonic_scenario)

non_monotonic_Bay <- print_result_Bay_max_eff(results_non_monotonic_Bay, c(non_monotonic_scenario, "average"), expectedEffect)
  
non_monotonic_Bay$kable_result
BayesianMCPMod
Power results different expected effects
variable linear emax sigEmax quadratic betaMod average
0.05 0.10 0.09 0.12 0.09 0.09 0.098
0.10 0.27 0.32 0.41 0.28 0.27 0.310
0.20 0.75 0.79 0.87 0.70 0.66 0.754
0.30 0.96 0.95 1.00 0.91 0.90 0.944
0.50 1.00 1.00 1.00 1.00 1.00 1.000

4.3.2 varying sample size

setup
# Initialize list to store results
results_list_nsample_non_monotonic_Bay <- list()

non_monotonic_models <- Mods(linear=NULL, 
                             emax = emax.g,
                             sigEmax = sigEmax.g,
                             quadratic = quad.g,
                             betaMod = beta.g,
                             doses = doses.sim,
                             placEff = plc.guess, 
                             maxEff = expectedEffect_fix,
                             direction = "increasing",
                             addArgs = list(scal=9.6))
 
  #optimal contrasts
contM <- getContr(mods = non_monotonic_models, 
                  dose_levels = doses.sim, 
                  prior_list = uninf_prior_list, 
                  dose_weights = c(1,1,1,1,1))

contM_2 <- getContr(mods = non_monotonic_models, 
                  dose_levels = doses.sim, 
                  prior_list = uninf_prior_list, 
                  dose_weights = c(2,1,1,1,2))

contrasts_list = list(contM_2, contM, contM_2, contM)
  


chunks <- chunkVector(seq_along(nsample_list), getDoParWorkers())
  
results_list_nsample_non_monotonic_Bay <- foreach(k = chunks, .combine = c, .export = c(as.character(non_monotonic_models), as.character(contrasts_list)))  %dorng% {
  
  lapply(k, function (i) {
    
    success_probabilities_non_monotonic <- assessDesign(
    n_patients  = nsample_list[[i]],
    mods        = non_monotonic_models,
    prior_list  = uninf_prior_list,
    sd          = sd.sim,
    n_sim       = n_sim,
    alpha_crit_val = alpha,
    contr = contrasts_list[[i]])
  })
}
  

results_non_monotonic_Bay_nsample <- extract_success_rates_nsample(results_list_nsample_non_monotonic_Bay, non_monotonic_scenario)


non_monotonic_nsample_Bay <- print_result_Bay_nsample(results_non_monotonic_Bay_nsample, c(non_monotonic_scenario, "average"), nsample_vector)
non_monotonic_nsample_Bay$kable_result
BayesianMCPMod
Success probability results different sample sizes
variable linear emax sigEmax quadratic betaMod average
(40,20,20,20,40) 0.73 0.76 0.81 0.58 0.50 0.676
(30,30,30,30,30) 0.63 0.64 0.77 0.52 0.49 0.610
(48,24,24,24,48) 0.74 0.75 0.83 0.65 0.59 0.712
(36,36,36,36,36) 0.67 0.67 0.78 0.59 0.56 0.654

4.4 Variability scenario

4.4.1 varying expected effect for maximum dose

setup
# prior 
var_uninf_prior_list <- list(
                   Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, n = 1), sigma = sd.sim_var, param = "mn"),
                   DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 0, n = 1), sigma = sd.sim_var, param = "mn"),
                   DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 0, n = 1), sigma = sd.sim_var, param = "mn"),
                   DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 0, n = 1), sigma = sd.sim_var, param = "mn"))

results_list_Bay_var <- list()



list_max_eff <- as.list(c(0.05,0.1,0.2,0.3,0.5))
chunks <- chunkVector(seq_along(list_max_eff), getDoParWorkers())


  
results_list_Bay_var <- foreach(k = chunks, .combine = c) %dorng% {
  
  lapply(k, function (i) {
    
      var_models <- Mods(linear = NULL, 
                     exponential = exp.g,
                     emax = emax.g,
                     doses = doses_var,
                     placEff = plc.guess,
                     maxEff = list_max_eff[[i]],
                     direction = "increasing")

  contM <- getContr(mods = var_models,
                    dose_levels = doses_var,
                    prior_list = var_uninf_prior_list,
                    dose_weights = c(1,1,1,1))

  success_probabilities_var <- assessDesign(
    n_patients  = Nsample_var,
    mods        = var_models,
    prior_list  = var_uninf_prior_list,
    sd          = sd.sim_var,
    n_sim       = n_sim,
    alpha_crit_val = alpha,
    contr = contM)

  })
  
}


results_variability_Bay <- extract_success_rates(results_list_Bay_var, variability_scenario)

variability_Bay <- print_result_Bay_max_eff(results_variability_Bay, c(variability_scenario, "average"), expectedEffect)

variability_Bay$kable_result
BayesianMCPMod
Power results different expected effects
variable linear exponential emax average
0.05 0.12 0.12 0.09 0.1100000
0.10 0.22 0.21 0.23 0.2200000
0.20 0.57 0.56 0.58 0.5700000
0.30 0.83 0.86 0.87 0.8533333
0.50 0.99 1.00 0.99 0.9933333

4.4.2 varying sample size

setup
results_list_nsample_Bay_var <- list()

var_models <- Mods(linear=NULL,
                   exponential = exp.g,
                   emax = emax.g,
                   doses = doses_var,
                   placEff = plc.guess, 
                   maxEff = expectedEffect_fix,
                   direction = "increasing")

#optimal contrasts
contM <- getContr(mods = var_models, 
                  dose_levels = doses_var, 
                  prior_list = var_uninf_prior_list, 
                  dose_weights = c(1,1,1,1))

contM_2 <- getContr(mods = var_models, 
                  dose_levels = doses_var, 
                  prior_list = var_uninf_prior_list, 
                  dose_weights = c(2,1,1,2))
contrasts_list = list(contM_2, contM, contM_2, contM)


nsample_list_var <- list(20*c(2,1,1,2),30*c(1,1,1,1), 24*c(2,1,1,2),36*c(1,1,1,1))

chunks <- chunkVector(seq_along(nsample_list_var), getDoParWorkers())
  
results_list_nsample_Bay_var <- foreach(k = chunks, .combine = c, .export = c(as.character(var_models), as.character(contrasts_list)))  %dorng% {
  
  lapply(k, function (i) {
    success_probabilities_var <- assessDesign(
    n_patients  = nsample_list_var[[i]],
    mods        = var_models,
    prior_list  = var_uninf_prior_list,
    sd          = sd.sim_var,
    n_sim       = n_sim,
    alpha_crit_val = alpha,
    contr = contrasts_list[[i]])
   
  })
}
  

results_variability_Bay_nsample <- extract_success_rates_nsample(results_list_nsample_Bay_var, variability_scenario)


var_nsample_Bay <-  print_result_Bay_nsample(results_variability_Bay_nsample, c(variability_scenario, "average"), as.vector(nsample_list_var))
var_nsample_Bay$kable_result
BayesianMCPMod
Success probability results different sample sizes
variable linear exponential emax average
40, 20, 20, 40 0.73 0.74 0.74 0.7366667
30, 30, 30, 30 0.54 0.54 0.55 0.5433333
48, 24, 24, 48 0.83 0.83 0.82 0.8266667
36, 36, 36, 36 0.72 0.70 0.68 0.7000000

5 Comparison

In the following, the comparisons between the success probabilities (i.e.power values for frequentist set-up) of various scenarios and different parameters are visualized.

The following plots show the difference between the results from MCPModPack and BayesianMCPMod. The results of MCPModPack are shown as a line and the difference to the result with BayesianMCPMod is presented as a bar. The results for the different assumed true dose-response models (that were the basis for simulating the data) are shown in different colours.

5.1 varying expected effect for maximum dose

5.1.1 minimal scenario

setup
#table with results
kable(results_min_MCP)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different expected effects " = 5), font_size = 15, bold = TRUE)%>%
    add_header_above(c("MCPModPack " = 5), font_size = 15, bold = TRUE)
MCPModPack
Power results different expected effects
max_eff linear exponential emax Average
0.05 0.13 0.16 0.17 0.1533333
0.10 0.29 0.43 0.31 0.3433333
0.20 0.80 0.77 0.82 0.7966667
0.30 0.99 0.97 0.98 0.9800000
0.50 1.00 1.00 1.00 1.0000000
setup
minimal_Bay$kable_result
BayesianMCPMod
Power results different expected effects
variable linear exponential emax average
0.05 0.09 0.09 0.05 0.0766667
0.10 0.28 0.26 0.30 0.2800000
0.20 0.80 0.81 0.78 0.7966667
0.30 0.97 0.97 0.99 0.9766667
0.50 1.00 1.00 1.00 1.0000000
setup
data_plot_eff_min <- data.frame(
  max_eff = expectedEffect,
  max_eff_num = c(1, 2, 3, 4, 5),
  start_linear = results_min_MCP$linear,
  end_linear = minimal_Bay$result_table$linear,
  start_exp =  results_min_MCP$exponential,
  end_exp = minimal_Bay$result_table$exponential,
  start_emax = results_min_MCP$emax,
  end_emax = minimal_Bay$result_table$emax
)




# Create the plot
ggplot(data = data_plot_eff_min, aes(x = max_eff_num)) +
  geom_segment(aes(x =  max_eff_num - 0.3, xend = max_eff_num - 0.3, y = start_linear, yend = end_linear, color = "linear", size = "BayesianMCPMod")) +
  geom_segment(aes(x =  max_eff_num - 0.4, xend = max_eff_num - 0.2 , y = start_linear, yend = start_linear, color = "linear", size = "MCPModPack")) +
  geom_segment(aes(x =  max_eff_num - 0.05, xend = max_eff_num - 0.05 , y = start_exp, yend = end_exp, color = "exponential", size = "BayesianMCPMod")) +
  geom_segment(aes(x =  max_eff_num - 0.15, xend = max_eff_num + 0.05, y = start_exp, yend = start_exp, color = "exponential", size = "MCPModPack")) +
  geom_segment(aes(x =  max_eff_num + 0.2, xend = max_eff_num + 0.2, y = start_emax, yend = end_emax, color = "emax", size = "BayesianMCPMod")) +
  geom_segment(aes(x =  max_eff_num + 0.1, xend = max_eff_num  + 0.3, y = start_emax, yend = start_emax, color = "emax", size = "MCPModPack")) +
  scale_x_continuous(breaks = data_plot_eff_min$max_eff_num, labels = data_plot_eff_min$max_eff) +
  scale_color_manual(name = "Assumed true model", values = c("linear" = "red", "exponential" = "blue", "emax" = "darkgreen", "logistic" = "orange", "sigemax" = "purple", "beta" = "deepskyblue", "quadratic" = "deeppink"))+
  scale_size_manual( name = "Package", values = c("MCPModPack" = 0.5, "BayesianMCPMod" = 3))+
   theme_minimal() +
  ylab("Power") +
  theme_minimal() +
  ylab("Power") +
  xlab("expected effect for maximum dose")+
  ggtitle("Power values different expected effect for maximum dose") +
  geom_vline(xintercept = data_plot_eff_min$max_eff_num + 0.5, linetype="dashed", color = "black") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow= 2 , byrow= TRUE ), size = guide_legend(nrow= 2 , byrow= TRUE ))

5.1.2 monotonic scenario

setup
#table with results
kable(results_monotonic_MCP)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different expected effects " = 7), font_size = 15, bold = TRUE)%>%
    add_header_above(c("MCPModPack " = 7), font_size = 15, bold = TRUE)
MCPModPack
Power results different expected effects
max_eff linear exp emax logistic sigEmax Average
0.05 0.23 0.08 0.11 0.22 0.13 0.154
0.10 0.35 0.28 0.32 0.34 0.31 0.320
0.20 0.74 0.78 0.76 0.92 0.92 0.824
0.30 0.98 0.98 0.97 1.00 1.00 0.986
0.50 1.00 1.00 1.00 1.00 1.00 1.000
setup
monotonic_Bay$kable_result
BayesianMCPMod
Power results different expected effects
variable linear exponential emax logistic sigEmax average
0.05 0.15 0.15 0.15 0.16 0.15 0.152
0.10 0.30 0.31 0.32 0.37 0.35 0.330
0.20 0.78 0.79 0.79 0.91 0.87 0.828
0.30 0.99 0.97 1.00 1.00 1.00 0.992
0.50 1.00 1.00 1.00 1.00 1.00 1.000
setup
data_plot_eff_monotonic <- data.frame(
  max_eff = expectedEffect,
  max_eff_num = c(1, 2, 3, 4, 5),
  start_linear =  results_monotonic_MCP$linear,
  end_linear =  monotonic_Bay$result_table$linear,
  start_exp = results_monotonic_MCP$exp,
  end_exp =  monotonic_Bay$result_table$exponential,
  start_emax = results_monotonic_MCP$emax,
  end_emax =   monotonic_Bay$result_table$emax,
  start_logistic = results_monotonic_MCP$logistic,
  end_logistic =  monotonic_Bay$result_table$logistic,
  start_sigemax = results_monotonic_MCP$sigEmax,
  end_sigemax =  monotonic_Bay$result_table$sigEmax
)


# Create the plot
ggplot(data = data_plot_eff_monotonic, aes(x = max_eff_num)) +
  geom_segment(aes(x = max_eff_num - 0.4, xend = max_eff_num - 0.4, y = start_linear, yend = end_linear, color = "linear", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num - 0.45, xend = max_eff_num - 0.35 , y = start_linear, yend = start_linear, color = "linear", size = "MCPModPack")) +
  geom_segment(aes(x = max_eff_num - 0.2, xend = max_eff_num - 0.2 , y = start_exp, yend = end_exp, color = "exponential", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num- 0.25, xend = max_eff_num - 0.15, y = start_exp, yend = start_exp, color = "exponential", size = "MCPModPack")) +
  geom_segment(aes(x = max_eff_num, xend = max_eff_num, y = start_emax, yend = end_emax, color = "emax", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num - 0.05, xend = max_eff_num + 0.05, y = start_emax, yend = start_emax, color = "emax", size = "MCPModPack")) +
  geom_segment(aes(x = max_eff_num + 0.2, xend = max_eff_num + 0.2, y = start_logistic, yend = end_logistic, color = "logistic", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num + 0.15, xend = max_eff_num + 0.25, y = start_logistic, yend = start_logistic, color = "logistic", size = "MCPModPack")) +
  geom_segment(aes(x = max_eff_num + 0.4, xend = max_eff_num + 0.4, y = start_sigemax, yend = end_sigemax, color = "sigemax", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num + 0.35, xend = max_eff_num + 0.45, y = start_sigemax, yend = start_sigemax, color = "sigemax", size = "MCPModPack")) +
  scale_x_continuous(breaks = data_plot_eff_monotonic$max_eff_num, labels = data_plot_eff_monotonic$max_eff) +
  scale_color_manual(name = "Assumed true model", values = c("linear" = "red", "exponential" = "blue", "emax" = "darkgreen", "logistic" = "orange", "sigemax" = "purple", "beta" = "deepskyblue", "quadratic" = "deeppink"))+
  scale_size_manual( name = "Package", values = c("MCPModPack" = 0.5, "BayesianMCPMod" = 2))+
  theme_minimal() +
  ylab("Power") +
  xlab("expected effect for maximum dose")+
  ggtitle("Power values different expected effect for maximum dose")+
  geom_vline(xintercept = data_plot_eff_monotonic$max_eff_num + 0.5, linetype="dashed", color = "black") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow= 2 , byrow= TRUE ), size = guide_legend(nrow= 2 , byrow= TRUE ))

5.1.3 non - monotonic scenario

setup
# table results
kable(results_non_monotonic_MCP)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different expected effects" = 7), font_size = 15, bold = TRUE)%>%
    add_header_above(c("MCPModPack" = 7), font_size = 15, bold = TRUE)
MCPModPack
Power results different expected effects
Max_eff linear emax sigEmax quadratic beta average
0.05 0.1361826 0.1438824 0.1720734 0.1306460 0.1272818 0.1420132
0.10 0.2985730 0.3174089 0.4040673 0.2768840 0.2649841 0.3123835
0.20 0.7408952 0.7645582 0.8819700 0.6924798 0.6605961 0.7480998
0.30 0.9677230 0.9740876 0.9957499 0.9464217 0.9301074 0.9628180
0.50 0.9999913 0.9999952 1.0000000 0.9999419 0.9998595 0.9999576
setup
non_monotonic_Bay$kable_result
BayesianMCPMod
Power results different expected effects
variable linear emax sigEmax quadratic betaMod average
0.05 0.10 0.09 0.12 0.09 0.09 0.098
0.10 0.27 0.32 0.41 0.28 0.27 0.310
0.20 0.75 0.79 0.87 0.70 0.66 0.754
0.30 0.96 0.95 1.00 0.91 0.90 0.944
0.50 1.00 1.00 1.00 1.00 1.00 1.000
setup
data_plot_eff_non_monotonic <- data.frame(
  max_eff = expectedEffect,
  max_eff_num = c(1, 2, 3, 4, 5),
  start_linear =  results_non_monotonic_MCP$linear,
  end_linear =  non_monotonic_Bay$result_table$linear,
  start_emax = results_non_monotonic_MCP$emax,
  end_emax = non_monotonic_Bay$result_table$emax,
  start_sigemax =  results_non_monotonic_MCP$sigEmax,
  end_sigemax =non_monotonic_Bay$result_table$sigEmax,
  start_quadratic = results_non_monotonic_MCP$quadratic,
  end_quadratic = non_monotonic_Bay$result_table$quadratic,
  start_beta = results_non_monotonic_MCP$beta,
  end_beta = non_monotonic_Bay$result_table$betaMod
)


# Create the plot
ggplot(data = data_plot_eff_non_monotonic, aes(x = max_eff_num)) +
  geom_segment(aes(x = max_eff_num - 0.4, xend = max_eff_num - 0.4, y = start_linear, yend = end_linear, color = "linear", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num - 0.45, xend = max_eff_num - 0.35 , y = start_linear, yend = start_linear, color = "linear", size = "DoseFinding")) +
  geom_segment(aes(x = max_eff_num - 0.2, xend = max_eff_num - 0.2 , y = start_emax, yend = end_emax, color = "emax", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num - 0.25, xend = max_eff_num - 0.15, y = start_emax, yend = start_emax, color = "emax", size = "DoseFinding")) +
  geom_segment(aes(x = max_eff_num, xend = max_eff_num, y = start_sigemax, yend = end_sigemax, color = "sigemax", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num - 0.05, xend = max_eff_num + 0.05, y = start_sigemax, yend = start_sigemax, color = "sigemax", size = "DoseFinding")) +
  geom_segment(aes(x = max_eff_num + 0.2, xend = max_eff_num + 0.2, y = start_quadratic, yend = end_quadratic, color = "quadratic", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num + 0.15, xend = max_eff_num  + 0.25, y = start_quadratic, yend = start_quadratic, color = "quadratic", size = "DoseFinding")) +
  geom_segment(aes(x = max_eff_num + 0.4, xend = max_eff_num + 0.4, y = start_beta, yend = end_beta, color = "beta", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num + 0.35, xend = max_eff_num  + 0.45, y = start_beta, yend = start_beta, color = "beta", size = "DoseFinding")) +
  scale_x_continuous(breaks = data_plot_eff_non_monotonic$max_eff_num, labels = data_plot_eff_non_monotonic$max_eff) +
  scale_color_manual(name = "Assumed true model", values = c("linear" = "red", "exponential" = "blue", "emax" = "darkgreen", "logistic" = "orange", "sigemax" = "purple", "beta" = "deepskyblue", "quadratic" = "deeppink"))+
   scale_size_manual( name = "Package", values = c("DoseFinding" = 0.5, "BayesianMCPMod" = 2))+
  theme_minimal() +
  ylab("Power") +
  xlab("expected effect for maximum dose")+
  ggtitle("Power values for different expected effect for maximum dose")+
  geom_vline(xintercept = data_plot_eff_non_monotonic$max_eff_num + 0.5, linetype="dashed", color = "black") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow= 2 , byrow= TRUE ), size = guide_legend(nrow= 2 , byrow= TRUE ))

5.1.4 variability scenario

setup
#table with results
kable(results_var_MCP)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different expected effects" = 5), font_size = 15, bold = TRUE)%>%
    add_header_above(c("MCPModPack" = 5), font_size = 15, bold = TRUE)
MCPModPack
Power results different expected effects
Max_eff linear exponential emax average
0.05 0.12 0.12 0.12 0.1200000
0.10 0.19 0.28 0.20 0.2233333
0.20 0.60 0.65 0.52 0.5900000
0.30 0.82 0.87 0.84 0.8433333
0.50 1.00 1.00 1.00 1.0000000
setup
variability_Bay$kable_result
BayesianMCPMod
Power results different expected effects
variable linear exponential emax average
0.05 0.12 0.12 0.09 0.1100000
0.10 0.22 0.21 0.23 0.2200000
0.20 0.57 0.56 0.58 0.5700000
0.30 0.83 0.86 0.87 0.8533333
0.50 0.99 1.00 0.99 0.9933333
setup
data_plot_eff_var <- data.frame(
  max_eff = expectedEffect,
  max_eff_num = c(1, 2, 3, 4, 5),
  start_linear =  results_var_MCP$linear,
  end_linear = variability_Bay$result_table$linear,
  start_exp = results_var_MCP$exponential,
  end_exp = variability_Bay$result_table$exponential,
  start_emax = results_var_MCP$emax,
  end_emax = variability_Bay$result_table$emax
)


# Create the plot
ggplot(data = data_plot_eff_var, aes(x = max_eff_num)) +
  geom_segment(aes(x = max_eff_num - 0.3, xend = max_eff_num - 0.3, y = start_linear, yend = end_linear, color = "linear", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num - 0.4, xend = max_eff_num - 0.2 , y = start_linear, yend = start_linear, color = "linear", size = "MCPModPack")) +
  geom_segment(aes(x = max_eff_num - 0.05, xend = max_eff_num - 0.05, y = start_exp, yend = end_exp, color = "exponential", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num - 0.15, xend = max_eff_num + 0.05 , y = start_exp, yend = start_exp, color = "exponential", size = "MCPModPack")) +
  geom_segment(aes(x = max_eff_num + 0.2, xend = max_eff_num + 0.2, y = start_emax, yend = end_emax, color = "emax", size = "BayesianMCPMod")) +
  geom_segment(aes(x = max_eff_num + 0.1, xend = max_eff_num + 0.3 , y = start_emax, yend = start_emax, color = "emax", size = "MCPModPack")) +
  scale_x_continuous(breaks = data_plot_eff_var$max_eff_num, labels = data_plot_eff_var$max_eff) +
  scale_color_manual(name = "Assumed true model", values = c("linear" = "red", "exponential" = "blue", "emax" = "darkgreen", "logistic" = "orange", "sigemax" = "purple", "beta" = "deepskyblue", "quadratic" = "deeppink"))+
   scale_size_manual( name = "Package", values = c("MCPModPack" = 0.5, "BayesianMCPMod" = 3))+
  theme_minimal() +
  ylab("Power") +
  xlab("expected effect for maximum dose")+
  ggtitle("Power values for different expected effect for maximum dose")+
  geom_vline(xintercept = data_plot_eff_var$max_eff_num + 0.5, linetype="dashed", color = "black") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow= 2 , byrow= TRUE ), size = guide_legend(nrow= 2 , byrow= TRUE ))

5.2 varying sample size

setup
nsample_vector <- as.vector(c("(40, 20, 20, 20, 40)", "(30, 30, 30, 30, 30)", "(48, 24, 24, 24, 48)", "(36, 36, 36, 36, 36)"))

5.2.1 minimal scenario

setup
kable(results_min_MCP_nsample)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different sample sizes " = 5), font_size = 15, bold = TRUE)%>%
    add_header_above(c("MCPModPack" = 5), font_size = 15, bold = TRUE)
MCPModPack
Power results different sample sizes
N_sample Linear Exponential Emax Average
(40,20,20,20,40) 0.67 0.71 0.70 0.6933333
(30,30,30,30,30) 0.68 0.72 0.65 0.6833333
(48,24,24,24,48) 0.77 0.89 0.80 0.8200000
(36,36,36,36,36) 0.69 0.78 0.72 0.7300000
setup
minimal_nsample_Bay$kable_result
BayesianMCPMod
Success probability results different sample sizes
variable linear exponential emax average
(40,20,20,20,40) 0.67 0.73 0.67 0.6900000
(30,30,30,30,30) 0.66 0.64 0.66 0.6533333
(48,24,24,24,48) 0.81 0.81 0.81 0.8100000
(36,36,36,36,36) 0.67 0.70 0.73 0.7000000
setup
data_plot_nsample_min <- data.frame(
  sample_sizes = nsample_vector,
  sample_sizes_num = c(1, 2, 3, 4),
  start_linear = results_min_MCP_nsample$Linear,
  end_linear = minimal_nsample_Bay$result_table$linear,
  start_exp =  results_min_MCP_nsample$Exponential,
  end_exp = minimal_nsample_Bay$result_table$exponential,
  start_emax = results_min_MCP_nsample$Emax,
  end_emax = minimal_nsample_Bay$result_table$emax
)




# Create the plot
ggplot(data = data_plot_nsample_min, aes(x = sample_sizes_num)) +
  geom_segment(aes(x = sample_sizes_num - 0.3, xend = sample_sizes_num - 0.3, y = start_linear, yend = end_linear, color = "linear", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num - 0.4, xend = sample_sizes_num - 0.2 , y = start_linear, yend = start_linear, color = "linear", size = "MCPModPack")) +
  geom_segment(aes(x = sample_sizes_num - 0.05, xend = sample_sizes_num - 0.05 , y = start_exp, yend = end_exp, color = "exponential", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num - 0.15, xend = sample_sizes_num + 0.05, y = start_exp, yend = start_exp, color = "exponential", size = "MCPModPack")) +
  geom_segment(aes(x = sample_sizes_num + 0.2, xend = sample_sizes_num + 0.2, y = start_emax, yend = end_emax, color = "emax", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num + 0.1, xend = sample_sizes_num  + 0.3, y = start_emax, yend = start_emax, color = "emax", size = "MCPModPack")) +
  scale_x_continuous(breaks = data_plot_nsample_min$sample_sizes_num, labels = data_plot_nsample_min$sample_sizes) +
  scale_color_manual(name = "Assumed true model", values = c("linear" = "red", "exponential" = "blue", "emax" = "darkgreen", "logistic" = "orange", "sigemax" = "purple", "beta" = "deepskyblue", "quadratic" = "deeppink"))+
  scale_size_manual( name = "Package", values = c("MCPModPack" = 0.5, "BayesianMCPMod" = 3))+
   theme_minimal() +
  ylab("Power") +
  theme_minimal() +
  ylab("Power") +
  ylim(c(0.25,1))+
  xlab("sample sizes")+
  ggtitle("Power values different sample sizes") +
  geom_vline(xintercept = data_plot_nsample_min$sample_sizes_num + 0.5, linetype="dashed", color = "black") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow= 2 , byrow= TRUE ), size = guide_legend(nrow= 2 , byrow= TRUE ))

5.2.2 monotonic scenario

setup
kable(results_monotonic_MCP_nsample)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different sample sizes" = 7), font_size = 15, bold = TRUE)%>%
    add_header_above(c("MCPModPack" = 7), font_size = 15, bold = TRUE)
MCPModPack
Power results different sample sizes
N_sample Linear Exponential Emax Logistic sigEmax Average
(40,20,20,20,40) 0.73 0.74 0.80 0.81 0.76 0.768
(30,30,30,30,30) 0.57 0.66 0.57 0.84 0.71 0.670
(48,24,24,24,48) 0.76 0.79 0.83 0.91 0.89 0.836
(36,36,36,36,36) 0.58 0.71 0.74 0.90 0.80 0.746
setup
monotonic_nsample_Bay$kable_result
BayesianMCPMod
Success probability results different sample sizes
variable linear exponential emax logistic sigEmax average
(40,20,20,20,40) 0.80 0.80 0.77 0.88 0.84 0.818
(30,30,30,30,30) 0.65 0.70 0.63 0.79 0.75 0.704
(48,24,24,24,48) 0.84 0.85 0.78 0.87 0.86 0.840
(36,36,36,36,36) 0.77 0.76 0.74 0.91 0.85 0.806
setup
data_plot_nsample_monotonic <- data.frame(
  sample_sizes = nsample_vector,
  sample_sizes_num = c(1, 2, 3, 4),
  start_linear =  results_monotonic_MCP_nsample$Linear,
  end_linear =  monotonic_nsample_Bay$result_table$linear,
  start_exp = results_monotonic_MCP_nsample$Exponential,
  end_exp =  monotonic_nsample_Bay$result_table$exponential,
  start_emax = results_monotonic_MCP_nsample$Emax,
  end_emax =   monotonic_nsample_Bay$result_table$emax,
  start_logistic = results_monotonic_MCP_nsample$Logistic,
  end_logistic =  monotonic_nsample_Bay$result_table$logistic,
  start_sigemax = results_monotonic_MCP_nsample$sigEmax,
  end_sigemax =  monotonic_nsample_Bay$result_table$sigEmax
)


# Create the plot
ggplot(data = data_plot_nsample_monotonic, aes(x = sample_sizes_num)) +
  geom_segment(aes(x = sample_sizes_num - 0.4, xend = sample_sizes_num - 0.4, y = start_linear, yend = end_linear, color = "linear", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num - 0.45, xend = sample_sizes_num - 0.35 , y = start_linear, yend = start_linear, color = "linear", size = "MCPModPack")) +
  geom_segment(aes(x = sample_sizes_num - 0.2, xend = sample_sizes_num - 0.2 , y = start_exp, yend = end_exp, color = "exponential", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num - 0.25, xend = sample_sizes_num - 0.15, y = start_exp, yend = start_exp, color = "exponential", size = "MCPModPack")) +
  geom_segment(aes(x = sample_sizes_num, xend = sample_sizes_num, y = start_emax, yend = end_emax, color = "emax", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num - 0.05, xend = sample_sizes_num + 0.05, y = start_emax, yend = start_emax, color = "emax", size = "MCPModPack")) +
  geom_segment(aes(x = sample_sizes_num + 0.2, xend = sample_sizes_num + 0.2, y = start_logistic, yend = end_logistic, color = "logistic", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num + 0.15, xend = sample_sizes_num  + 0.25, y = start_logistic, yend = start_logistic, color = "logistic", size = "MCPModPack")) +
  geom_segment(aes(x = sample_sizes_num + 0.4, xend = sample_sizes_num + 0.4, y = start_sigemax, yend = end_sigemax, color = "sigemax", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num + 0.35, xend = sample_sizes_num  + 0.45, y = start_sigemax, yend = start_sigemax, color = "sigemax", size = "MCPModPack")) +
  scale_x_continuous(breaks = data_plot_nsample_monotonic$sample_sizes_num, labels = data_plot_nsample_monotonic$sample_sizes) +
  scale_color_manual(name = "Assumed true model", values = c("linear" = "red", "exponential" = "blue", "emax" = "darkgreen", "logistic" = "orange", "sigemax" = "purple", "beta" = "deepskyblue", "quadratic" = "deeppink"))+
  scale_size_manual( name = "Package", values = c("MCPModPack" = 0.5, "BayesianMCPMod" = 3))+
  theme_minimal() +
  ylab("Power") +
  ylim(c(0.25,1))+
  xlab("sample sizes")+
  ggtitle("Power values different sample sizes")+
  geom_vline(xintercept = data_plot_nsample_monotonic$sample_sizes_num + 0.5, linetype="dashed", color = "black") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow= 2 , byrow= TRUE ), size = guide_legend(nrow= 2 , byrow= TRUE ))

5.2.3 non-monotonic scenario

setup
kable(results_non_monotonic_MCP_nsample)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different sample sizes" = 6), font_size = 15, bold = TRUE)%>%
    add_header_above(c("MCPModPack" = 6), font_size = 15, bold = TRUE)
MCPModPack
Power results different sample sizes
N_sample Linear Emax sigEmax quadratic beta
(40,20,20,20,40) 0.7182897 0.7452513 0.8118555 0.5991157 0.5492169
(30,30,30,30,30) 0.6265819 0.6539679 0.7872977 0.5801699 0.5521925
(48,24,24,24,48) 0.7869752 0.8122048 0.8716953 0.6689755 0.6177505
(36,36,36,36,36) 0.6994504 0.7248996 0.8504362 0.6509720 0.6195153
setup
non_monotonic_nsample_Bay$kable_result
BayesianMCPMod
Success probability results different sample sizes
variable linear emax sigEmax quadratic betaMod average
(40,20,20,20,40) 0.73 0.76 0.81 0.58 0.50 0.676
(30,30,30,30,30) 0.63 0.64 0.77 0.52 0.49 0.610
(48,24,24,24,48) 0.74 0.75 0.83 0.65 0.59 0.712
(36,36,36,36,36) 0.67 0.67 0.78 0.59 0.56 0.654
setup
data_plot_nsample_non_monotonic <- data.frame(
  sample_sizes = nsample_vector,
  sample_sizes_num = c(1, 2, 3, 4),
  start_linear =  results_non_monotonic_MCP_nsample$Linear,
  end_linear =  non_monotonic_nsample_Bay$result_table$linear,
  start_emax = results_non_monotonic_MCP_nsample$Emax,
  end_emax = non_monotonic_nsample_Bay$result_table$emax,
  start_sigemax =  results_non_monotonic_MCP_nsample$sigEmax,
  end_sigemax =non_monotonic_nsample_Bay$result_table$sigEmax,
  start_quadratic = results_non_monotonic_MCP_nsample$quadratic,
  end_quadratic = non_monotonic_nsample_Bay$result_table$quadratic,
  start_beta = results_non_monotonic_MCP_nsample$beta,
  end_beta = non_monotonic_nsample_Bay$result_table$betaMod
)


# Create the plot
ggplot(data = data_plot_nsample_non_monotonic, aes(x = sample_sizes_num)) +
  geom_segment(aes(x = sample_sizes_num - 0.4, xend = sample_sizes_num - 0.4, y = start_linear, yend = end_linear, color = "linear", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num - 0.45, xend = sample_sizes_num - 0.35 , y = start_linear, yend = start_linear, color = "linear", size = "DoseFinding")) +
  geom_segment(aes(x = sample_sizes_num - 0.2, xend = sample_sizes_num - 0.2 , y = start_emax, yend = end_emax, color = "emax", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num - 0.25, xend = sample_sizes_num - 0.15, y = start_emax, yend = start_emax, color = "emax", size = "DoseFinding")) +
  geom_segment(aes(x = sample_sizes_num, xend = sample_sizes_num, y = start_sigemax, yend = end_sigemax, color = "sigemax", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num - 0.05, xend = sample_sizes_num + 0.05, y = start_sigemax, yend = start_sigemax, color = "sigemax", size = "DoseFinding")) +
  geom_segment(aes(x = sample_sizes_num + 0.1, xend = sample_sizes_num + 0.1, y = start_quadratic, yend = end_quadratic, color = "quadratic", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num + 0.05, xend = sample_sizes_num  + 0.15, y = start_quadratic, yend = start_quadratic, color = "quadratic", size = "DoseFinding")) +
  geom_segment(aes(x = sample_sizes_num + 0.3, xend = sample_sizes_num + 0.3, y = start_beta, yend = end_beta, color = "beta", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num + 0.25, xend = sample_sizes_num  + 0.35, y = start_beta, yend = start_beta, color = "beta", size = "DoseFinding")) +
  scale_x_continuous(breaks = data_plot_nsample_non_monotonic$sample_sizes_num, labels = data_plot_nsample_non_monotonic$sample_sizes) +
  scale_color_manual(name = "Assumed true model", values = c("linear" = "red", "exponential" = "blue", "emax" = "darkgreen", "logistic" = "orange", "sigemax" = "purple", "beta" = "deepskyblue", "quadratic" = "deeppink"))+
   scale_size_manual( name = "Package", values = c("DoseFinding" = 0.5, "BayesianMCPMod" = 3))+
  theme_minimal() +
  ylab("Power") +
  ylim(c(0.25,1))+
  xlab("sample sizes")+
  ggtitle("Power values different sample sizes")+
  geom_vline(xintercept = data_plot_nsample_non_monotonic$sample_sizes_num + 0.5, linetype="dashed", color = "black") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow= 2 , byrow= TRUE ), size = guide_legend(nrow= 2 , byrow= TRUE ))

5.2.4 variability scenario

setup
kable(results_var_MCP_nsample)%>%
  kable_classic(full_width = TRUE)%>%
    add_header_above(c("Power results different sample sizes" = 5), font_size = 15, bold = TRUE)%>%
    add_header_above(c("Variability scenario" = 5), font_size = 15, bold = TRUE)
Variability scenario
Power results different sample sizes
N_sample Linear Exponential Emax Average
(40,20,20,20,40) 0.83 0.71 0.69 0.7433333
(30,30,30,30,30) 0.64 0.62 0.68 0.6466667
(48,24,24,24,48) 0.81 0.76 0.85 0.8066667
(36,36,36,36,36) 0.76 0.60 0.62 0.6600000
setup
var_nsample_Bay$kable_result
BayesianMCPMod
Success probability results different sample sizes
variable linear exponential emax average
40, 20, 20, 40 0.73 0.74 0.74 0.7366667
30, 30, 30, 30 0.54 0.54 0.55 0.5433333
48, 24, 24, 48 0.83 0.83 0.82 0.8266667
36, 36, 36, 36 0.72 0.70 0.68 0.7000000
setup
data_plot_nsample_var <- data.frame(
  sample_sizes = nsample_vector,
  sample_sizes_num = c(1, 2, 3, 4),
  start_linear =  results_var_MCP_nsample$Linear,
  end_linear = var_nsample_Bay$result_table$linear,
  start_exp = results_var_MCP_nsample$Exponential,
  end_exp = var_nsample_Bay$result_table$exponential,
  start_emax = results_var_MCP_nsample$Emax,
  end_emax = var_nsample_Bay$result_table$emax
)


# Create the plot
ggplot(data = data_plot_nsample_var, aes(x = sample_sizes_num)) +
  geom_segment(aes(x = sample_sizes_num - 0.3, xend = sample_sizes_num - 0.3, y = start_linear, yend = end_linear, color = "linear", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num - 0.4, xend = sample_sizes_num - 0.2 , y = start_linear, yend = start_linear, color = "linear", size = "MCPModPack")) +
  geom_segment(aes(x = sample_sizes_num - 0.05, xend = sample_sizes_num - 0.05, y = start_exp, yend = end_exp, color = "exponential", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num - 0.15, xend = sample_sizes_num + 0.05 , y = start_exp, yend = start_exp, color = "exponential", size = "MCPModPack")) +
  geom_segment(aes(x = sample_sizes_num + 0.2, xend = sample_sizes_num + 0.2, y = start_emax, yend = end_emax, color = "emax", size = "BayesianMCPMod")) +
  geom_segment(aes(x = sample_sizes_num + 0.1, xend = sample_sizes_num + 0.3 , y = start_emax, yend = start_emax, color = "emax", size = "MCPModPack")) +
  scale_x_continuous(breaks = data_plot_nsample_var$sample_sizes_num, labels = data_plot_nsample_var$sample_sizes) +
  scale_color_manual(name = "Assumed true model", values = c("linear" = "red", "exponential" = "blue", "emax" = "darkgreen", "logistic" = "orange", "sigemax" = "purple", "beta" = "deepskyblue", "quadratic" = "deeppink"))+
   scale_size_manual( name = "Package", values = c("MCPModPack" = 0.5, "BayesianMCPMod" = 3))+
  theme_minimal() +
  ylab("Power") +
  ylim(c(0.25,1))+
  xlab("sample sizes")+
  ggtitle("Power values different sample sizes")+
  geom_vline(xintercept = data_plot_nsample_var$sample_sizes_num + 0.5, linetype="dashed", color = "black") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow= 2 , byrow= TRUE ), size = guide_legend(nrow= 2 , byrow= TRUE ))

As expected for all the different scenarios operating characteristics for BayesianMCPMod with non-informative prior match with operating characteristics for frequentist MCPMod.

5.3 convergence of power values

In the following simulations, we examine the convergence of power values for an increasing number of simulations. We are considering the following number of simulations: 100, 500, 1000, 2500, 5000, 10000.

5.3.1 minimal BayesianMCPMod

setup
# Define a vector of different n_sim values
n_sim_values_list <- as.list(c(10, 50, 100, 250, 500, 1000)) #as.list(c(100, 500, 1000, 2500, 5000, 10000))
n_sim_values <- c(10, 50, 100, 250, 500, 1000) #c(100, 500, 1000, 2500, 5000, 10000)

# Initialize a list to store the results
results_list_nsim_Bay <- list()

min_models <- Mods(linear=NULL,
                     exponential = exp.g,
                     emax = emax.g,
                     doses = doses.sim,
                     placEff = plc.guess,
                     maxEff = expectedEffect_fix,
                     direction = "increasing")
#optimal contrasts
contM <- getContr(mods = min_models,
                    dose_levels = doses.sim,
                    prior_list = uninf_prior_list,
                    dose_weights = c(1,1,1,1,1))

chunks <- chunkVector(seq_along(n_sim_values_list), getDoParWorkers())

  
results_list_nsim_Bay <- foreach(k = chunks, .combine = c, .export = c(as.character(min_models), as.character(contM))) %dorng% {

  lapply(k, function (i) {


 # Simulation step
  success_probabilities_min <- assessDesign(
    n_patients  = Nsample,
    mods        = min_models,
    prior_list  = uninf_prior_list,
    sd          = sd.sim,
    n_sim       = n_sim_values_list[[i]],
    alpha_crit_val = alpha,
    contr = contM)

  })

}

results_nsim_Bay <- extract_success_rates_nsim(results_list_nsim_Bay, min_scenario, n_sim_values)

5.3.2 minimal MCPModPack

setup
# assumed dose - response model
sim_models_min_linear = list(linear = NA,
                  max_effect = expectedEffect_fix,
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)

sim_models_min_exp = list(exponential = 4.447149,
                  max_effect = expectedEffect_fix,
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)

sim_models_min_emax = list( emax = 0.6666667,
                  max_effect = expectedEffect_fix,
                  sd = rep(sd.sim, length(doses.sim)),
                  placebo_effect = plc.guess)




list_models <- list(sim_models_min_linear, sim_models_min_linear, sim_models_min_linear, sim_models_min_linear, sim_models_min_linear, sim_models_min_linear,
                    sim_models_min_exp, sim_models_min_exp, sim_models_min_exp, sim_models_min_exp, sim_models_min_exp, sim_models_min_exp, 
                    sim_models_min_emax, sim_models_min_emax, sim_models_min_emax, sim_models_min_emax, sim_models_min_emax, sim_models_min_emax)
n_sim_values_list <- as.list(c(10, 50, 100, 250, 500, 1000, 
                               10, 50, 100, 250, 500, 1000, 
                               10, 50, 100, 250, 500, 1000))
#as.list(c(100, 500, 1000, 2500, 5000, 10000, 
    #                           100, 500, 1000, 2500, 5000, 10000, 
                 #              100, 500, 1000, 2500, 5000, 10000))

chunks <- chunkVector(seq_along(list_models), getDoParWorkers())

results_list_nsim_MCP <- foreach(k = chunks, .combine = c, .export = c(as.character(min_modelsPack)))  %dorng% {
  
  lapply(k, function (i) {
    # Simulation parameters
  sim_parameters = list(n = Nsample,
                      doses = doses.sim,
                      dropout_rate = 0.0,
                      go_threshold = 0.1,
                      nsims = n_sim_values_list[[i]])
    
    func_sim(min_modelsPack, list_models[[i]], sim_parameters) 
    
  })
  
}

5.3.3 monotonic BayesianMCPMod

setup
# Initialize a list to store the results
 results_list_nsim_Bay_monotonic <- list()
 
monotonic_models <- Mods(linear=NULL,
                         exponential = exp.g,
                         emax = emax.g,
                         logistic=logit.g,
                         sigEmax = sigEmax.g,
                         doses = doses.sim,
                         placEff = plc.guess,
                         maxEff = expectedEffect_fix,
                         direction = "increasing")

#optimal contrasts
contM <- getContr(mods = monotonic_models,
                  dose_levels = doses.sim,
                  prior_list = uninf_prior_list,
                  dose_weights = c(1,1,1,1,1))

n_sim_values_list <- as.list(c(10, 50, 100, 250, 500, 1000)) #as.list(c(100, 500, 1000, 2500, 5000, 10000))

chunks <- chunkVector(seq_along(n_sim_values_list), getDoParWorkers())


results_list_nsim_Bay_monotonic <- foreach(k = chunks, .combine = c, .export = c(as.character(monotonic_models), as.character(contM))) %dorng% {

  lapply(k, function (i) {


 # Simulation step
  success_probabilities_monotonic <- assessDesign(
    n_patients  = Nsample,
    mods        = monotonic_models,
    prior_list  = uninf_prior_list,
    sd          = sd.sim,
    n_sim       = n_sim_values_list[[i]],
    alpha_crit_val = alpha,
    contr = contM)

  })

}



results_nsim_Bay_monotonic <- extract_success_rates_nsim(results_list_nsim_Bay_monotonic, monotonic_scenario, n_sim_values)

5.3.4 monotonic MCPModPack

setup
sim_models_monotonic_linear$max_effect <- expectedEffect_fix
sim_models_monotonic_exp$max_effect <- expectedEffect_fix
sim_models_monotonic_emax$max_effect <- expectedEffect_fix
sim_models_monotonic_logistic$max_effect <- expectedEffect_fix
sim_models_monotonic_sigemax$max_effect <- expectedEffect_fix

list_models <- list(sim_models_monotonic_linear, sim_models_monotonic_linear, sim_models_monotonic_linear, sim_models_monotonic_linear, sim_models_monotonic_linear, sim_models_monotonic_linear, 
                    sim_models_monotonic_exp, sim_models_monotonic_exp, sim_models_monotonic_exp, sim_models_monotonic_exp, sim_models_monotonic_exp, sim_models_monotonic_exp, 
                    sim_models_monotonic_emax, sim_models_monotonic_emax, sim_models_monotonic_emax, sim_models_monotonic_emax, sim_models_monotonic_emax, sim_models_monotonic_emax,
                    sim_models_monotonic_logistic, sim_models_monotonic_logistic, sim_models_monotonic_logistic, sim_models_monotonic_logistic, sim_models_monotonic_logistic, sim_models_monotonic_logistic,
                    sim_models_monotonic_sigemax, sim_models_monotonic_sigemax, sim_models_monotonic_sigemax, sim_models_monotonic_sigemax, sim_models_monotonic_sigemax, sim_models_monotonic_sigemax
                    )
n_sim_values_list <- as.list(c(10, 50, 100, 250, 500, 1000,
                               10, 50, 100, 250, 500, 1000,
                               10, 50, 100, 250, 500, 1000,
                               10, 50, 100, 250, 500, 1000,
                               10, 50, 100, 250, 500, 1000))
  
  #as.list(c(100, 500, 1000, 2500, 5000, 10000,
              #                 100, 500, 1000, 2500, 5000, 10000,
                 #              100, 500, 1000, 2500, 5000, 10000,
                  #             100, 500, 1000, 2500, 5000, 10000,
                    #           100, 500, 1000, 2500, 5000, 10000))

chunks <- chunkVector(seq_along(list_models), getDoParWorkers())

results_list_nsim_MCP_monotonic <- foreach(k = chunks, .combine = c, .export = c(as.character(monotonic_modelsPack)))  %dorng% {
  
  lapply(k, function (i) {
    # Simulation parameters
  sim_parameters = list(n = Nsample,
                      doses = doses.sim,
                      dropout_rate = 0.0,
                      go_threshold = 0.1,
                      nsims = n_sim_values_list[[i]])
    
    func_sim(monotonic_modelsPack, list_models[[i]], sim_parameters) 
    
  })
  
}

5.3.5 non-monotonic BayesianMCPMod

setup
# Initialize list to store results
results_list_nsim_Bay_non_monotonic <- list()

non_monotonic_models <- Mods(linear=NULL,
                             emax = emax.g,
                             sigEmax = sigEmax.g,
                             quadratic = quad.g,
                             betaMod = beta.g,
                             doses = doses.sim,
                             placEff = plc.guess,
                             maxEff = expectedEffect_fix,
                             direction = "increasing",
                             addArgs = list(scal=9.6))

  #optimal contrasts
contM <- getContr(mods = non_monotonic_models,
                  dose_levels = doses.sim,
                  prior_list = uninf_prior_list,
                  dose_weights = c(1,1,1,1,1))

n_sim_values_list <- as.list(c(10, 50, 100, 250, 500, 1000))#as.list(c(100, 500, 1000, 2500, 5000, 10000))

chunks <- chunkVector(seq_along(n_sim_values_list), getDoParWorkers())


results_list_nsim_Bay_non_monotonic <- foreach(k = chunks, .combine = c, .export = c(as.character(non_monotonic_models), as.character(contM))) %dorng% {

  lapply(k, function (i) {


 # Simulation step
  success_probabilities_non_monotonic <- assessDesign(
    n_patients  = Nsample,
    mods        = non_monotonic_models,
    prior_list  = uninf_prior_list,
    sd          = sd.sim,
    n_sim       = n_sim_values_list[[i]],
    alpha_crit_val = alpha,
    contr = contM)

  })

}



results_nsim_Bay_non_monotonic <- extract_success_rates_nsim(results_list_nsim_Bay_non_monotonic, non_monotonic_scenario, n_sim_values)

5.3.6 non-monotonic MCPModPack

setup
# linear with DoseFinding package
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(linear = NULL)
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(200))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
results_nsim_non_monotonic_linear <- result$linear
setup
# emax with DoseFinding package
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(emax = 0.6666667)
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(200))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
results_nsim_non_monotonic_emax <- result$`emax (ED50=0.6666667)`
setup
# sigemax with DoseFinding package
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(sigEmax = rbind(c(1.528629, 4.087463)))
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(200))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
results_nsim_non_monotonic_sigemax <- result$`sigEmax (ED50=1.528629,h=4.087463)`
setup
#quadratic with DoseFinding package
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(quadratic = -0.09688711)
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(200))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)

results_nsim_non_monotonic_quadratic <- result$`quadratic (delta=-0.09688711)`
setup
# beta with DoseFinding package
mvt_control <- mvtnorm.control(maxpts = 30000, abseps = 0.001, releps = 0)
doses <- c(0,1,2,4,8)
allocation <- c(1,1,1,1,1)
alpha <- 0.05
addArgs <- list(off = 0.08, scal = 9.6)

resp_mod_list <- list(betaMod = rbind(c(1.396434, 1.040978)))
cand_mod_list <- list(linear = NULL, emax = 0.6666667, quadratic = -0.09688711, betaMod = rbind(c(1.396434, 1.040978)), sigEmax = rbind(c(1.528629, 4.087463)))
mods <- do.call(Mods, append(cand_mod_list, 
        list(maxEff = 1, doses = doses, placEff = 0, addArgs = addArgs)))
cont_mat <- optContr(mods, w = allocation)

grd <- expand.grid(max_eff = c(0.2), sd = c(0.4), n_total = c(200))
power_list <- vector("list", nrow(grd))
dat_n <- t(sapply(grd$n_total, function(n) round(n*allocation/sum(allocation))))
colnames(dat_n) <- paste0("D", doses)
for (i in 1:nrow(grd)) {
  ######
  mods <- do.call(Mods, append(resp_mod_list, list(maxEff = grd$max_eff[i], doses = doses,
  placEff = 0, addArgs = addArgs)))
  n <- dat_n[i, ]
  power_list[[i]] <- powMCT(cont_mat, alpha, mods, n, grd$sd[i], control = mvt_control)
}
power_df <- as.data.frame(do.call("rbind", power_list))
mu <- getResp(mods)
parList <- attr(mu, 'parList')
mod_nams <- DoseFinding:::getModNams(parList)
colnames(power_df) <- mod_nams
power_df[['mean power']] <- apply(power_df, 1, mean)
result <- cbind(grd, power_df, dat_n)
results_nsim_non_monotonic_beta <- result$`betaMod (delta1=1.396434,delta2=1.040978,scal=9.6)`
setup
results_nsim_MCP_non_monotonic <- data.table(
  linear = results_nsim_non_monotonic_linear,
  emax = results_nsim_non_monotonic_emax,
  sigemax = results_nsim_non_monotonic_sigemax,
  quadratic = results_nsim_non_monotonic_quadratic,
  beta = results_nsim_non_monotonic_beta
)

5.3.7 variability BayesianMCPMod

setup
# Initialize list to store results
results_list_nsim_Bay_var <- list()

var_models <- Mods(linear=NULL,
                   exponential = exp.g,
                   emax = emax.g,
                   doses = doses_var,
                   placEff = plc.guess,
                   maxEff = expectedEffect_fix,
                   direction = "increasing")

contM <- getContr(mods = var_models,
                  dose_levels = doses_var,
                  prior_list = var_uninf_prior_list,
                  dose_weights = c(1,1,1,1))

n_sim_values_list <- as.list(c(10, 50, 100, 250, 500, 1000))#as.list(c(100, 500, 1000, 2500, 5000, 10000))

chunks <- chunkVector(seq_along(n_sim_values_list), getDoParWorkers())


results_list_nsim_Bay_var <- foreach(k = chunks, .combine = c, .export = c(as.character(var_models), as.character(contM))) %dorng% {

  lapply(k, function (i) {


 # Simulation step
  success_probabilities_var <- assessDesign(
    n_patients  = Nsample_var,
    mods        = var_models,
    prior_list  = var_uninf_prior_list,
    sd          = sd.sim_var,
    n_sim       = n_sim_values_list[[i]],
    alpha_crit_val = alpha,
    contr = contM)

  })

}



results_nsim_Bay_var <- extract_success_rates_nsim(results_list_nsim_Bay_var, variability_scenario, n_sim_values)

5.3.8 variability MCPModPack

setup
results_list_nsim_MCP_var = list(
  power_var_nsim_linear = vector(),
  power_var_nsim_exp = vector(),
  power_var_nsim_emax= vector())

# assume maximum effect = 0.2
sim_models_var_linear$max_effect <- expectedEffect_fix
sim_models_var_exp$max_effect <- expectedEffect_fix
sim_models_var_emax$max_effect <- expectedEffect_fix

sim_parameters_var$n <- Nsample_var


list_models <- list(sim_models_var_linear, sim_models_var_linear, sim_models_var_linear, sim_models_var_linear, sim_models_var_linear, sim_models_var_linear,
                    sim_models_var_exp, sim_models_var_exp, sim_models_var_exp, sim_models_var_exp, sim_models_var_exp, sim_models_var_exp,
                    sim_models_var_emax, sim_models_var_emax, sim_models_var_emax, sim_models_var_emax, sim_models_var_emax, sim_models_var_emax)
n_sim_values_list <- as.list(c(10, 50, 100, 250, 500, 1000,
                               10, 50, 100, 250, 500, 1000,
                               10, 50, 100, 250, 500, 1000))

#  as.list(c(100, 500, 1000, 2500, 5000, 10000,
     #                          100, 500, 1000, 2500, 5000, 10000,
           #                    100, 500, 1000, 2500, 5000, 10000))

chunks <- chunkVector(seq_along(list_models), getDoParWorkers())

results_list_nsim_MCP_var <- foreach(k = chunks, .combine = c, .export = c(as.character(var_modelsPack)))  %dorng% {
  
  lapply(k, function (i) {
    # Simulation parameters
  sim_parameters = list(n = Nsample_var,
                      doses = doses_var,
                      dropout_rate = 0.0,
                      go_threshold = 0.1,
                      nsims = n_sim_values_list[[i]])
    
    func_sim(var_modelsPack, list_models[[i]], sim_parameters) 
    
  })
  
}

5.4 Results

The following plots show the result of the convergence test of the power values for an increasing number of simulations. The difference between the success probabilities of the frequency and Bayesian simulations is shown. In the non-monotonous scenario, the DoseFinding package was used instead of the MCPModPack. This package does not simulate but calculates the success probabilities (i.e. no specific number of iterations needs to be defined). The difference between the Bayesian results and the value of the DoseFinding calculation is also shown.

5.4.1 minimal scenario

setup
#safe results in data.table for plot
results_nsim <- data.table(
  MCP_linear = c(results_list_nsim_MCP[[1]]$sim_results$power, results_list_nsim_MCP[[2]]$sim_results$power, results_list_nsim_MCP[[3]]$sim_results$power, results_list_nsim_MCP[[4]]$sim_results$power, results_list_nsim_MCP[[5]]$sim_results$power,  results_list_nsim_MCP[[6]]$sim_results$power),
  MCP_exp = c(results_list_nsim_MCP[[7]]$sim_results$power, results_list_nsim_MCP[[8]]$sim_results$power, results_list_nsim_MCP[[9]]$sim_results$power, results_list_nsim_MCP[[10]]$sim_results$power, results_list_nsim_MCP[[11]]$sim_results$power, results_list_nsim_MCP[[12]]$sim_results$power),
  MCP_emax = c(results_list_nsim_MCP[[13]]$sim_results$power, results_list_nsim_MCP[[14]]$sim_results$power, results_list_nsim_MCP[[15]]$sim_results$power, results_list_nsim_MCP[[16]]$sim_results$power, results_list_nsim_MCP[[17]]$sim_results$power,  results_list_nsim_MCP[[18]]$sim_results$power),
  Bay_linear = results_nsim_Bay$Bay_linear,
  Bay_exp = results_nsim_Bay$Bay_exponential,
  Bay_emax = results_nsim_Bay$Bay_emax)

results_nsim_diff <- data.table(
  linear = results_nsim$Bay_linear - results_nsim$MCP_linear,
  exponential = results_nsim$Bay_exp - results_nsim$MCP_exp,
  emax = results_nsim$Bay_emax - results_nsim$MCP_emax,
  n_sim = n_sim_values
)

results_nsim_diff <- melt(results_nsim_diff, id.vars = "n_sim")

plot_nsim <- plot_power_deviation(results_nsim_diff, results_nsim_diff$n_sim, "nsim")
plot_nsim

5.4.2 monotonic scenario

setup
#safe results in data.table for plot
results_nsim_monotonic <- data.table(
  MCP_linear = c(results_list_nsim_MCP_monotonic[[1]]$sim_results$power, results_list_nsim_MCP_monotonic[[2]]$sim_results$power, results_list_nsim_MCP_monotonic[[3]]$sim_results$power, results_list_nsim_MCP_monotonic[[4]]$sim_results$power, results_list_nsim_MCP_monotonic[[5]]$sim_results$power, results_list_nsim_MCP_monotonic[[6]]$sim_results$power),
  MCP_exp = c(results_list_nsim_MCP_monotonic[[7]]$sim_results$power, results_list_nsim_MCP_monotonic[[8]]$sim_results$power, results_list_nsim_MCP_monotonic[[9]]$sim_results$power, results_list_nsim_MCP_monotonic[[10]]$sim_results$power, results_list_nsim_MCP_monotonic[[11]]$sim_results$power, results_list_nsim_MCP_monotonic[[12]]$sim_results$power),
  MCP_emax = c(results_list_nsim_MCP_monotonic[[13]]$sim_results$power, results_list_nsim_MCP_monotonic[[14]]$sim_results$power, results_list_nsim_MCP_monotonic[[15]]$sim_results$power, results_list_nsim_MCP_monotonic[[16]]$sim_results$power, results_list_nsim_MCP_monotonic[[17]]$sim_results$power, results_list_nsim_MCP_monotonic[[18]]$sim_results$power),
  MCP_logistic = c(results_list_nsim_MCP_monotonic[[19]]$sim_results$power, results_list_nsim_MCP_monotonic[[20]]$sim_results$power, results_list_nsim_MCP_monotonic[[21]]$sim_results$power, results_list_nsim_MCP_monotonic[[22]]$sim_results$power, results_list_nsim_MCP_monotonic[[23]]$sim_results$power, results_list_nsim_MCP_monotonic[[24]]$sim_results$power),
  MCP_sigemax = c(results_list_nsim_MCP_monotonic[[25]]$sim_results$power, results_list_nsim_MCP_monotonic[[26]]$sim_results$power, results_list_nsim_MCP_monotonic[[27]]$sim_results$power, results_list_nsim_MCP_monotonic[[28]]$sim_results$power, results_list_nsim_MCP_monotonic[[29]]$sim_results$power, results_list_nsim_MCP_monotonic[[30]]$sim_results$power),
  Bay_linear = results_nsim_Bay_monotonic$Bay_linear,
  Bay_exp = results_nsim_Bay_monotonic$Bay_exponential,
  Bay_emax = results_nsim_Bay_monotonic$Bay_emax,
  Bay_logistic = results_nsim_Bay_monotonic$Bay_logistic,
  Bay_sigemax = results_nsim_Bay_monotonic$Bay_sigEmax)

results_nsim_diff_monotonic <- data.table(
  linear = results_nsim_monotonic$Bay_linear - results_nsim_monotonic$MCP_linear,
  exponential = results_nsim_monotonic$Bay_exp - results_nsim_monotonic$MCP_exp,
  emax = results_nsim_monotonic$Bay_emax - results_nsim_monotonic$MCP_emax,
  logistic = results_nsim_monotonic$Bay_logistic - results_nsim_monotonic$MCP_logistic,
  sigemax = results_nsim_monotonic$Bay_sigemax - results_nsim_monotonic$MCP_sigemax,
  n_sim = n_sim_values
)

results_nsim_diff_monotonic <- melt(results_nsim_diff_monotonic, id.vars = "n_sim")

plot_nsim_monotonic <-  plot_power_deviation(results_nsim_diff_monotonic, results_nsim_diff_monotonic$n_sim, "nsim")
plot_nsim_monotonic

5.4.3 non - monotonic scenario

setup
#safe results in data.table for plot
results_nsim_non_monotonic <- data.table(
  MCP_linear = replicate(6,results_nsim_MCP_non_monotonic$linear),
  MCP_emax = replicate(6,results_nsim_MCP_non_monotonic$emax),
  MCP_sigemax = replicate(6,results_nsim_MCP_non_monotonic$sigemax),
  MCP_quadratic = replicate(6,results_nsim_MCP_non_monotonic$quadratic),
  MCP_beta = replicate(6,results_nsim_MCP_non_monotonic$beta),
  Bay_linear = results_nsim_Bay_non_monotonic$Bay_linear,
  Bay_emax = results_nsim_Bay_monotonic$Bay_emax,
  Bay_sigemax = results_nsim_Bay_non_monotonic$Bay_sigEmax,
  Bay_quadratic = results_nsim_Bay_non_monotonic$Bay_quadratic,
  Bay_beta = results_nsim_Bay_non_monotonic$Bay_betaMod)

results_nsim_diff_non_monotonic <- data.table(
  linear = results_nsim_non_monotonic$Bay_linear - results_nsim_non_monotonic$MCP_linear,
  emax = results_nsim_non_monotonic$Bay_emax - results_nsim_non_monotonic$MCP_emax,
  sigemax =results_nsim_non_monotonic$Bay_sigemax - results_nsim_non_monotonic$MCP_sigemax,
  quadratic = results_nsim_non_monotonic$Bay_quadratic - results_nsim_non_monotonic$MCP_quadratic,
  beta = results_nsim_non_monotonic$Bay_beta -  results_nsim_non_monotonic$MCP_beta,
  n_sim = n_sim_values
)


results_nsim_diff_non_monotonic <- melt(results_nsim_diff_non_monotonic, id.vars = "n_sim")

plot_nsim_non_monotonic <-  plot_power_deviation(results_nsim_diff_non_monotonic, results_nsim_diff_non_monotonic$n_sim, "nsim")
plot_nsim_non_monotonic

5.4.4 variability scenario

setup
#safe results in data.table for plot
results_nsim_var <- data.table(
  MCP_linear = c(results_list_nsim_MCP_var[[1]]$sim_results$power, results_list_nsim_MCP_var[[2]]$sim_results$power, results_list_nsim_MCP_var[[3]]$sim_results$power, results_list_nsim_MCP_var[[4]]$sim_results$power, results_list_nsim_MCP_var[[5]]$sim_results$power,  results_list_nsim_MCP_var[[6]]$sim_results$power),
  MCP_exp = c(results_list_nsim_MCP_var[[7]]$sim_results$power, results_list_nsim_MCP_var[[8]]$sim_results$power, results_list_nsim_MCP_var[[9]]$sim_results$power, results_list_nsim_MCP_var[[10]]$sim_results$power, results_list_nsim_MCP_var[[11]]$sim_results$power,  results_list_nsim_MCP_var[[12]]$sim_results$power),
  MCP_emax = c(results_list_nsim_MCP_var[[13]]$sim_results$power, results_list_nsim_MCP_var[[14]]$sim_results$power, results_list_nsim_MCP_var[[15]]$sim_results$power, results_list_nsim_MCP_var[[16]]$sim_results$power, results_list_nsim_MCP_var[[17]]$sim_results$power,  results_list_nsim_MCP_var[[18]]$sim_results$power),
  Bay_linear = results_nsim_Bay_var$Bay_linear,
  Bay_exp = results_nsim_Bay_var$Bay_exponential,
  Bay_emax = results_nsim_Bay_var$Bay_emax
)

results_nsim_diff_var <- data.table(
  linear = results_nsim_var$Bay_linear - results_nsim_var$MCP_linear,
  exponential = results_nsim_var$Bay_exp - results_nsim_var$MCP_exp,
  emax = results_nsim_var$Bay_emax - results_nsim_var$MCP_emax,
  n_sim = n_sim_values
  )

results_nsim_diff_var <- melt(results_nsim_diff_var, id.vars = "n_sim")

plot_nsim_var <-  plot_power_deviation(results_nsim_diff_var, results_nsim_diff_var$n_sim, "nsim")
plot_nsim_var

setup
# save( results_min_MCP, minimal_Bay, results_monotonic_MCP, monotonic_Bay, results_non_monotonic_MCP, non_monotonic_Bay, results_var_MCP, variability_Bay, 
# results_min_MCP_nsample, minimal_nsample_Bay, results_monotonic_MCP_nsample, monotonic_nsample_Bay, results_non_monotonic_MCP_nsample, non_monotonic_nsample_Bay, results_var_MCP_nsample, 
# var_nsample_Bay, results_list_nsim_MCP, results_nsim_Bay, results_list_nsim_MCP_monotonic, results_nsim_Bay_monotonic, results_nsim_MCP_non_monotonic, results_nsim_Bay_non_monotonic, results_list_nsim_MCP_var, results_nsim_Bay_var
#      file = "simulation_data.RData")