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-2
 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,19 @@
 # Generated by roxygen2: do not edit by hand
 
-export(BMCPMod)
-export(doFit)
-export(estimateModels)
-export(getGenAICs)
-export(postShape)
+S3method(plot,modelFits)
+S3method(predict,ModelFits)
+S3method(print,BayesianMCP)
+S3method(print,BayesianMCPMod)
+S3method(print,modelFits)
+S3method(print,postList)
+S3method(summary,postList)
+export(assessDesign)
+export(getBootstrapBands)
+export(getContrMat)
+export(getCritProb)
+export(getModelFits)
+export(getPosterior)
+export(getPriorList)
+export(performBayesianMCP)
+export(performBayesianMCPMod)
 export(simulateData)

R/BMCPMod.R

@@ -1,26 +1,277 @@
-#' @title BMCPMod
-#' @param ancova1 tbd
-#'
-#' @param cont_Mat1 tbd
+#' @title assessDesign
+#' 
+#' @param n_patients tbd
+#' @param mods tbd
+#' @param prior_list tbd
+#' @param n_sim tbd
+#' @param alpha_crit_val tbd
+#' @param simple tbd
+#' 
+#' @export
+assessDesign <- function (
+
+  n_patients,
+  mods,
+  prior_list,
+
+  n_sim          = 1e3,
+  alpha_crit_val = 0.05,
+  simple         = TRUE
+
+) {
+
+  dose_levels <- attr(prior_list, "dose_levels")
+
+  data <- simulateData(
+    n_patients  = n_patients,
+    dose_levels = dose_levels,
+    sd          = attr(prior_list, "sd_tot"),
+    mods        = mods,
+    n_sim       = n_sim)
+
+  model_names <- names(mods)
+
+  eval_design <- lapply(model_names, function (model_name) {
+
+    posterior_list <- getPosterior(
+      data       = getModelData(data, model_name),
+      prior_list = prior_list)
+
+    crit_pval <- getCritProb(
+      mods           = mods,
+      dose_levels    = dose_levels,
+      dose_weights   = n_patients,
+      alpha_crit_val = alpha_crit_val)
+
+    contr_mat_prior <- getContrMat(
+      mods           = mods,
+      dose_levels    = dose_levels,
+      dose_weights   = n_patients,
+      prior_list     = prior_list)
+
+    b_mcp_mod <- performBayesianMCPMod(
+      posteriors_list = posterior_list,
+      contr_mat       = contr_mat_prior,
+      crit_prob       = crit_pval,
+      simple          = simple)
+
+  })
+
+  names(eval_design) <- model_names
+
+  return (eval_design)
+
+}
+
+#' @title getContrMat
+#' 
+#' @param mods tbd
+#' @param dose_levels tbd
+#' @param dose_weights tbd
+#' @param prior_list tbd
+#' 
+#' @export
+getContrMat <- function (
+
+  mods,
+  dose_levels,
+  dose_weights,
+  prior_list
+
+) {
+
+  ess_prior <- suppressMessages(round(unlist(lapply(prior_list, RBesT::ess))))
+
+  contr_mat <- DoseFinding::optContr(
+    models = mods,
+    doses  = dose_levels,
+    w      = dose_weights + ess_prior)
+
+  return (contr_mat)
+
+}
+
+#' @title getCritProb
+#' 
+#' @param mods tbd
+#' @param dose_levels tbd
+#' @param dose_weights tbd
+#' @param alpha_crit_val tbd
+#' 
+#' @export
+getCritProb <- function (
+
+  mods,
+  dose_levels,
+  dose_weights,
+  alpha_crit_val
+
+) {
+
+  contr_mat <- DoseFinding::optContr(
+    models = mods,
+    doses  = dose_levels,
+    w      = dose_weights)
+
+  crit_pval <- pnorm(DoseFinding:::critVal(
+    corMat      = contr_mat$corMat,
+    alpha       = alpha_crit_val,
+    df          = 0,
+    alternative = "one.sided"))
+
+  return (crit_pval)
+
+}
+
+#' @title performBayesianMCPMod
+#' 
+#' @param posteriors_list tbd
+#' @param contr_mat tbd
 #' @param crit_prob tbd
-#' @param n_simulations tbd
-#'
+#' @param simple tbd
+#' 
 #' @export
-BMCPMod <- function(
-    ancova1,
-    cont_Mat1,
-    crit_prob,
-    n_simulations) {
-  adj1_p <- list()
+performBayesianMCPMod <- function (
+
+  posteriors_list,
+  contr_mat,
+  crit_prob,
+  simple = FALSE
+
+) {
+
+  if (class(posteriors_list) == "postList") {
+
+    posteriors_list <- list(posteriors_list)
+
+  }
+
+  b_mcp <- performBayesianMCP(
+    posteriors_list = posteriors_list,
+    contr_mat       = contr_mat,
+    crit_prob       = crit_prob)
+
+  model_shapes <- colnames(contr_mat$contMat)
+  dose_levels  <- as.numeric(rownames(contr_mat$contMat))
+
+  fits_list <- lapply(seq_along(posteriors_list), function (i) {
+
+    if (b_mcp[i, 1]) {
+
+      sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "crit_prob")
+
+      model_fits  <- getModelFits(
+        models      = model_shapes,
+        dose_levels = dose_levels,
+        posterior   = posteriors_list[[i]],
+        simple      = simple)
+
+      model_fits  <- addSignificance(model_fits, sign_models)
+
+    } else {
+
+      NULL
+
+    }
+
+  })
+
+  bmcpmod        <- list(BayesianMCP = b_mcp, Mod = fits_list)
+  class(bmcpmod) <- "BayesianMCPMod"
+
+  return (bmcpmod)
+
+}
+
+addSignificance <- function (
+
+  model_fits,
+  sign_models
+
+) {
+
+  names(sign_models) <- NULL
+
+  model_fits_out <- lapply(seq_along(model_fits), function (i) {
+
+    c(model_fits[[i]], significant = sign_models[i])
+
+  })
+
+  attributes(model_fits_out) <- attributes(model_fits)
+
+  return (model_fits_out)
+
+}
 
-  for (i in 1:n_simulations) {
-    ancova <- ancova1[[i]]
-    adj1_p[[i]] <- BayesMCPMod(
-      ancova,
-      cont_Mat1,
-      crit_prob
-    )
+#' @title BayesianMCP
+#' 
+#' @param posteriors_list tbd
+#' @param contr_mat tbd
+#' @param crit_prob tbd
+#' 
+#' @export
+performBayesianMCP <- function(
+
+  posteriors_list,
+  contr_mat,
+  crit_prob
+
+) {
+
+  if (class(posteriors_list) == "postList") {
+
+    posteriors_list <- list(posteriors_list)
+
   }
+
+  b_mcp <- t(sapply(posteriors_list, BayesMCPi, contr_mat, crit_prob))
+
+  attr(b_mcp, "crit_prob") <- crit_prob
+  class(b_mcp)             <- "BayesianMCP"
+
+  return (b_mcp)
+
+}
 
-  return(adj1_p)
+BayesMCPi <- 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)
+
+  }
+
+  post_combs_i <- getPostCombsI(posterior_i)
+  post_probs   <- apply(contr_mat$contMat, 2, getPostProb, post_combs_i)
+
+  res <- c(sign       = ifelse(max(post_probs) > crit_prob, 1, 0),
+           p_val      = max(post_probs),
+           post_probs = post_probs)
+
+  return (res)
+
 }