R/assessDesign.R

#' @title assessDesign
#'
#' @description This function performs simulation based trial design evaluations for a set of specified dose-response models
#'
#' @param n_patients Vector specifying the planned number of patients per dose group. A minimum of 2 patients are required in each group.
#' @param mods An object of class `Mods` as specified in the `DoseFinding` package.
#' @param prior_list A prior_list object specifying the utilized prior for the different dose groups
#' @param sd A positive value, specification of assumed sd. Not required if either `data_sim` or `estimates_sim` is provided. Also not required in case of binary endpoint. Default NULL
#' @param contr An object of class `optContr` as created by the `DoseFinding::getContr` function. Allows specification of a fixed contrasts matrix. Default NULL.
#' @param dr_means A vector, allows specification of individual (not model based) assumed effects per dose group. Default NULL.
#' @param data_sim An optional data frame for custom simulated data. Must follow the data structure as provided by `simulateData()`. Default NULL.
#' @param estimates_sim An optional named list of 1) list of vectors for the estimated means per dose group (`estimates_sim$mu_hats`) and 2) of list of matrices for the covariance matrices specifying the (estimated) variabilities (`estimates_sim$S_hats`). Dimensions of entries must match the number of dose levels. Default NULL.
#' @param n_sim Number of simulations to be performed
#' @param alpha_crit_val (Un-adjusted) Critical value to be used for the MCP testing step. Passed to the `getCritProb` function for the calculation of adjusted critical values (on the probability scale). Default 0.05.
#' @param modeling Boolean variable defining whether the Mod part of Bayesian MCP-Mod will be performed in the assessment. More heavy on resources. Default FALSE.
#' @param simple Boolean variable defining whether simplified fit will be applied, see `?getModelFits`. Set automatically to TRUE if argument `delta` is provided. Passed to `getModelFits` Default TRUE.
#' @param avg_fit Boolean variable, defining whether an average fit (based on generalized AIC weights) should be performed in addition to the individual models. Default TRUE.
#' @param reestimate Boolean variable defining whether critical value should be calculated with re-estimated contrasts (see `getCritProb` function for more details). Default FALSE.
#' @param delta A numeric value for the threshold Delta for the MED assessment. If NULL, no MED assessment is performed. Default NULL.
#' @param evidence_level A numeric value between 0 and 1 for the evidence level gamma for the MED assessment. Only required for Bayesian MED assessment, see `?getMED` for details. If NULL, MED assessment will be performed on the fitted model according to the argument `med_selection`.  Default NULL.
#' @param n_bs_samples Number of bootstrap samples for the MED assessment if `evidence_level` is provided. Note that more extreme quantiles (i.e., quantiles closer to 0 or 1) tend to require more bootstrap samples to maintain precision. Default 1000.
#' @param med_selection A string, either `"avgFit"` or `"bestFit"`, for the method of MED selection. Default `"avgFit"`.
#' @param probability_scale A boolean to specify if the trial has a continuous or a binary outcome. Setting to TRUE will transform calculations from the logit scale to the probability scale, which can be desirable for a binary outcome. Default FALSE.
#'
#' @return Returns success probabilities for the different assumed dose-response shapes, attributes also includes information around average success rate (across all assumed models) and prior Effective sample size.
#'
#' @examples
#' mods <- DoseFinding::Mods(linear      = NULL,
#'                           emax        = c(0.5, 1.2),
#'                           exponential = 2,
#'                           betaMod     = c(1, 1),
#'                           doses       = c(0, 0.5, 2,4, 8),
#'                           maxEff      = 6)
#'                           
#' sd <- 12
#' prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2),
#'                    DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
#'                    DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,
#'                    DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
#'                    DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2))
#' n_patients <- c(40, 60, 60, 60, 60)
#' dose_levels  <- c(0, 0.5, 2, 4, 8)
#' 
#' success_probabilities <- assessDesign(
#'   n_patients  = n_patients,
#'   mods        = mods,
#'   prior_list  = prior_list,
#'   sd          = sd,
#'   n_sim       = 1e2) # speed up example run time
#'
#' success_probabilities
#' 
#' ## Analysis with custom data
#' data_sim <- simulateData(
#'   n_patients        = n_patients,
#'   dose_levels       = dose_levels,
#'   sd                = sd,
#'   mods              = mods,
#'   n_sim             = 10)
#'
#' success_probabilities_cd <- assessDesign(
#'   n_patients  = n_patients,
#'   mods        = mods,
#'   prior_list  = prior_list,
#'   data_sim    = data_sim,
#'   sd          = sd,
#'   n_sim       = 1e2) # speed up example run time
#'
#' success_probabilities_cd
#' 
#' ## Analysis with custom dose response relationship
#' custom_dr_means <- c(1, 2, 3, 4, 5)
#' 
#' success_probs_custom_dr <- assessDesign(
#'   n_patients  = n_patients,
#'   mods        = mods,
#'   prior_list  = prior_list,
#'   dr_means    = custom_dr_means,
#'   sd          = sd,
#'   n_sim       = 1e2) # speed up example run time
#'
#' success_probs_custom_dr
#' 
#' ## Analysis with custom estimates for means and variabilies
#' ## No simulated data, only simulated model estimates
#' estimates_sim <- list(mu_hats = replicate(100, list(c(1, 2, 3, 4, 5) + rnorm(5, 0, 1))),
#'                       S_hats  = list(diag(1, 5)))
#' 
#' success_probs_custom_est <- assessDesign(
#'   n_patients    = n_patients,
#'   mods          = mods,
#'   prior_list    = prior_list,
#'   estimates_sim = estimates_sim)
#'
#' success_probs_custom_est
#'
#' if (interactive()) { # takes typically > 5 seconds
#'
#' # with MED estimation without bootstrapping
#' # see ?getMED for details
#'
#' success_probabilities <- assessDesign(
#'   n_patients     = n_patients,
#'   mods           = mods,
#'   prior_list     = prior_list,
#'   sd             = sd,
#'   modeling       = TRUE,
#'   n_sim          = 10, # speed up example run time
#'   delta          = 7)
#'
#'   success_probabilities
#'
#' # with MED estimation with bootstrapping
#'
#' success_probabilities <- assessDesign(
#'   n_patients     = n_patients,
#'   mods           = mods,
#'   prior_list     = prior_list,
#'   sd             = sd,
#'   modeling       = TRUE,
#'   n_sim          = 10, # speed up example run time
#'   delta          = 7,
#'   evidence_level = 0.8)
#'
#'   success_probabilities
#'
#' }
#'
#' @export
assessDesign <- function (
    
  n_patients,
  mods,
  prior_list,
  
  sd                = NULL,
  
  contr             = NULL,
  
  dr_means          = NULL,
  data_sim          = NULL,
  estimates_sim     = NULL,
  
  n_sim             = 1e3,
  alpha_crit_val    = 0.05,
  modeling          = FALSE,
  simple            = TRUE,
  avg_fit           = TRUE,
  reestimate        = FALSE,
  
  delta             = NULL,
  evidence_level    = NULL,
  n_bs_samples      = 1000,
  med_selection     = c("avgFit", "bestFit"),
  
  probability_scale = FALSE
  
) {
  
  checkmate::assert_vector(n_patients, len = length(attr(mods, "doses")), any.missing = FALSE)
  checkmate::assert_double(n_patients, lower = 2, upper = Inf)
  checkmate::assert_class(mods, classes = "Mods")
  checkmate::assert_list(prior_list, names = "named", len = length(attr(mods, "doses")), any.missing = FALSE)
  checkmate::assert_list(estimates_sim, null.ok = TRUE)
  # sensitive to how DoseFinding labels their attributes for "Mods" class
  checkmate::assert_integerish(n_sim, lower = 1, upper = Inf)
  checkmate::assert_double(alpha_crit_val, lower = 0, upper = 1)
  checkmate::assert_flag(modeling)
  checkmate::assert_flag(probability_scale)
  
  # TODO: check that prior_list has 'sd_tot' attribute, and that it's numeric # this is not applicable at the moment
  
  stopifnot("Cannot specify both 'data_sim' and 'estimates_sim' simultaneously" = 
              is.null(data_sim)  & is.null(estimates_sim) |
              !is.null(data_sim) & is.null(estimates_sim) |
              is.null(data_sim)  & !is.null(estimates_sim))
  
  modeling <- ifelse(!is.null(delta), TRUE, modeling)
  
  dose_levels <- attr(mods, "doses")
  
  if (!is.null(estimates_sim)) {
    
    stopifnot("'estimates_sim' must be a named list of lenght 2 with list items called 'mu_hats' and 'S_hats'" =
                length(estimates_sim) ==  2L &
                all(names(estimates_sim) %in% c("mu_hats", "S_hats")))
    stopifnot("The length of 'estimates_sim$S_hats' must match 'estimates_sim$mu_hats' or 1." =
                length(estimates_sim$S_hats) == length(estimates_sim$mu_hats) |
                length(estimates_sim$S_hats) == 1L)
    stopifnot("The length of each entry of 'estimates_sim$mu_hats' must match the number of dose levels." =
                length(estimates_sim$mu_hats[[1]]) == length(dose_levels))
    stopifnot("The dimensions of each entry of 'estimates_sim$S_hats' must match the number of dose levels." =
                ncol(estimates_sim$S_hats[[1]]) == length(dose_levels))
    
  } else if (!is.null(data_sim)) {
    
    ## lazily simulate data anyway ...
    data <- simulateData(
      n_patients        = n_patients,
      dose_levels       = dose_levels,
      sd                = 1,
      mods              = mods,
      n_sim             = floor(nrow(data_sim) / sum(n_patients)), # adjust the number of simulations to data_sim
      dr_means          = dr_means,
      probability_scale = probability_scale)
    
    ## ... and then check if formats of data_sim and data match
    stopifnot("'data_sim' must follow the data structure as provided by simulateData()" =
                isTRUE(all.equal(data_sim[, 1:3], data[, 1:3])))
    stopifnot("data_sim must have named columns" = !is.null(colnames(data_sim)))
    checkmate::assert_data_frame(data_sim, types = "numeric", any.missing = FALSE)
    
    data <- data_sim
    
  } else {
    
    data <- simulateData(
      n_patients        = n_patients,
      dose_levels       = dose_levels,
      sd                = sd,
      mods              = mods,
      n_sim             = n_sim,
      dr_means          = dr_means,
      probability_scale = probability_scale)
    
  }
  
  crit_prob_adj <- getCritProb(
    mods           = mods,
    dose_levels    = dose_levels,
    dose_weights   = n_patients,
    alpha_crit_val = alpha_crit_val)
  
  if (!reestimate & is.null(contr)) {
    
    if (!is.null(data_sim) | !is.null(estimates_sim)) {
      
      message("Consider to provide 'contr' for your custom simulated data or analysis results.")
      
    }
    
    contr <- getContr(
      mods           = mods,
      dose_levels    = dose_levels,
      dose_weights   = n_patients,
      prior_list     = prior_list)
    
  }
  
  # get model names of true underlying models
  if (is.null(estimates_sim))  {
    
    model_names <- colnames(data)[-c(1, 2, 3)]
    
  } else {
    
    model_names <- "estimates_sim"
    
  }
  
  eval_design <- lapply(model_names, function (model_name) {
    
    if (is.null(estimates_sim)) {
      
      posterior_list <- getPosterior(
        data              = getModelData(data, model_name),
        prior_list        = prior_list,
        probability_scale = probability_scale)
      
    } else {
      
      posterior_list <- mapply(getPosterior,
             mu_hat   = estimates_sim$mu_hats,
             S_hat    = estimates_sim$S_hats,
             MoreArgs = list(prior_list = prior_list),
             SIMPLIFY = FALSE)
      
      names(posterior_list) <- seq_along(estimates_sim$mu_hats)
      
    }
    
    if (reestimate & is.null(contr)) {
      
      post_sds <- sapply(posterior_list, function (post) summary(post)[, 2])
      
      contr <- apply(post_sds, 2, function (post_sd) getContr(
        mods          = mods,
        dose_levels   = dose_levels,
        cov_posterior = diag(post_sd^2)))
      
    }
    
    if (modeling) {
      
      attr(contr, "direction") <- attr(mods,"direction")
      
      b_mcp_mod <- performBayesianMCPMod(
        posterior_list    = posterior_list,
        contr             = contr,
        crit_prob_adj     = crit_prob_adj,
        simple            = simple,
        avg_fit           = avg_fit,
        delta             = delta,
        evidence_level    = evidence_level,
        med_selection     = med_selection,
        n_samples         = n_bs_samples,
        probability_scale = probability_scale) # should this also be transferred via attribute?
      
    } else {
      
      b_mcp <- performBayesianMCP(
        posterior_list = posterior_list,
        contr          = contr,
        crit_prob_adj  = crit_prob_adj)
      
    }
    
  })
  
  names(eval_design) <- model_names
  
  avg_success_rate <- mean(sapply(eval_design, function (x) {
    
    ifelse(identical(modeling, FALSE),
           attr(x, "successRate"),
           attr(x$BayesianMCP, "successRate"))
    
  }))
  
  attr(eval_design, "avgSuccessRate") <- avg_success_rate
  
  if (modeling & !is.null(delta)) {
    
    avg_med_id_rate <- mean(sapply(eval_design, function (x) {
      
      mean(attr(x, "MED")[, "med_reached"])
      
    }))
    
    attr(eval_design, "avgMEDIdentificationRate") <- avg_med_id_rate
    
  }
  
  attr(eval_design, "placEff")        <- ifelse(test = is.null(dr_means),
                                                yes  = attr(mods, "placEff"),
                                                no   = dr_means[1])
  attr(eval_design, "maxEff")         <- ifelse(test = is.null(dr_means),
                                                yes  = attr(mods, "maxEff"), # this is on logit scale in case of binary endpoint
                                                no   = diff(range(dr_means)))
  attr(eval_design, "sampleSize")     <- n_patients
  attr(eval_design, "priorESS")       <- round(getESS(prior_list), 1)
  
  return (eval_design)
  
}

Try the BayesianMCPMod package in your browser

Any scripts or data that you put into this service are public.

BayesianMCPMod documentation built on Feb. 23, 2026, 5:06 p.m.