Merge branch 'main' into negative_effect_MED_fixes

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: BayesianMCPMod
 Title: Simulate, Evaluate, and Analyze Dose Finding Trials with Bayesian MCPMod
-Version: 1.1.0.3
+Version: 1.2.0
 Authors@R: c(
     person("Boehringer Ingelheim Pharma GmbH & Co. KG", role = c("cph", "fnd")),
     person("Stephan", "Wojciekowski", , "[email protected]", role = c("aut", "cre")),

NEWS.md

@@ -1,3 +1,7 @@
+## BayesianMCPMod 1.2.0 (TODO: DATE)
+
+- 
+
 ## BayesianMCPMod 1.1.0 (07-Mar-2025)
 
 - Fixed a bug in plot.modelFits() that would plot credible bands based on incorrectly selected bootstrapped quantiles

R/bootstrapping.R

@@ -74,6 +74,8 @@ getBootstrapQuantiles <- function (
                     sample_type == "sample" ~ "abs",
                     sample_type == "sample_diff" ~ "diff"))
 
+  attr(bs_quantiles, "direction") <- attr(model_fits, "direction")
+
   return (bs_quantiles)
 
 }

R/modeling.R

@@ -18,7 +18,8 @@
 #' }
 #' where \eqn{Q} denotes the number of models included in the averaging procedure.
 #' @references Schorning K, Bornkamp B, Bretz F, Dette H. 2016. Model selection versus model averaging in dose finding studies. Stat Med; 35; 4021-4040.
-#' @param models List (or vector) of model names for which a fit will be performed.
+#' @param models A Mods object as created with `DoseFinding::Mods()` or a vector
+#' of model names for which a fit will be performed.
 #' Implemented model shapes are `"linear"`, `"exponential"`, `"logistic"`,
 #' `"emax"`, `"sigEmax"`, `"quadratic"`, and `"betaMod"`.
 #' @param dose_levels A vector containing the different dosage levels.
@@ -60,24 +61,31 @@ getModelFits <- function (
 
 ) {
 
-  if (inherits(models, "character")) {
-    models <- stats::setNames(as.list(models), models)
-  } else if (inherits(models, "Mods")) {
-    models <- stats::setNames(as.list(names(models)), names(models))
-  }
-
-  checkmate::assert_list(models, any.missing = FALSE, types = "character")
+  checkmate::assert(
+    checkmate::check_class(models, "Mods"),
+    checkmate::check_character(models, any.missing = FALSE)
+  )
   checkmate::assert_double(dose_levels, lower = 0, any.missing = FALSE)
   checkmate::assert(
     checkmate::check_class(posterior, "postList"),
     checkmate::check_list(posterior) && all(sapply(posterior, inherits, c("normMix", "mix")))
   )
   checkmate::assert_logical(avg_fit)
   checkmate::assert_logical(simple)
+
+  if (inherits(models, "character")) {
+
+    model_names <- models
+
+  } else if (inherits(models, "Mods")) {
+
+    model_names <- names(models)
+
+  }
+
+  model_names <- sort(unique(gsub("\\d", "", model_names)))
 
-  model_names <- sort(unique(gsub("\\d", "", names(models))))
-
-  getModelFit <- ifelse(simple, getModelFitSimple, getModelFitOpt)
+  getModelFit <- if (simple) getModelFitSimple else getModelFitOpt
 
   model_fits  <- lapply(model_names, getModelFit, dose_levels, posterior, list("scal" = attr(models, "scal")))
   model_fits  <- addModelWeights(model_fits)
@@ -90,13 +98,57 @@ getModelFits <- function (
 
   }
 
+  attr(model_fits, "direction") <- getDirection(models, model_fits)
   attr(model_fits, "posterior") <- posterior
   class(model_fits)             <- "modelFits"
 
   return (model_fits)
 
 }
 
