merge

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: BayesianMCPMod
 Title: Bayesian MCPMod
-Version: 0.1.3-5
+Version: 0.1.3-6
 Authors@R: c(
     person("Sebastian", "Bossert",
            role = "aut",

R/BMCPMod.R

@@ -16,6 +16,8 @@ assessDesign <- function (
   mods,
   prior_list,
 
+  sd             = NULL,
+
   n_sim          = 1e3,
   alpha_crit_val = 0.05,
   simple         = TRUE
@@ -33,11 +35,16 @@ assessDesign <- function (
   # TODO: check that prior_list has 'sd_tot' attribute, and that it's numeric
 
   dose_levels <- attr(prior_list, "dose_levels")
+  sd          <- ifelse(is.null(sd), attr(prior_list, "sd_tot"), sd)
+
+  stopifnot(
+    "sd length must coincide with number of dose levels" =
+      length(sd) == length(dose_levels))
 
   data <- simulateData(
     n_patients  = n_patients,
     dose_levels = dose_levels,
-    sd          = attr(prior_list, "sd_tot"),
+    sd          = sd,
     mods        = mods,
     n_sim       = n_sim)
 
@@ -55,7 +62,7 @@ assessDesign <- function (
       dose_weights   = n_patients,
       alpha_crit_val = alpha_crit_val)
 
-    contr_mat_prior <- getContrMat(
+    contr_mat_prior <- getContr(
       mods           = mods,
       dose_levels    = dose_levels,
       dose_weights   = n_patients,
@@ -75,7 +82,7 @@ assessDesign <- function (
 
 }
 
-#' @title getContrMat
+#' @title getContr
 #' 
 #' @description This function calculates contrast vectors that are optimal for detecting certain alternatives. More information and link to publication will be added.
 #' 
@@ -88,11 +95,12 @@ assessDesign <- function (
 #' 
 #' @export
 getContrMat <- function (
-
   mods,
   dose_levels,
-  dose_weights,
-  prior_list
+  dose_weights = NULL,
+  prior_list   = NULL,
+  se_new_trial = NULL,
+  sd_posterior = NULL
 
 ) {
 
@@ -102,13 +110,55 @@ getContrMat <- function (
   checkmate::check_list(prior_list, names = "named", len = length(attr(prior_list, "dose_levels")), any.missing = FALSE)
 
   ess_prior <- suppressMessages(round(unlist(lapply(prior_list, RBesT::ess))))
+
+  if (is.null(prior_list)) { # frequentist
+
+    if (!is.null(se_new_trial)) { # re-estimate, se_new_trial
+
+      w <- NULL
+      S <- diag((se_new_trial)^2)
+
+    } else { # do not re-estimate, dose_weights
+
+      w <- dose_weights
+      S <- NULL
+
+    }
+
+  } else { # Bayes
+
+    if (!is.null(sd_posterior)) { # re-estimate, sd_posterior
+
+      w <- NULL
+      S <- diag((sd_posterior)^2)
+
+    } else { # do not re-estimate, dose_weights + prior_list
+
+      w <- dose_weights +
+        suppressMessages(round(unlist(lapply(prior_list, RBesT::ess))))
+      S <- NULL
+
+    }
+
+  }
 
-  contr_mat <- DoseFinding::optContr(
-    models = mods,
-    doses  = dose_levels,
-    w      = dose_weights + ess_prior)
+  if (is.null(w)) {
+
+    contr <- DoseFinding::optContr(
+      models = mods,
+      doses  = dose_levels,
+      S      = S)
+
+  } else {
+
+    contr <- DoseFinding::optContr(
+      models = mods,
+      doses  = dose_levels,
+      w      = w)
+
+  }
 
-  return (contr_mat)
+  return (contr)
 
 }
 
@@ -287,6 +337,7 @@ addSignificance <- function (
 
 }
 
+
 #' @title getPostProb
 #' 
 #' @description calculates posterior probabilities for the models.
@@ -344,7 +395,8 @@ BayesMCPi <- function (
 
   res <- c(sign       = ifelse(max(post_probs) > crit_prob, 1, 0),
            p_val      = max(post_probs),
-           post_probs = post_probs)
+           post_probs = post_probs,
+           crit_prob  = crit_prob) # TODO  attr crit_prob??
 
   return (res)
 

R/posterior.R

@@ -10,9 +10,8 @@
 #' @param robustify_weight needs to be provided as a numeric
 #' value for the weight of the robustification component
 #'
-#' @export
 getPriorList <- function (
-    
+
   hist_data,
   dose_levels,
   dose_names       = NULL,
@@ -81,8 +80,9 @@ getPriorList <- function (
 getPosterior <- function(
   data = NULL,
   prior_list,
+  data   = NULL,
   mu_hat = NULL,
-  sd_hat = NULL
+  se_hat = NULL
 
 ) {
   checkmate::check_data_frame(data, null.ok = TRUE)
@@ -101,8 +101,8 @@ getPosterior <- function(
                            prior_list = prior_list,
                            mu_hat     = mu_hat,
                            sd_hat     = sd_hat)
+
   }
-
   if (length(posterior_list) == 1) {
 
     posterior_list <- posterior_list[[1]]
@@ -115,10 +115,10 @@ getPosterior <- function(
 
 getPosteriorI <- function(
 
-  data_i,
+  data_i = NULL,
   prior_list,
   mu_hat = NULL,
-  sd_hat = NULL
+  se_hat = NULL
 
 ) {
 
@@ -136,14 +136,23 @@ getPosteriorI <- function(
     # posterior <- getPosterior(data = simulateData(4, dose_levels, new_trial$sd, mods), prior = prior_list,
     #                           mu_hat = NULL,
     #                           sd_hat = NULL)
-
     anova_res <- stats::lm(data_i$response ~ factor(data_i$dose) - 1)
     mu_hat    <- summary(anova_res)$coefficients[, 1]
-    sd_hat    <- summary(anova_res)$coefficients[, 2]
+    se_hat    <- summary(anova_res)$coefficients[, 2]
+
+  } else if (!is.null(mu_hat) && !is.null(se_hat)) {
+
+    stopifnot("m_hat length must match number of dose levels" = 
+                length(prior_list) == length(mu_hat),
+              "se_hat length must match number of dose levels" = 
+                length(prior_list) == length(se_hat))
+
+  } else {
 
+    stop ("Both mu_hat and se_hat must be provided.")
   }
 
-  post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = sd_hat)
+  post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = se_hat)
 
   if (is.null(names(prior_list))) {
 

R/s3methods.R

@@ -53,6 +53,12 @@ print.BayesianMCP <- function (
   cat("  Estimated Success Rate: ", power, "\n")
   cat("  N Simulations:          ", n_sim)
 
+  ## TODO if n_nim == 1
+  # c(sign       = ifelse(max(post_probs) > crit_prob, 1, 0),
+  #   p_val      = max(post_probs),
+  #   post_probs = post_probs,
+  #   crit_prob  = crit_prob)
+
 }
 
 ## ModelFits ----------------------------------------------

R/simulation.R

@@ -75,4 +75,4 @@ getModelData <- function (
 
   return (model_data)
 
-}
+}

vignettes/Simulation_Example.Rmd

@@ -26,7 +26,7 @@ set.seed(7015)
 # Background and data
 
 In this vignette we will show the use of the Bayesian MCPMod package for trial planning for continuous distributed data. 
-As in [link other vignette] we focus on the indication MDD and make use of historical data that is included in the clinDR package. 
+As in [link to other vignette] we focus on the indication MDD and make use of historical data that is included in the clinDR package. 
 More specifically trial results for BRINTELLIX will be utilized to establish an informative prior for the control group.
 
 # Calculation of a MAP prior