Update for Analysis Vignette

DESCRIPTION

@@ -38,17 +38,19 @@ Imports:
     RBesT,
     stats
 Suggests:
+    reactable,
+    tibble,
+    quarto,
     clinDR,
     dplyr,
     knitr,
     rmarkdown,
     spelling,
     testthat (>= 3.0.0)
-VignetteBuilder: 
-    knitr
+VignetteBuilder: quarto
 Config/testthat/edition: 3
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 Language: en-US

R/BMCPMod.R

@@ -322,7 +322,7 @@ getCritProb <- function (
 
 #' @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 modelling step
+#' @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
 #' 
 #' @param posterior_list An object of class 'postList' as created by getPosterior() containing information about the (mixture) posterior distribution per dose group
 #' @param contr An object of class 'optContr' as created by the getContr() function. It contains the contrast matrix to be used for the testing step.
@@ -361,7 +361,7 @@ getCritProb <- function (
 #'                       crit_prob_adj  = critVal,
 #'                       simple         = FALSE)
 #' 
-#' @return Bayesian MCP test result as well as modelling result.
+#' @return Bayesian MCP test result as well as modeling result.
 #' 
 #' @export
 performBayesianMCPMod <- function (
@@ -405,7 +405,7 @@ performBayesianMCPMod <- function (
     stop ("Argument 'contr' must be of type 'optContr'")
 
   }
-  
+
   b_mcp <- performBayesianMCP(
     posterior_list = posterior_list,
     contr          = contr,
@@ -414,7 +414,7 @@ performBayesianMCPMod <- function (
   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(

R/bootstrapping.R

@@ -34,13 +34,13 @@
 #' 
 #' @export
 getBootstrapQuantiles <- function (
-
+    
   model_fits,
   quantiles,
   n_samples = 1e3,
   doses     = NULL,
   avg_fit   = TRUE
-
+  
 ) {
 
   mu_hat_samples <- sapply(attr(model_fits, "posterior"),
@@ -105,13 +105,10 @@ getBootstrapQuantiles <- function (
   cr_bounds_data <- cbind(
     doses  = doses,
     models = rep(
-      x    = factor(model_names,
-                    levels = c("linear", "emax", "exponential",
-                               "sigEmax", "logistic", "quadratic",
-                               "avgFit")),
+      x    = factor(model_names, levels = model_names),
       each = length(doses)),
     as.data.frame(quant_mat))
-
+  
   return (cr_bounds_data)
-
+  
 }

R/modelling.R → R/modeling.R

@@ -50,24 +50,26 @@ getModelFits <- function (
   simple = FALSE
 
 ) {
-
+  if (inherits(models, "character")) {
+    models <- stats::setNames(as.list(models), models)
+  }
   checkmate::check_list(models, any.missing = FALSE)
   checkmate::check_double(dose_levels, lower = 0, any.missing = FALSE, len = length(models))
   checkmate::check_class(posterior, "postList")
   checkmate::check_logical(simple)
 
-  models      <- unique(gsub("\\d", "", models))
-  
+  model_names <- unique(gsub("\\d", "", names(models)))
+
   getModelFit <- ifelse(simple, getModelFitSimple, getModelFitOpt)
-  model_fits  <- lapply(models, getModelFit, dose_levels, posterior)
-
+
+  model_fits  <- lapply(model_names, getModelFit, dose_levels, posterior, list("scal" = attr(models, "scal")))
   model_fits  <- addModelWeights(model_fits)
 
-  names(model_fits)             <- models
+  names(model_fits)             <- model_names
   attr(model_fits, "posterior") <- posterior
   class(model_fits)             <- "modelFits"
 
-
+  
   return (model_fits)
 
 }
@@ -77,7 +79,8 @@ getModelFitSimple <- function (
 
   model,
   dose_levels,
-  posterior
+  posterior,
+  addArgs = NULL
 
 ) {
 
@@ -110,45 +113,64 @@ getModelFitOpt <- function (
 
   model,
   dose_levels,
-  posterior
+  posterior,
+  addArgs = NULL
 
 ) {
 
   getOptParams <- function (
 
     model,
     dose_levels,
-    posterior
+    posterior,
+    addArgs = NULL
 
   ) {
-
-    switch (model,
-            "emax" = {
-              lb     <- c(-Inf, -Inf, 0.001 * max(dose_levels))
-              ub     <- c(Inf, Inf, 1.5 * max(dose_levels))
-              expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + (theta[2] * dose_levels) / (theta[3] + dose_levels)))^2 / (post_combs$vars[i, ])))},
-            "sigEmax" = {
-              lb     <- c(-Inf, -Inf, 0.001 * max(dose_levels), 0.5)
-              ub     <- c(Inf, Inf, 1.5 * max(dose_levels), 10)
-              expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + (theta[2] * dose_levels^theta[4]) / (theta[3]^theta[4] + dose_levels^theta[4])))^2 / (post_combs$vars[i, ])))},
-            "exponential" = {
-              lb     <- c(-Inf, -Inf, 0.1 * max(dose_levels))
-              ub     <- c(Inf, Inf, 2 * max(dose_levels))
-              expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] * (exp(dose_levels / theta[3]) - 1)))^2 / (post_combs$vars[i, ])))},
-            "quadratic" = {
-              lb     <- NULL
-              ub     <- c(Inf, Inf, 0)
-              expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] * dose_levels + theta[3] * dose_levels^2))^2 / (post_combs$vars[i, ])))},
-            "linear" = {
-              lb     <- NULL
-              ub     <- NULL
-              expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] * dose_levels))^2 / (post_combs$vars[i, ])))},
-            "logistic" = {
-              lb     <- c(-Inf, -Inf, 0.001 * max(dose_levels), 0.01 * max(dose_levels))
-              ub     <- c(Inf, Inf, 1.5 * max(dose_levels), 0.5 * max(dose_levels))
-              expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] / (1 + exp((theta[3] - dose_levels) / theta[4]))))^2 / (post_combs$vars[i, ])))},
-            {
-              stop ("error")})
+    switch (
+      model,
+      "emax" = {
+        lb     <- c(-Inf, -Inf, 0.001 * max(dose_levels))
+        ub     <- c(Inf, Inf, 1.5 * max(dose_levels))
+        expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + (theta[2] * dose_levels) / (theta[3] + dose_levels)))^2 / (post_combs$vars[i, ])))
+      },
+      "sigEmax" = {
+        lb     <- c(-Inf, -Inf, 0.001 * max(dose_levels), 0.5)
+        ub     <- c(Inf, Inf, 1.5 * max(dose_levels), 10)
+        expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + (theta[2] * dose_levels^theta[4]) / (theta[3]^theta[4] + dose_levels^theta[4])))^2 / (post_combs$vars[i, ])))
+      },
+      "exponential" = {
+        lb     <- c(-Inf, -Inf, 0.1 * max(dose_levels))
+        ub     <- c(Inf, Inf, 2 * max(dose_levels))
+        expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] * (exp(dose_levels / theta[3]) - 1)))^2 / (post_combs$vars[i, ])))
+      },
+      "quadratic" = {
+        lb     <- NULL
+        ub     <- c(Inf, Inf, Inf)
+        expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] * dose_levels + theta[3] * dose_levels^2))^2 / (post_combs$vars[i, ])))
+      },
+      "linear" = {
+        lb     <- NULL
+        ub     <- NULL
+        expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] * dose_levels))^2 / (post_combs$vars[i, ])))
+      },
+      "logistic" = {
+        lb     <- c(-Inf, -Inf, 0.001 * max(dose_levels), 0.01 * max(dose_levels))
+        ub     <- c(Inf, Inf, 1.5 * max(dose_levels), 0.5 * max(dose_levels))
+        expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] / (1 + exp((theta[3] - dose_levels) / theta[4]))))^2 / (post_combs$vars[i, ])))
+      },
+      "betaMod" = {
+        lb     <- c(-Inf, -Inf, 0.05, 0.05)
+        ub     <- c(Inf, Inf, 4, 4)
+        scal   <- ifelse(is.null(addArgs), 1.2 * max(dose_levels), addArgs[["scal"]]) #for betaMod shape
+        expr_i <- substitute(
+          sum(
+              ((post_combs$means[i, ] - (theta[1] + theta[2] * (((theta[3] + theta[4])^(theta[3] + theta[4])) / ((theta[3] ^ theta[3]) * (theta[4]^theta[4]))) * (dose_levels / scal)^(theta[3]) * ((1 - dose_levels / scal)^(theta[4]))))^2 / (post_combs$vars[i, ])))
+          , 
+          list(scal = scal)
+        )
+      },
+      stop(paste("error:", model, "shape not yet implemented"))
+    )
 
     simple_fit <- getModelFitSimple(
       model       = model,
@@ -157,7 +179,7 @@ getModelFitOpt <- function (
 
     param_list <- list(
       expr_i  = expr_i,
-      opts    = list("algorithm" = "NLOPT_LN_NELDERMEAD", maxeval = 1e3),
+      opts    = list("algorithm" = "NLOPT_LN_NELDERMEAD", maxeval = 1e3), #,stopval=0
       params  = list(
         x0 = simple_fit$coeffs,
         lb = lb,
@@ -169,24 +191,21 @@ getModelFitOpt <- function (
   }
 
   optFun <- function (
-
     theta,
-
     dose_levels,
     post_combs,
     expr_i
-
   ) {
 
     comps <- sapply(seq_along(post_combs$weights), function (i) eval(expr_i))
 
     return (sum(post_combs$weights * comps))
 
   }
-  
+
   param_list <- getOptParams(model, dose_levels, posterior)
   post_combs <- getPostCombsI(posterior)
-  
+
   fit <- nloptr::nloptr(
     eval_f = optFun,
     x0     = param_list$params$x0,
@@ -195,7 +214,8 @@ getModelFitOpt <- function (
     opts   = param_list$opts,
     dose_levels = dose_levels,
     post_combs  = post_combs,
-    expr_i      = param_list$expr_i)
+    expr_i      = param_list$expr_i
+  )
 
   names(fit$solution) <- param_list$c_names
 
@@ -204,7 +224,7 @@ getModelFitOpt <- function (
     coeffs      = fit$solution,
     dose_levels = dose_levels)
 
-  model_fit$pred_values <- predictModelFit(model_fit)
+  model_fit$pred_values <- predictModelFit(model_fit = model_fit, doses = model_fit$dose_levels, addArgs = addArgs)
   model_fit$max_effect  <- max(model_fit$pred_values) - min(model_fit$pred_values)
   model_fit$gAIC        <- getGenAIC(model_fit, post_combs)
 
@@ -215,16 +235,13 @@ getModelFitOpt <- function (
 predictModelFit <- function (
 
   model_fit,
-  doses = NULL
+  doses = NULL,
+  addArgs = NULL
 
 ) {
-
   if (is.null(doses)) {
-
     doses <- model_fit$dose_levels
-
   }
-
   pred_vals <- switch (
     model_fit$model,
     "emax" = {DoseFinding::emax(
@@ -258,7 +275,14 @@ predictModelFit <- function (
       model_fit$coeffs["eMax"],
       model_fit$coeffs["ed50"],
       model_fit$coeffs["delta"])},
-    {stop("error")})
+    "betaMod" = {DoseFinding::betaMod(
+      doses,
+      model_fit$coeffs["e0"],
+      model_fit$coeffs["eMax"],
+      model_fit$coeffs["delta1"],
+      model_fit$coeffs["delta2"],
+      ifelse(is.null(addArgs) | is.null(addArgs[["scal"]]), 1.2 * max(doses), addArgs[["scal"]]))},
+    {stop(paste("error:", model_fit$model, "shape not yet implemented"))})
 
   return (pred_vals)
 
@@ -285,7 +309,6 @@ addModelWeights <- function (
   model_fits
 
 ) {
-
   g_aics        <- sapply(model_fits, function (model_fit) model_fit$gAIC)
   exp_values    <- exp(-0.5 * g_aics)
   model_weights <- exp_values / sum(exp_values)

R/plot.R

@@ -87,11 +87,8 @@ plot.modelFits <- function (
   gg_data <- data.frame(
     doses  = rep(dose_seq, length(model_names)),
     fits   = as.vector(preds_models),
-    models = rep(factor(model_names,
-                           levels = c("linear", "emax", "exponential",
-                                      "sigEmax", "logistic", "quadratic",
-                                      "avgFit")),
-                    each = plot_res))
+    models = rep(factor(model_names, levels = model_names),
+                 each = plot_res))
 
   if (gAIC) {
 
@@ -106,7 +103,9 @@ plot.modelFits <- function (
         gsub("quadratic",  "quad", x = _) |>
         gsub("linear",      "lin", x = _) |>
         gsub("logistic",    "log", x = _) |>
-        gsub("sigEmax",    "sigE", x = _)
+        gsub("sigEmax",    "sigE", x = _) |>
+        gsub("betaMod",    "betaM", x = _)|>
+        gsub("quadratic",    "quad", x = _)
 
       label_avg  <- paste0(paste_names, "=", round(mod_weights, 1),
                            collapse = ", ")
@@ -150,7 +149,7 @@ plot.modelFits <- function (
         x = -Inf, y = Inf, hjust = "inward", vjust = "inward",
         size = 3)
 
-      }
+    }
     }
 
   if (cr_bands) {
@@ -195,7 +194,7 @@ plot.modelFits <- function (
         width   = 0,
         alpha   = 0.5)
 
-      }
+    }
     } + 
     ## Posterior Medians
     ggplot2::geom_point(