From a17f3e63a06917d2f159aa27f6cc9caa1dc27445 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Tue, 9 Apr 2024 15:08:40 +0200 Subject: [PATCH] - Fixed a bug that would lead to success rates for the models not printed correctly with print.BayesianMCPMod() in case more than 1 model of the same family was tested, e.g. 2 different shapes of sigEmax --- DESCRIPTION | 2 +- R/BMCPMod.R | 23 ++++++++++++------ R/s3methods.R | 46 ++++++++++------------------------- vignettes/analysis_normal.Rmd | 2 ++ 4 files changed, 32 insertions(+), 41 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index ba222ee..8a05aeb 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: BayesianMCPMod Title: Simulate, Evaluate, and Analyze Dose Finding Trials with Bayesian MCPMod -Version: 1.0.1 +Version: 1.0.2 Authors@R: c( person("Boehringer Ingelheim Pharma GmbH & Co. KG", role = c("cph", "fnd")), person("Stephan", "Wojciekowski", , "stephan.wojciekowski@boehringer-ingelheim.com", role = c("aut", "cre")), diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 3c86b1a..8d1c161 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -384,12 +384,6 @@ performBayesianMCPMod <- function ( } - if (inherits(posterior_list, "postList")) { - - posterior_list <- list(posterior_list) - - } - if (inherits(contr, "optContr")) { model_shapes <- colnames(contr$contMat) @@ -548,11 +542,26 @@ performBayesianMCP <- function( }))) - return (b_mcp) } +getModelSuccesses <- function (b_mcp) { + + stopifnot(inherits(b_mcp, "BayesianMCP")) + + model_indices <- grepl("post_probs.", colnames(b_mcp)) + model_names <- colnames(b_mcp)[model_indices] |> + sub(pattern = "post_probs.", replacement = "", x = _) + + model_successes <- colMeans(b_mcp[, model_indices] > b_mcp[, "crit_prob_adj"]) + + names(model_successes) <- model_names + + return (model_successes) + +} + BayesMCPi <- function ( posterior_i, diff --git a/R/s3methods.R b/R/s3methods.R index 7c763c6..272183c 100644 --- a/R/s3methods.R +++ b/R/s3methods.R @@ -8,40 +8,7 @@ print.BayesianMCPMod <- function ( ) { - model_names <- colnames(x$BayesianMCP)[ - grepl("post_probs.", colnames(x$BayesianMCP))] |> - sub(pattern = "post_probs.", replacement = "", x = _) |> - gsub(pattern = "\\d", replacement = "", x = _) |> - unique(x = _) - - model_success <- colMeans(do.call(rbind, lapply(x$Mod, function (y) { - - if (!is.null(y)) { - - sapply(y, function (z) z$significant) - - } else { - - model_signs <- rep(FALSE, length(model_names)) - names(model_signs) <- model_names - - return (model_signs) - - } - - }))) - print(x$BayesianMCP) - cat("\n") - cat("Model Significance Frequencies\n") - print(model_success, ...) - - if (any(!is.na(attr(x$BayesianMCP, "ess_avg")))) { - - cat("Average Posterior ESS\n") - print(attr(x$BayesianMCP, "ess_avg"), ...) - - } } @@ -57,6 +24,8 @@ print.BayesianMCP <- function ( n_sim <- nrow(x) + model_successes <- getModelSuccesses(x) + cat("Bayesian Multiple Comparison Procedure\n") if (n_sim == 1L) { @@ -74,6 +43,17 @@ print.BayesianMCP <- function ( } + cat("\n") + cat("Model Significance Frequencies\n") + print(model_successes, ...) + + if (any(!is.na(attr(x, "ess_avg")))) { + + cat("Average Posterior ESS\n") + print(attr(x, "ess_avg"), ...) + + } + } ## ModelFits ---------------------------------------------- diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index c39b46b..af6159f 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -252,6 +252,7 @@ It is possible to plot the fitted dose response models and an AIC based average plot(fit_simple) plot(fit) ``` + To assess the uncertainty, one can additionally visualize credible bands (orange shaded areas, default levels are 50% and 95%). These credible bands are calculated with a bootstrap method as follows: Samples from the posterior distribution are drawn and for every sample the simplified fitting step and a prediction is performed. @@ -259,6 +260,7 @@ These predictions are then used to identify and visualize the specified quantile ```{r Plot with bootstrap} plot(fit, cr_bands = TRUE) ``` + The bootstrap based quantiles can also be directly calculated via the getBootstrapQuantiles() function. For this example, only 6 quantiles are bootstrapped for each model fit. ```{r Bootstrap}