adjustment based on function changes

NAMESPACE

@@ -9,7 +9,6 @@ S3method(print,postList)
 S3method(summary,postList)
 export(assessDesign)
 export(getBootstrapQuantiles)
-export(getBootstrapSamples)
 export(getContr)
 export(getCritProb)
 export(getESS)

R/bootstrapping.R

@@ -6,7 +6,9 @@
 #' This approach can be considered as the Bayesian equivalent of the frequentist bootstrap approach described in O'Quigley et al. (2017).
 #' Instead of drawing n bootstrap samples from the sampling distribution of the trial dose-response estimates, here the samples are directly taken from the posterior distribution.
 #' @references O'Quigley J, Iasonos A, Bornkamp B. 2017. Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials (1st ed.). Chapman and Hall/CRC. doi:10.1201/9781315151984
+#'
 #' @param model_fits An object of class modelFits, i.e. information about fitted models & corresponding model coefficients as well as the posterior distribution that was the basis for the model fitting
+#' @param quantiles A vector of quantiles that should be evaluated
 #' @param n_samples Number of samples that should be drawn as basis for the bootstrapped quantiles
 #' @param doses A vector of doses for which a prediction should be performed
 #' @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.
@@ -25,15 +27,17 @@
 #'                                simple      = TRUE)
 #'
 #' getBootstrapSamples(model_fits = fit,
+#'                     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
 #'
 #' @export
