- added parallelization with future.apply for bootstrapping and

DESCRIPTION

@@ -34,6 +34,7 @@ Depends:
 Imports:
     checkmate,
     DoseFinding (>= 1.1-1),
+    future.apply,
     ggplot2,
     methods,
     nloptr,

R/BMCPMod.R

@@ -1,3 +1,24 @@
+addSignificance <- function (
+
+  model_fits,
+  sign_models
+
+) {
+
+  names(sign_models) <- NULL
+
+  model_fits_out <- lapply(seq_along(model_fits), function (i) {
+
+    c(model_fits[[i]], significant = sign_models[i])
+
+  })
+
+  attributes(model_fits_out) <- attributes(model_fits)
+
+  return (model_fits_out)
+
+}
+
 #' @title assessDesign
 #'
 #' @description This function performs simulation based trial design evaluations for a set of specified dose-response models
@@ -38,9 +59,19 @@
 #'   mods        = mods,
 #'   prior_list  = prior_list,
 #'   sd          = sd,
-#'   n_sim       = 1e2) # speed up exammple run time
+#'   n_sim       = 1e2) # speed up example run time
 #'
 #' success_probabilities
+#' 
+#' success_probabilities <- assessDesign(
+#'   n_patients  = n_patients,
+#'   mods        = mods,
+#'   prior_list  = prior_list,
+#'   sd          = sd,
+#'   modeling    = TRUE,
+#'   n_sim       = 1e2) # speed up example run time
+#'   
+#'   success_probabilities
 #'
 #' }
 #'
@@ -119,7 +150,7 @@ assessDesign <- function (
         sd_posterior   = post_sd))
 
     }
-    
+
     if (identical(modeling, FALSE)) {
 
       b_mcp <- performBayesianMCP(
@@ -163,6 +194,49 @@ assessDesign <- function (
 
 }
 
+BayesMCPi <- function (
+
+  posterior_i,
+  contr,
+  crit_prob_adj
+
+) {
+
+  getPostProb <- function (
+
+    contr_j,     # j: dose level
+    post_combs_i # i: simulation outcome
+
+  ) {
+
+    ## Test statistic = sum over all components of
+    ## posterior weight * normal probability distribution of
+    ## critical values for doses * estimated mean / sqrt(product of critical values for doses)
+
+    ## Calculation for each component of the posterior
+    contr_theta   <- apply(post_combs_i$means, 1, `%*%`, contr_j)
+    contr_var     <- apply(post_combs_i$vars, 1, `%*%`, contr_j^2)
+    contr_weights <- post_combs_i$weights
+
+    ## P(c_m * theta > 0 | Y = y) for a shape m (and dose j)
+    post_probs <- sum(contr_weights * stats::pnorm(contr_theta / sqrt(contr_var)))
+
+    return (post_probs)
+
+  }
+
+  post_combs_i <- getPostCombsI(posterior_i)
+  post_probs   <- apply(contr$contMat, 2, getPostProb, post_combs_i)
+
+  res <- c(sign          = ifelse(max(post_probs) > crit_prob_adj, 1, 0),
+           crit_prob_adj = crit_prob_adj,
+           max_post_prob = max(post_probs),
+           post_probs    = post_probs)
+
+  return (res)
+
+}
+
 #' @title getContr
 #'
 #' @description This function calculates contrast vectors that are optimal for detecting certain alternatives via applying the function optContr() of the DoseFinding package.
@@ -339,6 +413,22 @@ getCritProb <- function (
 
 }
 
+getModelSuccesses <- function (b_mcp) {
+
+  stopifnot(inherits(b_mcp, "BayesianMCP"))
+
+  model_indices <- grepl("post_probs.", colnames(b_mcp))
+  model_names   <- colnames(b_mcp)[model_indices] |>
+    sub(pattern = "post_probs.", replacement = "", x = _)
+
+  model_successes <- colMeans(b_mcp[, model_indices] > b_mcp[, "crit_prob_adj"])
+
+  names(model_successes) <- model_names
+
+  return (model_successes)
+
+}
+
 #' @title performBayesianMCPMod
 #'
 #' @description Performs Bayesian MCP Test step and modeling in a combined fashion. See performBayesianMCP() function for MCP Test step and getModelFits() for the modeling step
@@ -423,57 +513,21 @@ performBayesianMCPMod <- function (
     posterior_list = posterior_list,
     contr          = contr,
     crit_prob_adj  = crit_prob_adj)
+
+  b_mod <- performBayesianMod(
+    b_mcp          = b_mcp,
+    posterior_list = posterior_list,
+    model_shapes   = model_shapes,
+    dose_levels    = dose_levels,
+    simple         = simple)
 
-  fits_list <- lapply(seq_along(posterior_list), function (i) {
-
-    if (b_mcp[i, 1]) {
-
-      sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "crit_prob_adj")
-
-      model_fits  <- getModelFits(
-        models      = model_shapes,
-        dose_levels = dose_levels,
-        posterior   = posterior_list[[i]],
-        simple      = simple)
-
-      model_fits  <- addSignificance(model_fits, sign_models)
-
-    } else {
-
-      NULL
-
-    }
-
-  })
-
-  bmcpmod        <- list(BayesianMCP = b_mcp, Mod = fits_list)
+  bmcpmod        <- list(BayesianMCP = b_mcp, Mod = b_mod)
   class(bmcpmod) <- "BayesianMCPMod"
 
   return (bmcpmod)
 
 }
 
-addSignificance <- function (
-
-  model_fits,
-  sign_models
-
-) {
-
-  names(sign_models) <- NULL
-
-  model_fits_out <- lapply(seq_along(model_fits), function (i) {
-
-    c(model_fits[[i]], significant = sign_models[i])
-
-  })
-
-  attributes(model_fits_out) <- attributes(model_fits)
-
-  return (model_fits_out)
-
-}
-
 #' @title performBayesianMCP
 #'
 #' @description Performs Bayesian MCP Test step, as described in Fleischer et al. (2022).