+getDirection <- function (
+
+  models,
+  model_fits
+
+) {
+
+  if (inherits(models, "Mods")) {
+
+    direction   <- attr(models, "direction")
+
+  } else if (inherits(models, "character")) {
+
+    ## Try guessing the direction from the fitted model
+    ## in case a character vector is provided
+
+    preds <- model_fits[[1]]$pred_values
+
+    if (min(preds[-1]) < preds[1]) {
+      ## min non-placebo effect less than placebo
+
+      direction <- "decreasing"
+
+    } else if (max(preds[-1]) > preds[1]) {
+      ## max non-placebo effect greater than placebo
+
+      direction <- "increasing"
+
+    } else {
+
+      direction <- character(0)
+
+      warning(paste0("The direction of the beneficial direction ",
+                     "with increasing dose levels could not be determined."))
+
+    }
+
+  }
+
+  return (direction)
+
+}
+
 addAvgFit <- function (
 
   model_fits,
@@ -470,6 +522,8 @@ getMED <- function (
   checkmate::assert_double(dose_levels, lower = 0, any.missing = FALSE, null.ok = TRUE)
 
   if (!is.null(model_fits)) {
+    ## Direction not needed, because mean value
+    ## (approx. evidence_level = 0.5) is used.
 
     checkmate::assert_class(model_fits, "modelFits")
 
@@ -486,6 +540,13 @@ getMED <- function (
 
   } else {
 
+    if (attr(bs_quantiles, "direction") == "decreasing") {
+
+      bs_quantiles <- bs_quantiles |>
+        dplyr::mutate(q_prob = 1 - q_prob)
+
+    }
+
     if (is.null(dose_levels)) {
 
       dose_levels <- unique(bs_quantiles$dose)
@@ -510,7 +571,8 @@ getMED <- function (
                     sample_type == "diff") |>
       tidyr::pivot_wider(names_from = dose, values_from = q_val) |>
       dplyr::select(-model, -sample_type, -q_prob) |>
-      as.matrix()
+      as.matrix() |>
+      abs()
 
     rownames(abs_diffs) <- unique(bs_quantiles$model)
 

R/posterior.R

@@ -128,7 +128,7 @@ getPosteriorI <- function(
   checkmate::check_vector(se_hat, any.missing = FALSE, null.ok = TRUE)
   checkmate::check_double(se_hat, null.ok = TRUE, lower = 0, upper = Inf)
 
-  if (is.null(mu_hat) && is.null(se_hat)) {
+  if (is.null(mu_hat) && is.null(se_hat) && !is.null(data_i)) {
 
     checkmate::check_data_frame(data_i, null.ok = FALSE)
     checkmate::assert_names(names(data_i), must.include = "response")
@@ -138,7 +138,7 @@ getPosteriorI <- function(
     mu_hat    <- summary(anova_res)$coefficients[, 1]
     se_hat    <- summary(anova_res)$coefficients[, 2]
 
-  } else if (!is.null(mu_hat) && !is.null(se_hat)) {
+  } else if (!is.null(mu_hat) && !is.null(se_hat) && is.null(data_i)) {
 
     stopifnot("m_hat length must match number of dose levels" =
                 length(prior_list) == length(mu_hat))
@@ -148,7 +148,7 @@ getPosteriorI <- function(
 
   } else {
 
-    stop ("Both mu_hat and se_hat must be provided.")
+    stop ("Both mu_hat and se_hat or data_i must be provided.")
 
   }
 
@@ -164,6 +164,8 @@ getPosteriorI <- function(
   names(post_list) <- names(prior_list)
   class(post_list) <- "postList"
 
+  attr(post_list, "ess") <- if (calc_ess) getESS(post_list) else numeric(0)
+
   attr(post_list, "posteriorInfo") <- priorList2priorMix(post_list)
 
   return (post_list)
@@ -287,53 +289,30 @@ postMix2posteriorList <- function (
     SIMPLIFY = FALSE)
 
   # create posterior list
-  posterior_list <- lapply(seq_along(combined_vectors), function (x)
+  post_list <- lapply(seq_along(combined_vectors), function (x)
     do.call(RBesT::mixnorm,
             c(combined_vectors[[x]], sigma = stats::sigma(prior_list[[x]]))))
 
   ## fix component names
-  names(posterior_list) <- names(prior_list)
+  names(post_list) <- names(prior_list)
   comp_names <- lapply(prior_list, colnames)
 
-  for (i in seq_along(posterior_list)) {
+  for (i in seq_along(post_list)) {
 
-    colnames(posterior_list[[i]]) <- comp_names[[i]]
+    colnames(post_list[[i]]) <- comp_names[[i]]
 
   }
   rm(i)
 
   ## set attributes
-  class(posterior_list) <- "postList"
+  class(post_list) <- "postList"
 
-  if (calc_ess) {
-
-    attr(posterior_list, "ess") <- calcEss(calc_ess, posterior_list)
-
-  }
+  attr(post_list, "ess") <- if (calc_ess) getESS(post_list) else numeric(0)
 
   names(post_mix) <- c("weights", "means", "covMats")
-  attr(posterior_list, "posteriorInfo") <- post_mix
-
-  return (posterior_list)
-
-}
+  attr(post_list, "posteriorInfo") <- post_mix
 
-calcEss <- function(calc_ess, posterior_output) {
-
-  checkmate::assert_logical(calc_ess, null.ok = FALSE, len = 1)
-  checkmate::assert_list(posterior_output, names = "named", any.missing = FALSE, null.ok = FALSE)
-
-  if (calc_ess) {
-
-    post_ESS <- getESS(posterior_output)
-
-  } else {
-
-    post_ESS <- numeric(0)
-
-  }
-
-  return(post_ESS)
+  return (post_list)
 
 }
 

cran_submission_script.R

@@ -7,8 +7,8 @@ attachment::att_amend_desc()
 
 # Run tests and examples
 devtools::document()
-devtools::test()
 devtools::run_examples()
+devtools::test()
 # autotest::autotest_package(test = TRUE)
 
 # Check package as CRAN

tests/testthat/setup.R

@@ -340,7 +340,7 @@ getPosteriorI_vec <- function(
   comp_indx    <- createMapping(prior_list)
   comp_mat_ind <- do.call("expand.grid", comp_indx)
 
-  attr(post_list, "ess") <- calcEss(calc_ess, post_list)
+  attr(post_list, "ess") <- if (calc_ess) getESS(post_list) else numeric(0)
 
   diagonals <- lapply(seq_along(comp_mat_ind[, 1]), function(x) {
 

tests/testthat/test-posterior.R

@@ -81,20 +81,21 @@ test_that("getPriorList input parameters do work as intented", {
 })
 
 test_that("getPosteriorI works correctly", {
-  # Prepare test data and parameters
-  data_i <- data.frame(
-    dose = c(0, 1, 2, 3, 4),
-    response = c(10, 20, 30, 40, 50)
-  )
 
   #prior_list <- list(1, 2, 3, 4)
   #mu_hat <- c(10, 20, 30, 40)
   #se_hat <- matrix(c(1, 2, 3, 4), nrow = 4, ncol = 1)
 
   # Test getPosteriorI function
-  post_list <- getPosteriorI(data_i = data_i, prior_list = prior_list_matrix, mu_hat = mu_hat, se_hat = se_hat_vector)
+  post_list <- getPosteriorI(data_i = NULL, prior_list = prior_list_matrix, mu_hat = mu_hat, se_hat = se_hat_vector)
   expect_type(post_list, "list")
   expect_s3_class(post_list, "postList")
+
+  # Prepare test data and parameters
+  data_i <- data.frame(
+    dose = c(0, 1, 2, 3, 4),
+    response = c(10, 20, 30, 40, 50)
+  )
 
   # Test mu_hat and sd_hat both null branch
   post_list <- getPosteriorI(data_i, prior_list_matrix, NULL, NULL)

vignettes/analysis_normal.Rmd

@@ -471,9 +471,10 @@ while the bootstrap median is the median fit of random sampling.
 
 The Minimally Efficacious Dose (MED) per model shape can be assessed with the function `getMED()`.
 ```{r}
-getMED(delta       = 4,
-       model_fits  = model_fits,
-       dose_levels = seq(min(dose_levels), max(dose_levels), by = 0.01))
+getMED(
+  delta       = 4,
+  model_fits  = model_fits,
+  dose_levels = seq(min(dose_levels), max(dose_levels), by = 0.01))
 ```
 For an optional Bayesian decision rule for the MED assessment and further details, please see `?getMED()`.