-getBootstrapSamples <- function (
+getBootstrapQuantiles <- function (
 
   model_fits,
+  quantiles,
   n_samples = 1e3,
   doses     = NULL,
   avg_fit   = TRUE
@@ -46,6 +50,7 @@ getBootstrapSamples <- function (
 
   dose_levels    <- model_fits[[1]]$dose_levels
   model_names    <- names(model_fits)
+  quantile_probs <- sort(unique(quantiles))
 
   if (is.null(doses)) {
 
@@ -92,20 +97,21 @@ getBootstrapSamples <- function (
 
   }
 
-  bs_samples <- as.data.frame(preds[order(rep(seq_along(model_names), length(doses))), ])
-  colnames(bs_samples) <- paste0("preds_", seq_len(n_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))
+
 
-  bs_samples_data <- cbind(
+  cr_bounds_data <- cbind(
     doses  = doses,
     models = rep(
       x    = factor(model_names, levels = model_names),
-      each = length(doses))
-  )
-    # as.data.frame(preds))
+      each = length(doses)),
+    as.data.frame(quant_mat))
 
-  # tidyr::pivot_longer(bs_samples_data, cols = contains("preds"))
-
-  return (bs_samples_data)
+  return (cr_bounds_data)
 
 }
 
@@ -138,26 +144,26 @@ getBootstrapSamples <- function (
 #' @return  A data frame with entries doses, models, and quantiles
 #'
 #' @export
-getBootstrapQuantiles <- function (
-
-  bs_samples,
-  quantiles
-
-) {
-
-  quantile_probs <- sort(unique(quantiles))
-
-  bs_quantiles <- t(apply(X      = bs_samples[, -c(1, 2)],
-                          MARGIN = 1,
-                          FUN    = stats::quantile,
-                          probs  = quantile_probs))
-
-  bs_quantiles_data <- cbind(
-    bs_samples[, c(1, 2)],
-    as.data.frame(bs_quantiles))
-
-  return (bs_quantiles_data)
-
-}
+# getBootstrapQuantiles <- function (
+#
+#   bs_samples,
+#   quantiles
+#
+# ) {
+#
+#   quantile_probs <- sort(unique(quantiles))
+#
+#   bs_quantiles <- t(apply(X      = bs_samples[, -c(1, 2)],
+#                           MARGIN = 1,
+#                           FUN    = stats::quantile,
+#                           probs  = quantile_probs))
+#
+#   bs_quantiles_data <- cbind(
+#     bs_samples[, c(1, 2)],
+#     as.data.frame(bs_quantiles))
+#
+#   return (bs_quantiles_data)
+#
+# }
 
 

R/plot.R

@@ -1,11 +1,11 @@
 #' @title plot.modelFits
-#' 
+#'
 #' @description Plot function based on the ggplot2 package. Providing visualizations for each model and a average Fit.
 #' Black lines show the fitted dose response models and an AIC based average model.
-#' Dots indicate the posterior median and vertical lines show corresponding credible intervals (i.e. the variability of the posterior distribution of the respective dose group). 
+#' Dots indicate the posterior median and vertical lines show corresponding credible intervals (i.e. the variability of the posterior distribution of the respective dose group).
 #' To assess the uncertainty of the model fit one can in addition visualize credible bands (default coloring as orange shaded areas).
 #' The calculation of these bands is performed via the getBootstrapQuantiles() function.
-#' The default setting is that these credible bands are not calculated. 
+#' The default setting is that these credible bands are not calculated.
 #' @param x An object of type modelFits
 #' @param gAIC Logical value indicating whether gAIC values are shown in the plot. Default TRUE
 #' @param avg_fit Logical value indicating whether average fit is presented in the plot. Default TRUE
@@ -19,7 +19,7 @@
 #' @examples
 #' posterior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2),
 #'                        DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2),
-#'                        DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) ,  
+#'                        DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) ,
 #'                        DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) ,
 #'                        DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1), sigma = 2))
 #' models <- c("exponential", "linear")
@@ -28,12 +28,12 @@
 #'                     posterior   = posterior_list,
 #'                     dose_levels = dose_levels,
 #'                     simple      = TRUE)
-#'                            
-#' plot(fit) 
+#'
+#' plot(fit)
 #' @return A ggplot2 object
 #' @export
 plot.modelFits <- function (
-    
+
   x,
   gAIC      = TRUE,
   avg_fit   = TRUE,
@@ -44,12 +44,11 @@ plot.modelFits <- function (
   n_bs_smpl = 1e3,
   acc_color = "orange",
   ...
-  
+
 ) {
-
   ## R CMD --as-cran appeasement
   .data <- NULL
-  
+
   checkmate::check_class(x, "modelFits")
   checkmate::check_logical(gAIC)
   checkmate::check_logical(avg_fit)
@@ -59,44 +58,44 @@ plot.modelFits <- function (
   checkmate::check_double(alpha_CrB, lower = 0, upper = 1, len = 2)
   checkmate::check_integer(n_bs_smpl, lower = 1, upper = Inf)
   checkmate::check_string(acc_color, na.ok = TRUE)
-  
+
   plot_res   <- 1e2
   model_fits <- x
-  
+
   dose_levels  <- model_fits[[1]]$dose_levels
   post_summary <- summary.postList(
     object = attr(model_fits, "posterior"),
     probs  = c(alpha_CrI / 2, 0.5, 1 - alpha_CrI / 2))
   dose_seq <- seq(from       = min(dose_levels),
                   to         = max(dose_levels),
                   length.out = plot_res)
-  
+
   preds_models <- sapply(model_fits, predictModelFit, doses = dose_seq)
   model_names  <- names(model_fits)
-  
+
   if (avg_fit) {
-    
+
     mod_weights <- sapply(model_fits, function (x) x$model_weight)
     avg_preds   <- preds_models %*% mod_weights
-    
+
     preds_models <- cbind(preds_models, avg_preds)
     model_names  <- c(model_names, "avgFit")
-    
+
   }
-  
+
   gg_data <- data.frame(
     doses  = rep(dose_seq, length(model_names)),
     fits   = as.vector(preds_models),
     models = rep(factor(model_names, levels = model_names),
                  each = plot_res))
-  
+
   if (gAIC) {
-    
+
     g_AICs     <- sapply(model_fits, function (x) x$gAIC)
     label_gAUC <- paste("AIC:", round(g_AICs, digits = 1))
-    
+
     if (avg_fit) {
-      
+
       mod_weights <- sort(mod_weights, decreasing = TRUE)
       paste_names <- names(mod_weights) |>
         gsub("exponential", "exp", x = _) |>
@@ -106,88 +105,84 @@ plot.modelFits <- function (
         gsub("sigEmax",    "sigE", x = _) |>
         gsub("betaMod",    "betaM", x = _)|>
         gsub("quadratic",    "quad", x = _)
-      
+
       label_avg  <- paste0(paste_names, "=", round(mod_weights, 1),
                            collapse = ", ")
       label_gAUC <- c(label_gAUC, label_avg)
-      
+
     }
-    
+
   }
-  
+
   if (cr_bands) {
-    
-    bs_samples <- getBootstrapSamples(
+
+    crB_data <- getBootstrapQuantiles(
       model_fits = model_fits,
       n_samples  = n_bs_smpl,
+      quantiles  = c(0.5, sort(unique(c(alpha_CrB / 2, 1 - alpha_CrB / 2)))),
       avg_fit    = avg_fit,
       doses      = dose_seq)
-
-    crB_data <- getBootstrapQuantiles(
-      bs_samples = bs_samples,
-      quantiles  = c(0.5, sort(unique(c(alpha_CrB / 2, 1 - alpha_CrB / 2)))))
-    colnames(crB_data) <- c("models", colnames(crB_data)[-1])
-
+
     getInx <- function (alpha_CrB) {
       n        <- length(alpha_CrB)
-      inx_list <- lapply(seq_len(n), function (i) c(i, 2 * n - i + 2) + 2)
+      inx_list <- lapply(seq_len(n), function (i) c(i, 2 * n - i + 1) + 3)
       return (inx_list)}
-    
+
   }
-  
-  plts <- ggplot2::ggplot() + 
+
+  plts <- ggplot2::ggplot() +
     ## Layout etc.
-    ggplot2::theme_bw() + 
+    ggplot2::theme_bw() +
     ggplot2::labs(x = "Dose",
                   y = "Model Fits") +
     ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
                    panel.grid.minor = ggplot2::element_blank()) +
     ## gAIC
     {if (gAIC) {
-      
+
       ggplot2::geom_text(
         data    = data.frame(
           models = unique(gg_data$models),
           label  = label_gAUC),
         mapping = ggplot2::aes(label = label_gAUC),
         x = -Inf, y = Inf, hjust = "inward", vjust = "inward",
         size = 3)
-      
+
     }
     }
-  
+
   if (cr_bands) {
-    
+
     ## Bootstrapped Credible Bands
     for (inx in getInx(alpha_CrB)) {
-      
+
       loop_txt <- paste0(
         "ggplot2::geom_ribbon(
           data    = crB_data,
-          mapping = ggplot2::aes(x    = .data$dose,
+          mapping = ggplot2::aes(x    = doses,
                                  ymin = crB_data[, ", inx[1], "],
                                  ymax = crB_data[, ", inx[2], "]),
-          fill    = acc_color,                    
+          fill    = acc_color,
           alpha   = 0.25)")
-      
+
       plts <- plts + eval(parse(text = loop_txt))
-      
+
     }
     rm(inx)
-    
+
     ## Bootstrapped Fit
     plts <- plts +
       ggplot2::geom_line(
         data    = crB_data,
-        mapping = ggplot2::aes(.data$dose, .data$`50%`),
+        mapping = ggplot2::aes(.data$doses, .data$`50%`),
         color   = acc_color)
-    
+
   }
-  
+
   plts <- plts +
     ## Posterior Credible Intervals
     {if (cr_intv) {
-      
+
       ggplot2::geom_errorbar(
         data    = data.frame(x    = dose_levels,
                              ymin = post_summary[, 3],
@@ -197,9 +192,9 @@ plot.modelFits <- function (
                                ymax = .data$ymax),
         width   = 0,
         alpha   = 0.5)
-      
+
     }
-    } + 
+    } +
     ## Posterior Medians
     ggplot2::geom_point(
       data    = data.frame(
@@ -210,10 +205,10 @@ plot.modelFits <- function (
     ## Fitted Models
     ggplot2::geom_line(
       data    = gg_data,
-      mapping = ggplot2::aes(.data$doses, .data$fits)) + 
+      mapping = ggplot2::aes(.data$doses, .data$fits)) +
     ## Faceting
     ggplot2::facet_wrap(~ models)
-  
+
   return (plts)
-  
-}
+
+}