Merge pull request #24 from Boehringer-Ingelheim/avg_model_move

- made avgFit part of the getModelFits function

R/bootstrapping.R

@@ -10,8 +10,7 @@
 #' @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.
+#' @param doses A vector of doses for which a prediction should be performed. If NULL, the dose levels of the model_fits will be used. Default NULL.
 #'
 #' @examples
 #' posterior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2),
@@ -39,8 +38,7 @@ getBootstrapQuantiles <- function (
   model_fits,
   quantiles,
   n_samples = 1e3,
-  doses     = NULL,
-  avg_fit   = TRUE
+  doses     = NULL
 
 ) {
 
@@ -54,23 +52,30 @@ getBootstrapQuantiles <- function (
 
   if (is.null(doses)) {
 
-    doses <- seq(min(dose_levels), max(dose_levels), length.out = 100L)
+    doses <- dose_levels
 
   }
 
+  avg_fit <- "avgFit" %in% model_names
+
+  if (avg_fit) {
+
+    model_names <- setdiff(model_names, "avgFit")
+
+  }
+
   # 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,
+        dose  = dose_levels,
         resp  = mu_hat,
         S     = diag(sd_hat^2),
         model = model,
         type  = "general",
-        bnds  = DoseFinding::defBnds(
-          mD = max(model_fits[[1]]$dose_levels))[[model]])
+        bnds  = DoseFinding::defBnds(mD = max(dose_levels))[[model]])
 
       preds <- stats::predict(fit, doseSeq = doses, predType = "ls-means")
       attr(preds, "gAIC") <- DoseFinding::gAIC(fit)

R/modeling.R

@@ -18,9 +18,10 @@
 #' }
 #' 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 of model names for which a fit will be performed.
+#' @param models List (or vector) of model names for which a fit will be performed.
 #' @param dose_levels A vector containing the different dosage levels.
 #' @param posterior A getPosterior object, containing the (multivariate) posterior distribution per dosage level.
+#' @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 simple Boolean variable, defining whether simplified fit will be applied. Default FALSE.
 #' @examples
 #' posterior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2),
