Feature/restructure

.github/workflows/test-coverage.yaml

@@ -0,0 +1,38 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/master/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main]
+  pull_request:
+    branches: [main]
+  workflow_dispatch:
+
+name: test-coverage
+
+jobs:
+  test-coverage:
+    runs-on: ubuntu-latest
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      RENV_CONFIG_SANDBOX_ENABLED: false
+
+    steps:
+      - uses: actions/checkout@v2
+
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
+              with:
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: covr
+
+      - name: Install BayesianMCPMod
+        shell: bash
+        run: R CMD INSTALL --preclean .
+
+      - name: Test coverage
+        run: covr::codecov()
+        shell: Rscript {0}

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: BayesianMCPMod
 Title: Bayesian MCPMod
-Version: 0.1.3
+Version: 0.1.3-1
 Authors@R: c(
     person("Sebastian", "Bossert",
            role = "aut",
@@ -22,11 +22,19 @@ Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3
 Imports:
     DoseFinding,
+    ggplot2,
     stats,
     RBesT,
-    nloptr
+    nloptr,
+    clinDR,
+    knitr,
+    rmarkdown,
+    dplyr
 Suggests:
     testthat (>= 3.0.0)
-Config/testthat/edition: 3
 Depends:
   R (>= 4.1)
+VignetteBuilder: knitr
+Config/testthat/edition: 3
+URL: https://github.com/Boehringer-Ingelheim/BayesianMCPMod
+BugReports: https://github.com/Boehringer-Ingelheim/BayesianMCPMod/issues

NAMESPACE

@@ -1,8 +1,8 @@
 # Generated by roxygen2: do not edit by hand
 
 export(BMCPMod)
-export(doFit)
-export(estimateModels)
-export(getGenAICs)
-export(postShape)
+export(getBootsrapBands)
+export(getModelFits)
+export(getPosterior)
+export(plot_modelFits)
 export(simulateData)

R/BMCPMod.R

@@ -1,26 +1,61 @@
 #' @title BMCPMod
-#' @param ancova1 tbd
-#'
-#' @param cont_Mat1 tbd
+#' 
+#' @param posteriors_list tbd
+#' @param contr_mat tbd
 #' @param crit_prob tbd
-#' @param n_simulations tbd
-#'
+#' 
 #' @export
 BMCPMod <- function(
-    ancova1,
-    cont_Mat1,
-    crit_prob,
-    n_simulations) {
-  adj1_p <- list()
+
+  posteriors_list,
+  contr_mat,
+  crit_prob
+
+) {
+
+  lapply(posteriors_list, BayesMCPMod, contr_mat, crit_prob)
+
+}
 
-  for (i in 1:n_simulations) {
-    ancova <- ancova1[[i]]
-    adj1_p[[i]] <- BayesMCPMod(
-      ancova,
-      cont_Mat1,
-      crit_prob
-    )
+BayesMCPMod <- function (
+
+  posterior_i,
+  contr_mat,
+  crit_prob
+
+) {
+
+  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)
+
   }
-
-  return(adj1_p)
+
+  post_combs_i <- getPostCombsI(posterior_i)
+  post_probs   <- apply(contr_mat$contMat, 2, getPostProb, post_combs_i)
+
+  res_list <- list(
+    sign       = ifelse(max(post_probs) > crit_prob, 1, 0),
+    p_val      = max(post_probs),
+    post_probs = post_probs)
+
+  return (res_list)
+
 }

R/bootstrapping.R

@@ -0,0 +1,92 @@
+#' @title getBootsrapBands
+#'
+#' @param model_fits tbd
+#' @param n_samples tbd
+#' @param alpha tbd
+#' @param avg_fit tbd
+#' @param dose_seq tbd
+#'
+#' @return tbd
+#' @export
+getBootsrapBands <- function (
+
+  model_fits,
+  n_samples = 1e3,
+  alpha     = c(0.05, 0.5),
+  avg_fit   = TRUE,
+  dose_seq  = NULL
+
+) {
+
+  mu_hat_samples <- sapply(attr(model_fits, "posterior"),
+                           RBesT::rmix, n = n_samples)
+  sd_hat         <- summary.postList(attr(model_fits, "posterior"))[, 2]
+
+  dose_levels    <- model_fits[[1]]$dose_levels
+  model_names    <- names(model_fits)
+  quantile_probs <- sort(unique(c(alpha / 2, 1 - alpha / 2)))
+
+  if (is.null(dose_seq)) {
+
+    dose_seq <- seq(min(dose_levels), max(dose_levels), length.out = 100L)
+
+  }
+
+  preds <- 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,
+        resp  = mu_hat,
+        S     = diag(sd_hat^2),
+        model = model,
+        type  = "general",
+        bnds  = DoseFinding::defBnds(
+          mD = max(model_fits[[1]]$dose_levels))[[model]])
+
+      preds <- stats::predict(fit, doseSeq = dose_seq, predType = "ls-means")
+      attr(preds, "gAIC") <- DoseFinding::gAIC(fit)
+
+      return (preds)
+
+    }, simplify = FALSE)
+
+    preds_mu_mat <- do.call(rbind, preds_mu_hat)
+
+    if (avg_fit) {
+
+      avg_fit_indx <- which.min(sapply(preds_mu_hat, attr, "gAIC"))
+      preds_mu_mat <- rbind(preds_mu_mat, avgFit = preds_mu_mat[avg_fit_indx, ])
+
+    }
+
+    return (preds_mu_mat)
+
+  })
+
+  if (avg_fit) {
+
+    model_names <- c(model_names, "avgFit")
+
+  }
+
+  sort_indx <- order(rep(seq_along(model_names), length(dose_seq)))
+  quant_mat <- t(apply(X      = preds[sort_indx, ],
+                       MARGIN = 1,
+                       FUN    = stats::quantile,
+                       probs  = quantile_probs))
+
+  cr_bounds_data <- cbind(
+    dose_seqs = dose_seq,
+    models    = rep(
+      x    = factor(model_names,
+                    levels = c("linear", "emax", "exponential",
+                               "sigEmax", "logistic", "quadratic",
+                               "avgFit")),
+      each = length(dose_seq)),
+    as.data.frame(quant_mat))
+
+  return (cr_bounds_data)
+
+}