Unit tests

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: BayesianMCPMod
 Title: Bayesian MCPMod
-Version: 0.1.3-3
+Version: 0.1.4-3
 Authors@R: c(
     person("Boehringer Ingelheim Pharma GmbH & Co. KG",
            role = c("cph", "fnd")),
@@ -11,8 +11,8 @@ Authors@R: c(
            role = "aut",
            email = "[email protected]"),
     person("Steven", "Brooks",
-           role = "aut",
-           email = "lars.andersen@boehringer-ingelheim.com"),
+           role = "ctb",
+           email = "steven.brooks@boehringer-ingelheim.com"),
     person("Sebastian", "Bossert",
            role = "aut",
            email = "[email protected]")
@@ -32,7 +32,8 @@ Imports:
     clinDR,
     knitr,
     rmarkdown,
-    dplyr
+    dplyr,
+    checkmate
 Suggests:
     testthat (>= 3.0.0)
 Depends:

R/BMCPMod.R

@@ -1,17 +1,17 @@
 #' @title assessDesign
 #'.
 #' @description This function performs simulation based trial design evaluations for a set of specified dose-response models
 #'
 #' @param n_patients Vector specifying the planned number of patients per dose group
 #' @param mods An object of class "Mods" as specified in the DoseFinding package.
 #' @param prior_list a prior_list object specifying the utilized prior for the different dose groups 
 #' @param sd a positive value, specification of assumed sd 
 #' @param n_sim number of simulations to be performed
 #' @param alpha_crit_val (unadjusted) critical value to be used for the MCT testing step. Passed to the getCritProb function for the calculation of adjusted critical values (on the probability scale). Default is 0.05.
 #' @param simple boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.
 #' @param reestimate boolean variable, defining whether critical value should be calculated with re-estimated contrasts (see getCritProb function for more details). Default FALSE
 #' @param contr Allows specification of a fixed contrasts matrix. Default NULL
 #' @param dr_means a vector, allows specification of  individual (not model based) assumed effects per dose group. Default NULL
 #' 
 #' @return returns success probabilities for the different assumed dose-response shapes, attributes also includes information around average success rate (across all assumed models) and prior Effective sample size
 #' 
@@ -24,7 +24,6 @@
 #'                    DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,  
 #'                    DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
 #'                    DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13) ,sigma = 2))
-#' 
 #' n_patients <- c(40,60,60,60,60)
 #' success_probabilities <- assessDesign(
 #' n_patients  = n_patients,
@@ -53,6 +52,15 @@
 
 ) {
 
+  checkmate::assert_vector(n_patients, len = length(attr(mods, "doses")), any.missing = FALSE)
+  checkmate::check_class(mods, classes = "Mods")
+  checkmate::check_list(prior_list, names = "named", len = length(attr(mods, "doses")), any.missing = FALSE)
+  # sensitive to how DoseFinding labels their attributes for "Mods" class
+  checkmate::check_double(n_sim, lower = 1, upper = Inf)
+  checkmate::check_double(alpha_crit_val, lower = 0, upper = 1)
+  checkmate::check_logical(simple)
+  # TODO: check that prior_list has 'sd_tot' attribute, and that it's numeric
+
   dose_levels <- attr(mods, "doses")
 
   data <- simulateData(
@@ -172,6 +180,11 @@
 
 ) {
 
+  checkmate::check_class(mods, classes = "Mods")
+  checkmate::check_double(dose_levels, lower = 0, any.missing = FALSE, len = length(attr(mods, "doses")))
+  checkmate::check_double(dose_weights, any.missing = FALSE, len = length(attr(mods, "doses")))
+  checkmate::check_list(prior_list, names = "named", len = length(attr(mods, "doses")), any.missing = FALSE)
+
   # frequentist & re-estimation
   if (!is.null(se_new_trial) & 
       is.null(dose_weights) & is.null(prior_list) & is.null(sd_posterior)) {
@@ -267,6 +280,11 @@
 
 ) {
 
+  checkmate::check_class(mods, classes = "Mods")
+  checkmate::check_double(dose_levels, lower = 0, any.missing = FALSE, len = length(dose_weights))
+  checkmate::check_double(dose_weights, any.missing = FALSE, len = length(dose_levels))
+  checkmate::check_double(alpha_crit_val, lower = 0, upper = 1)
+
   contr <- getContr(mods           = mods,
                     dose_levels    = dose_levels ,
                     dose_weights   = dose_weights,
@@ -329,6 +347,17 @@
 
 ) {
 
+  checkmate::check_class(posterior_list, "postList")
+  checkmate::check_class(contr, "optContr")
+  checkmate::check_class(crit_prob_adj, "numeric")
+  checkmate::check_logical(simple)
+
+  if (inherits(posterior_list,  "postList")) {
+
+    posterior_list <- list(posterior_list)
+
+  }
+
   if (inherits(posterior_list,  "postList")) {
 
     posterior_list <- list(posterior_list)
@@ -453,6 +482,12 @@
   crit_prob_adj
 
 ) {
+
+  checkmate::check_class(posterior_list, "postList")
+  checkmate::check_class(contr, "optContr")
+  checkmate::check_class(crit_prob_adj, "numeric")
+  checkmate::check_numeric(crit_prob_adj, lower = 0, upper = Inf)
+
   if (inherits(posterior_list,  "postList")) {
 
     posterior_list <- list(posterior_list)

R/modelling.R

@@ -45,6 +45,10 @@ getModelFits <- function (
 
 ) {
 
+  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))
 
   getModelFit <- ifelse(simple, getModelFitSimple, getModelFitOpt)
@@ -137,7 +141,7 @@ getModelFitOpt <- function (
               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 (GENERAL$ERROR$MODEL_OPTIONS)})
+              stop ("error")})
 
     simple_fit <- getModelFitSimple(
       model       = model,
@@ -247,7 +251,7 @@ predictModelFit <- function (
       model_fit$coeffs["eMax"],
       model_fit$coeffs["ed50"],
       model_fit$coeffs["delta"])},
-    {stop(GENERAL$ERROR$MODEL_OPTIONS)})
+    {stop("error")})
 
   return (pred_vals)
 