@@ -565,61 +619,36 @@ performBayesianMCP <- function(
 
 }
 
-getModelSuccesses <- function (b_mcp) {
-
-  stopifnot(inherits(b_mcp, "BayesianMCP"))
-
-  model_indices <- grepl("post_probs.", colnames(b_mcp))
-  model_names   <- colnames(b_mcp)[model_indices] |>
-    sub(pattern = "post_probs.", replacement = "", x = _)
-
-  model_successes <- colMeans(b_mcp[, model_indices] > b_mcp[, "crit_prob_adj"])
-
-  names(model_successes) <- model_names
-
-  return (model_successes)
-
-}
-
-BayesMCPi <- function (
-
-  posterior_i,
-  contr,
-  crit_prob_adj
-
+performBayesianMod <- function (
+
+  b_mcp,
+  posterior_list,
+  model_shapes,
+  dose_levels,
+  simple
+
 ) {
-
-  getPostProb <- function (
-
-    contr_j,     # j: dose level
-    post_combs_i # i: simulation outcome
-
-  ) {
-
-    ## Test statistic = sum over all components of
-    ## posterior weight * normal probability distribution of
-    ## critical values for doses * estimated mean / sqrt(product of critical values for doses)
-
-    ## Calculation for each component of the posterior
-    contr_theta   <- apply(post_combs_i$means, 1, `%*%`, contr_j)
-    contr_var     <- apply(post_combs_i$vars, 1, `%*%`, contr_j^2)
-    contr_weights <- post_combs_i$weights
-
-    ## P(c_m * theta > 0 | Y = y) for a shape m (and dose j)
-    post_probs <- sum(contr_weights * stats::pnorm(contr_theta / sqrt(contr_var)))
-
-    return (post_probs)
-
-  }
-
-  post_combs_i <- getPostCombsI(posterior_i)
-  post_probs   <- apply(contr$contMat, 2, getPostProb, post_combs_i)
-
-  res <- c(sign          = ifelse(max(post_probs) > crit_prob_adj, 1, 0),
-           crit_prob_adj = crit_prob_adj,
-           max_post_prob = max(post_probs),
-           post_probs    = post_probs)
-
-  return (res)
-
+
+  fits_list <- future.apply::future_lapply(seq_along(posterior_list), function (i) {
+
+    if (b_mcp[i, 1]) {
+
+      sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "crit_prob_adj")
+
+      model_fits  <- getModelFits(
+        models      = model_shapes,
+        dose_levels = dose_levels,
+        posterior   = posterior_list[[i]],
+        simple      = simple) 
+
+      model_fits  <- addSignificance(model_fits, sign_models)
+
+    } else {
+
+      NULL
+
+    }
+
+  })
+
 }

R/bootstrapping.R

@@ -27,9 +27,9 @@
 #'                                simple      = TRUE)
 #'
 #' getBootstrapQuantiles(model_fits = fit,
-#'                     quantiles  = c(0.025, 0.5, 0.975),
-#'                     n_samples  = 10, # speeding up example run time
-#'                     doses      = c(0, 6, 8))
+#'                       quantiles  = c(0.025, 0.5, 0.975),
+#'                       n_samples  = 10, # speeding up example run time
+#'                       doses      = c(0, 6, 8))
 #'
 #' @return  A data frame with columns for model, dose, and bootstrapped samples
 #'
@@ -57,11 +57,12 @@ getBootstrapQuantiles <- function (
     doses <- seq(min(dose_levels), max(dose_levels), length.out = 100L)
 
   }
-
-  preds <- apply(mu_hat_samples, 1, function (mu_hat) {
-
+
+  # predictions[(dose1 model1, dose1 model2, ..., dose2 model1, ...), samples]
+  preds <- future.apply::future_apply(mu_hat_samples, 1, function (mu_hat) {
+
     preds_mu_hat <- sapply(model_names, function (model) {
-
+      
       fit <- DoseFinding::fitMod(
         dose  = model_fits[[1]]$dose_levels,
         resp  = mu_hat,
@@ -70,25 +71,26 @@ getBootstrapQuantiles <- function (
         type  = "general",
         bnds  = DoseFinding::defBnds(
           mD = max(model_fits[[1]]$dose_levels))[[model]])
-
+      
       preds <- stats::predict(fit, doseSeq = doses, predType = "ls-means")
       attr(preds, "gAIC") <- DoseFinding::gAIC(fit)
-
+      
       return (preds)
-
+      
     }, simplify = FALSE)
-
+    
     preds_mu_mat <- do.call(rbind, preds_mu_hat)
-
+    
     if (avg_fit) {
-
+      
       avg_fit_indx <- which.min(sapply(preds_mu_hat, attr, "gAIC"))
       preds_mu_mat <- rbind(preds_mu_mat, avgFit = preds_mu_mat[avg_fit_indx, ])
-
+      
     }
-
+
+    # predictions[models, doses]
     return (preds_mu_mat)
-
+    
   })
 
   if (avg_fit) {
@@ -97,13 +99,13 @@ getBootstrapQuantiles <- function (
 
   }
 
+  # predictions[(dose1 model1, dose2 model1, ..., dose1 model2, ...), samples]
   sort_indx <- order(rep(seq_along(model_names), length(doses)))
   quant_mat <- t(apply(X      = preds[sort_indx, ],
                        MARGIN = 1,
                        FUN    = stats::quantile,
                        probs  = quantile_probs))
 
-
   cr_bounds_data <- cbind(
     doses  = doses,
     models = rep(