@@ -47,16 +48,19 @@ getModelFits <- function (
   models,
   dose_levels,
   posterior,
-  simple = FALSE
+  avg_fit = TRUE,
+  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(avg_fit)
   checkmate::check_logical(simple)
 
   model_names <- unique(gsub("\\d", "", names(models)))
@@ -65,22 +69,84 @@ getModelFits <- function (
 
   model_fits  <- lapply(model_names, getModelFit, dose_levels, posterior, list("scal" = attr(models, "scal")))
   model_fits  <- addModelWeights(model_fits)
+
+  names(model_fits) <- model_names
+
+  if (avg_fit) {
+
+    model_fits <- addAvgFit(model_fits)
+
+  }
 
-  names(model_fits)             <- model_names
   attr(model_fits, "posterior") <- posterior
   class(model_fits)             <- "modelFits"
 
   return (model_fits)
 
 }
 
+addAvgFit <- function (
+
+  model_fits,
+  doses = NULL
+
+) {
+
+  pred_vals  <- predictAvgFit(model_fits = model_fits, doses = doses)
+
+  avg_fit  <- list(model        = "avgFit",
+                   coeff        = NA,
+                   dose_levels  = model_fits[[1]]$dose_levels,
+                   pred_values  = pred_vals,
+                   max_effect   = max(pred_vals) - min(pred_vals),
+                   gAIC         = NA,
+                   model_weight = NA)
+
+  model_fits <- c(model_fits, avgFit = list(avg_fit))
+
+  return (model_fits)
+
+}
+
+predictAvgFit <- function (
+
+  model_fits,
+  doses = NULL
+
+) {
+
+  model_fits_sub <- model_fits[names(model_fits) != "avgFit"]
+
+  if (is.null(doses)) {
+
+    dose_levels <- model_fits_sub[[1]]$dose_levels
+
+    stopifnot(all(sapply(model_fits_sub,
+                         function (x) identical(x$dose_levels, dose_levels))))
+
+    mod_preds <- sapply(model_fits_sub, function (x) x$pred_values)
+
+  } else {
+
+    mod_preds <- do.call(cbind, predict.modelFits(model_fits_sub, doses = doses))
+
+  }
+
+  mod_weights <- sapply(model_fits_sub, function (x) x$model_weight)
+
+  pred_vals   <- as.vector(mod_preds %*% mod_weights)
+
+  return (pred_vals)
+
+}
+
 ## ignoring mixture posterior
 getModelFitSimple <- function (
 
   model,
   dose_levels,
   posterior,
-  addArgs = NULL
+  add_args = NULL
 
 ) {
 
@@ -114,7 +180,7 @@ getModelFitOpt <- function (
   model,
   dose_levels,
   posterior,
-  addArgs = NULL
+  add_args = NULL
 
 ) {
 
@@ -123,7 +189,7 @@ getModelFitOpt <- function (
     model,
     dose_levels,
     posterior,
-    addArgs = NULL
+    add_args = NULL
 
   ) {
 
@@ -162,7 +228,7 @@ getModelFitOpt <- function (
       "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
+        scal   <- ifelse(is.null(add_args), 1.2 * max(dose_levels), add_args[["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))
       },
@@ -221,24 +287,29 @@ getModelFitOpt <- function (
     coeffs      = fit$solution,
     dose_levels = dose_levels)
 
-  model_fit$pred_values <- predictModelFit(model_fit = model_fit, doses = model_fit$dose_levels, addArgs = addArgs)
+  model_fit$pred_values <- predictModelFit(model_fit = model_fit, doses = model_fit$dose_levels, add_args = add_args)
   model_fit$max_effect  <- max(model_fit$pred_values) - min(model_fit$pred_values)
   model_fit$gAIC        <- getGenAIC(model_fit, post_combs)
 
   return (model_fit)
 
 }
 
+
+
 predictModelFit <- function (
 
   model_fit,
-  doses = NULL,
-  addArgs = NULL
+  doses    = NULL,
+  add_args = NULL
 
 ) {
+
   if (is.null(doses)) {
     doses <- model_fit$dose_levels
   }
+
+  ## cannot predict average values as model_fits required instead of model_fit
   pred_vals <- switch (
     model_fit$model,
     "emax" = {DoseFinding::emax(
@@ -278,7 +349,8 @@ predictModelFit <- function (
       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"]]))},
+      ifelse(is.null(add_args) | is.null(add_args[["scal"]]), 1.2 * max(doses), add_args[["scal"]]))},
+    "avgFit" = {stop("error: use predictAvgFit for the avgFit")},
     {stop(paste("error:", model_fit$model, "shape not yet implemented"))})
 
   return (pred_vals)
@@ -306,6 +378,7 @@ 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

@@ -8,14 +8,14 @@
 #' 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
 #' @param cr_intv Logical value indicating whether credible intervals are included in the plot. Default TRUE
 #' @param alpha_CrI Numerical value of the width of the credible intervals. Default is set to 0.05 (i.e 95% CI are shown).
 #' @param cr_bands Logical value indicating whether bootstrapped based credible bands are shown in the plot. Default FALSE
-#' @param alpha_CrB Numerical vector of the width of the credible bands. Default is set to 0.05 and 0.5 (i.e 95% CB and median are shown).
-#' @param n_bs_smpl Number of bootstrap samples being used. Default set to 1000.
-#' @param acc_color Color of the credible bands. Default set to "orange"
-#' @param ... optional parameter to be passed.
+#' @param alpha_CrB Numerical vector of the width of the credible bands. Default is set to 0.05 and 0.5 (i.e 95% CB and 50% CB  are shown).
+#' @param n_bs_smpl Number of bootstrap samples being used. Default 1000.
+#' @param acc_color Color of the credible bands. Default "orange".
+#' @param plot_res Number of plotted doses within the range of the dose levels, i.e., the resolution of the plot. Default 100.
+#' @param ... optional parameter to be passed to plot().
 #' @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),
@@ -36,30 +36,29 @@ plot.modelFits <- function (
 
   x,
   gAIC      = TRUE,
-  avg_fit   = TRUE,
   cr_intv   = TRUE,
   alpha_CrI = 0.05,
   cr_bands  = FALSE,
   alpha_CrB = c(0.05, 0.5),
   n_bs_smpl = 1e3,
   acc_color = "orange",
+  plot_res  = 1e2,
   ...
 
 ) {
+
   ## R CMD --as-cran appeasement
   .data <- NULL
 
-  checkmate::check_class(x, "modelFits")
   checkmate::check_logical(gAIC)
-  checkmate::check_logical(avg_fit)
   checkmate::check_logical(cr_intv)
   checkmate::check_double(alpha_CrI, lower = 0, upper = 1)
   checkmate::check_logical(cr_bands)
   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)
+  checkmate::check_integer(plot_res, lower = 1, upper = Inf)
 
-  plot_res   <- 1e2
   model_fits <- x
 
   dose_levels  <- model_fits[[1]]$dose_levels
@@ -70,18 +69,10 @@ plot.modelFits <- function (
                   to         = max(dose_levels),
                   length.out = plot_res)
 
-  preds_models <- sapply(model_fits, predictModelFit, doses = dose_seq)
+  preds_models <- do.call(c, predict.modelFits(model_fits, 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")
-
-  }
+
+  avg_fit <- "avgFit" %in% model_names
 
   gg_data <- data.frame(
     doses  = rep(dose_seq, length(model_names)),
@@ -95,6 +86,8 @@ plot.modelFits <- function (
     label_gAUC <- paste("AIC:", round(g_AICs, digits = 1))
 
     if (avg_fit) {
+
+      mod_weights <- sapply(model_fits, function (y) y$model_weight)[-length(model_names)]
 
       mod_weights <- sort(mod_weights, decreasing = TRUE)
       paste_names <- names(mod_weights) |>
@@ -108,7 +101,7 @@ plot.modelFits <- function (
 
       label_avg  <- paste0(paste_names, "=", round(mod_weights, 1),
                            collapse = ", ")
-      label_gAUC <- c(label_gAUC, label_avg)
+      label_gAUC <- c(label_gAUC[-length(label_gAUC)], label_avg)
 
     }
 
@@ -119,14 +112,19 @@ plot.modelFits <- function (
     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,
+      quantiles  = sort(unique(c(0.5, alpha_CrB / 2, 1 - alpha_CrB / 2))),
       doses      = dose_seq)
 
     getInx <- function (alpha_CrB) {
-      n        <- length(alpha_CrB)
-      inx_list <- lapply(seq_len(n), function (i) c(i, 2 * n - i + 1) + 3)
-      return (inx_list)}
+
+      vec_length  <- length(alpha_CrB) * 2 + 1
+      middle_indx <- (vec_length + 1) / 2
+      indx_list   <- lapply(seq_len(middle_indx - 1),
+                            function (i) c(i, vec_length - i + 1) + 2)
+
+      return (indx_list)
+
+    }
 
   }
 

R/posterior.R

@@ -47,11 +47,11 @@ getPosterior <- function(
   checkmate::check_vector(S_hat, any.missing = FALSE, null.ok = TRUE)
   checkmate::check_double(S_hat, null.ok = TRUE, lower = 0, upper = Inf)
 
-  is_matrix_S_hat <- FALSE
-
   stopifnot("prior_list must be an object of RBesT package" =
-              all(sapply(prior_list, function(x) methods::is(x, "normMix") |
-                           methods::is(x, "betaMix") | methods::is(x, "mix"))))
+              all(sapply(prior_list, function(x)
+                methods::is(x, "normMix") |
+                  methods::is(x, "betaMix") |
+                  methods::is(x, "mix"))))
 
   if (!is.null(mu_hat) && !is.null(S_hat) && is.null(data)) {
 
@@ -89,7 +89,7 @@ getPosterior <- function(
         prior_list     = prior_list,
         calc_ess       = calc_ess)
 
-    } else if (!is_matrix_S_hat) {
+    } else {
 
       posterior_list <- getPosteriorI(
         prior_list = prior_list,