R/plot.R

@@ -43,6 +43,15 @@ plot.modelFits <- function (
   ...
 
 ) {
+  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)
 
   plot_res   <- 1e2
   model_fits <- x

R/posterior.R

@@ -1,4 +1,3 @@
-
 #' @title getPosterior
 #' 
 #' @description Either the patient level data or both mu_hat as well as sd_hat must to be provided. If patient level data is provided mu_hat and se_hat are calculated within the function using a linear model.
@@ -36,6 +35,12 @@ getPosterior <- function(
   calc_ess = FALSE
 
 ) {
+  checkmate::check_data_frame(data, null.ok = TRUE)
+  checkmate::check_list(prior_list, names = "named", any.missing = FALSE)
+  checkmate::check_vector(mu_hat, any.missing = FALSE, null.ok = TRUE)
+  checkmate::check_double(mu_hat, null.ok = TRUE, lower = -Inf, upper = Inf)
+  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) && is.null(data)) {
 
@@ -55,7 +60,7 @@ getPosterior <- function(
     stop ("Either 'data' or 'mu_hat' and 'se_hat' must not be NULL.")
 
   }
-
+  
   if (length(posterior_list) == 1) {
 
     posterior_list <- posterior_list[[1]]
@@ -75,9 +80,19 @@ getPosteriorI <- function(
   calc_ess = FALSE
 
 ) {
+
+  checkmate::check_data_frame(data_i, null.ok = TRUE)
+  checkmate::check_list(prior_list, names = "named", any.missing = FALSE)
+  checkmate::check_vector(mu_hat, any.missing = FALSE, null.ok = TRUE)
+  checkmate::check_double(mu_hat, null.ok = TRUE, lower = -Inf, upper = Inf)
+  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)) {
-
+    checkmate::check_data_frame(data_i, null.ok = FALSE)
+    checkmate::assert_names(names(data_i), must.include = "response")
+    checkmate::assert_names(names(data_i), must.include = "dose")
+
     anova_res <- stats::lm(data_i$response ~ factor(data_i$dose) - 1)
     mu_hat    <- summary(anova_res)$coefficients[, 1]
     se_hat    <- summary(anova_res)$coefficients[, 2]
@@ -92,7 +107,6 @@ getPosteriorI <- function(
   } else {
 
     stop ("Both mu_hat and se_hat must be provided.")
-
   }
 
   post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = se_hat,

R/simulation.R

@@ -48,6 +48,13 @@ simulateData <- function(
 
 ) {
 
+  checkmate::check_vector(n_patients, any.missing = FALSE, len = length(dose_levels))
+  checkmate::check_double(dose_levels, lower = 0, any.missing = FALSE, len = length(n_patients))
+  checkmate::check_double(sd, len = 1, null.ok = FALSE, lower = 0, upper = Inf)
+  checkmate::check_class(mods, classes = "Mods")
+  checkmate::check_numeric(n_sim, lower = 0, upper = Inf, len = 1)
+  checkmate::check_string(true_model, null.ok = TRUE)
+
   if (!is.null(true_model)) {
 
     n_sim            <- 1
@@ -106,7 +113,7 @@ getModelData <- function (
 
 ) {
 
-  model_data              <- sim_data[, c("simulation", "dose", model_name)]
+  model_data <- sim_data[, c("simulation", "dose", model_name)]
   colnames(model_data)[3] <- "response"
 
   return (model_data)

tests/testthat/setup.R

@@ -0,0 +1,156 @@
+#' @title getPriorList
+#' 
+#' @param hist_data historical trial summary level data,
+#' needs to be provided as a dataframe. Including information of the
+#' estimates and variability.
+#' @param dose_levels vector of the different doseage levels
+#' @param dose_names character vector of dose levels,
+#' default NULL and will be automatically created
+#' based on the dose levels parameter.
+#' @param robust_weight needs to be provided as a numeric
+#' value for the weight of the robustification component
+#'
+getPriorList <- function (
+
+  hist_data,
+  dose_levels,
+  dose_names       = NULL,
+  robust_weight
+
+) {
+
+  checkmate::check_data_frame(hist_data)
+  checkmate::assert_double(dose_levels, lower = 0, any.missing = FALSE)
+  checkmate::check_string(dose_names, null.ok = TRUE)
+  checkmate::check_vector(dose_names, null.ok = TRUE, len = length(dose_levels))
+  checkmate::check_numeric(robust_weight, len = 1, null.ok = FALSE)
+
+  sd_tot <- with(hist_data, sum(sd * n) / sum(n))
+
+  gmap <- RBesT::gMAP(
+    formula    = cbind(est, se) ~ 1 | trial,
+    weights    = hist_data$n,
+    data       = hist_data,
+    family     = gaussian,
+    beta.prior = cbind(0, 100 * sd_tot),
+    tau.dist   = "HalfNormal",
+    tau.prior  = cbind(0, sd_tot / 4))
+
+  prior_ctr <- RBesT::automixfit(gmap)
+
+  prior_ctr <- suppressMessages(RBesT::robustify(
+    priormix = prior_ctr,
+    weight   = robust_weight,
+    sigma    = sd_tot))
+
+
+  prior_trt <- RBesT::mixnorm(
+    comp1 = c(w = 1, m = summary(prior_ctr)[1], n = 1),
+    sigma = sd_tot,
+    param = "mn")
+
+  prior_list <- c(list(prior_ctr),
+                  rep(x     = list(prior_trt),
+                      times = length(dose_levels[-1])))
+
+  if (is.null(dose_names)) {
+    dose_names <- c("Ctr", paste0("DG_", seq_along(dose_levels[-1])))
+  }
+
+  names(prior_list) <- dose_names
+
+  return (prior_list)
+
+}
+
+
+getPostProb <- function (
+
+  contr_j,     # j: dose level
+  post_combs_i # i: simulation outcome
+
+) {
+
+  ## Test statistic = sum over all components of
+  ## posterior weight * normal probability distribution of
+  ## critical values for doses * estimated mean / sqrt(product of critical values for doses)
+
+  ## Calculation for each component of the posterior
+  contr_theta   <- apply(post_combs_i$means, 1, `%*%`, contr_j)
+  contr_var     <- apply(post_combs_i$vars, 1, `%*%`, contr_j^2)
+  contr_weights <- post_combs_i$weights
+
+  ## P(c_m * theta > 0 | Y = y) for a shape m (and dose j)
+  post_probs <- sum(contr_weights * stats::pnorm(contr_theta / sqrt(contr_var)))
+
+  return (post_probs)
+
+}
+
+# Create minimal test case
+n_hist_trials = 2
+
+hist_data <- data.frame(
+  trial = seq(1, n_hist_trials, 1),
+  est   = rep(1, n_hist_trials),
+  se    = rep(1, n_hist_trials),
+  sd    = rep(1, n_hist_trials),
+  n     = rep(1, n_hist_trials)
+)
+
+n_patients <- c(2, 1)
+dose_levels <- c(0, 2.5)
+mean <- c(8, 12)
+sd <- c(0.5, 0.8)
+
+mods <- DoseFinding::Mods(
+  linear = NULL, 
+  doses = dose_levels
+)
+
+
+prior_list <- getPriorList(
+  hist_data   = hist_data,
+  dose_levels = dose_levels,
+  robust_weight = 0.5
+)
+
+n_sim = 1
+alpha_crit_val = 0.05
+simple = TRUE
+
+data <- simulateData(
+  n_patients  = n_patients,
+  dose_levels = dose_levels,
+  sd          = sd,
+  mods        = mods,
+  n_sim       = n_sim
+)
+
+posterior_list <- getPosterior(
+  data = getModelData(data, names(mods)[1]),
+  prior_list = prior_list
+)
+
+contr_mat = getContr(
+  mods = mods, 
+  dose_levels = dose_levels, 
+  dose_weights = n_patients,
+  prior_list = prior_list
+)
+
+crit_pval = getCritProb(
+  mods = mods, 
+  dose_levels = dose_levels, 
+  dose_weights = n_patients, 
+  alpha_crit_val = alpha_crit_val
+)
+
+# eval_design <- assessDesign(
+#   n_patients = n_patients, 
+#   mods = mods, 
+#   prior_list = prior_list,
+#   n_sim = n_sim,
+#   alpha_crit_val = alpha_crit_val,
+#   simple = TRUE
+# )