From a181b11fd646564a960942da19b6a6f2566737da Mon Sep 17 00:00:00 2001 From: Bossert Date: Mon, 4 Dec 2023 16:57:10 +0100 Subject: [PATCH 01/39] Minor update --- vignettes/analysis_normal.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index e84ebac..cdd3535 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -28,7 +28,7 @@ In this vignette we will show the use of the Bayesian MCPMod package for continu Hereby the focus is on the utilization of an informative prior and the BayesianMCPMod evaluation of a single trial. We will use data that is included in the clinDR package. More specifically trial results for BRINTELLIX will be used to illustrate the specification of an informative prior and the usage of such a prior for the bayesian evaluation of a (hypothetical) new trial. -More information around BRINTELLIX to be added... +More information around BRINTELLIX will be added... # Calculation of a MAP prior From f21d98be374d30cef769fda2d02a7e5da017654e Mon Sep 17 00:00:00 2001 From: Bossert Date: Tue, 5 Dec 2023 12:59:57 +0100 Subject: [PATCH 02/39] minor update --- vignettes/Simulation_Example.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/Simulation_Example.Rmd b/vignettes/Simulation_Example.Rmd index bcd8cef..4891721 100644 --- a/vignettes/Simulation_Example.Rmd +++ b/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 From 7a0dac782237d8d8b3c52ae42c7dfac6bc96fc94 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Thu, 7 Dec 2023 10:45:48 +0100 Subject: [PATCH 03/39] - updated getPosterior to handle data = NULL - updated vignette accordingly --- R/posterior.R | 48 ++++++++++++++++++++++------------- vignettes/analysis_normal.Rmd | 11 +++++--- 2 files changed, 38 insertions(+), 21 deletions(-) diff --git a/R/posterior.R b/R/posterior.R index df37f5e..04a88d4 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -5,9 +5,8 @@ #' @param dose_names prior_list #' @param robustify_weight tbd #' -#' @export getPriorList <- function ( - + hist_data, dose_levels, dose_names = NULL, @@ -66,21 +65,34 @@ getPriorList <- function ( #' @param data tbd #' @param prior_list prior_list #' @param mu_hat tbd -#' @param sd_hat tbd +#' @param se_hat tbd #' #' @export getPosterior <- function( - data, + prior_list, + data = NULL, mu_hat = NULL, - sd_hat = NULL + se_hat = NULL ) { - posterior_list <- lapply(split(data, data$simulation), getPosteriorI, - prior_list = prior_list, - mu_hat = mu_hat, - sd_hat = sd_hat) + if (!is.null(mu_hat) && !is.null(se_hat) && is.null(data)) { + + posterior_list <- getPosteriorI( + data_i = NULL, + prior_list = prior_list, + mu_hat = mu_hat, + se_hat = se_hat) + + } else { + + posterior_list <- lapply(split(data, data$simulation), getPosteriorI, + prior_list = prior_list, + mu_hat = mu_hat, + se_hat = se_hat) + + } if (length(posterior_list) == 1) { @@ -94,33 +106,33 @@ getPosterior <- function( getPosteriorI <- function( - data_i, + data_i = NULL, prior_list, mu_hat = NULL, - sd_hat = NULL + se_hat = NULL ) { - if (is.null(mu_hat) && is.null(sd_hat)) { + if (is.null(mu_hat) && is.null(se_hat)) { 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(sd_hat)) { + } 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), - "sd_hat length must match number of dose levels" = - length(prior_list) == length(sd_hat)) + "se_hat length must match number of dose levels" = + length(prior_list) == length(se_hat)) } else { - stop ("Both mu_hat and sd_hat must be provided.") + 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))) { diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index cdd3535..94542e7 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -101,9 +101,14 @@ data_emax <- simulateData( mods = mods, true_model = "emax") -posterior <- getPosterior(prior=prior_list,data=data_emax, - mu_hat = new_trial$rslt, - sd_hat = new_trial$se) +posterior <- getPosterior(prior = prior_list, + data = data_emax) + +posterior <- getPosterior(prior = prior_list, + mu_hat = new_trial$rslt, + se_hat = new_trial$se) + +summary(posterior) ``` From ff9a940b072d43dcf445d6c6944d4b8bcbc07a8b Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Thu, 7 Dec 2023 12:41:17 +0100 Subject: [PATCH 04/39] - updated getContr() --- R/BMCPMod.R | 80 +++++++++++++++++++++++++++++------ R/s3methods.R | 6 +++ R/simulation.R | 2 +- man/assessDesign.Rd | 1 + man/getContr.Rd | 31 ++++++++++++++ man/getContrMat.Rd | 20 --------- man/getPosterior.Rd | 8 ++-- man/performBayesianMCP.Rd | 4 +- vignettes/analysis_normal.Rmd | 2 +- 9 files changed, 112 insertions(+), 42 deletions(-) create mode 100644 man/getContr.Rd delete mode 100644 man/getContrMat.Rd diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 0a3f56c..cf46575 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -14,6 +14,8 @@ assessDesign <- function ( mods, prior_list, + sd = NULL, + n_sim = 1e3, alpha_crit_val = 0.05, simple = TRUE @@ -21,11 +23,16 @@ assessDesign <- function ( ) { 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) @@ -43,7 +50,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, @@ -63,31 +70,75 @@ assessDesign <- function ( } -#' @title getContrMat +#' @title getContr #' #' @param mods tbd #' @param dose_levels tbd #' @param dose_weights tbd #' @param prior_list tbd +#' @param se_new_trial tbd +#' @param sd_posterior tbd #' #' @export -getContrMat <- function ( +getContr <- function ( mods, dose_levels, - dose_weights, - prior_list + dose_weights = NULL, + prior_list = NULL, + se_new_trial = NULL, + sd_posterior = NULL ) { - 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) } @@ -204,7 +255,7 @@ addSignificance <- function ( } -#' @title BayesianMCP +#' @title performBayesianMCP #' #' @param posteriors_list tbd #' @param contr_mat tbd @@ -270,7 +321,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) diff --git a/R/s3methods.R b/R/s3methods.R index 14aa663..a9bb96f 100644 --- a/R/s3methods.R +++ b/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 ---------------------------------------------- diff --git a/R/simulation.R b/R/simulation.R index 7c60082..2bd89b4 100644 --- a/R/simulation.R +++ b/R/simulation.R @@ -61,4 +61,4 @@ getModelData <- function ( return (model_data) -} \ No newline at end of file +} diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd index 6883d23..5233e5b 100644 --- a/man/assessDesign.Rd +++ b/man/assessDesign.Rd @@ -8,6 +8,7 @@ assessDesign( n_patients, mods, prior_list, + sd = NULL, n_sim = 1000, alpha_crit_val = 0.05, simple = TRUE diff --git a/man/getContr.Rd b/man/getContr.Rd new file mode 100644 index 0000000..103d311 --- /dev/null +++ b/man/getContr.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/BMCPMod.R +\name{getContr} +\alias{getContr} +\title{getContr} +\usage{ +getContr( + mods, + dose_levels, + dose_weights = NULL, + prior_list = NULL, + se_new_trial = NULL, + sd_posterior = NULL +) +} +\arguments{ +\item{mods}{tbd} + +\item{dose_levels}{tbd} + +\item{dose_weights}{tbd} + +\item{prior_list}{tbd} + +\item{se_new_trial}{tbd} + +\item{sd_posterior}{tbd} +} +\description{ +getContr +} diff --git a/man/getContrMat.Rd b/man/getContrMat.Rd deleted file mode 100644 index d5393e3..0000000 --- a/man/getContrMat.Rd +++ /dev/null @@ -1,20 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/BMCPMod.R -\name{getContrMat} -\alias{getContrMat} -\title{getContrMat} -\usage{ -getContrMat(mods, dose_levels, dose_weights, prior_list) -} -\arguments{ -\item{mods}{tbd} - -\item{dose_levels}{tbd} - -\item{dose_weights}{tbd} - -\item{prior_list}{tbd} -} -\description{ -getContrMat -} diff --git a/man/getPosterior.Rd b/man/getPosterior.Rd index 501f5de..b404e81 100644 --- a/man/getPosterior.Rd +++ b/man/getPosterior.Rd @@ -4,16 +4,16 @@ \alias{getPosterior} \title{getPosterior} \usage{ -getPosterior(data, prior_list, mu_hat = NULL, sd_hat = NULL) +getPosterior(prior_list, data = NULL, mu_hat = NULL, se_hat = NULL) } \arguments{ -\item{data}{tbd} - \item{prior_list}{prior_list} +\item{data}{tbd} + \item{mu_hat}{tbd} -\item{sd_hat}{tbd} +\item{se_hat}{tbd} } \description{ getPosterior diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd index e6b3543..6493c31 100644 --- a/man/performBayesianMCP.Rd +++ b/man/performBayesianMCP.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/BMCPMod.R \name{performBayesianMCP} \alias{performBayesianMCP} -\title{BayesianMCP} +\title{performBayesianMCP} \usage{ performBayesianMCP(posteriors_list, contr_mat, crit_prob) } @@ -14,5 +14,5 @@ performBayesianMCP(posteriors_list, contr_mat, crit_prob) \item{crit_prob}{tbd} } \description{ -BayesianMCP +performBayesianMCP } diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 94542e7..610fdf0 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -123,7 +123,7 @@ crit_pval <- getCritProb( dose_weights = n_patients, alpha_crit_val = 0.1) -contr_mat_prior <- getContrMat( +contr_mat_prior <- getContr( mods = mods, dose_levels = dose_levels, dose_weights = n_patients, From 629c5ad457e9ad9f981f006c2f3941d777635e11 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Mon, 11 Dec 2023 15:58:40 +0100 Subject: [PATCH 05/39] - removed getPriorList from exported functions - added getPriorList to vignette - fixed s3 method for modelFits predictions - added model weights to s3 print function of modelFits - changed printing behaviour for single simulation BayesianMCP - some renaming of variables & functions - added Steven Brooks as author --- DESCRIPTION | 15 ++- NAMESPACE | 7 +- R/BMCPMod.R | 92 ++++++------- R/bootstrapping.R | 28 ++-- R/plot.R | 26 ++-- R/posterior.R | 73 +---------- R/s3methods.R | 30 +++-- ...strapBands.Rd => getBootstrapQuantiles.Rd} | 10 +- man/performBayesianMCP.Rd | 4 +- man/performBayesianMCPMod.Rd | 9 +- vignettes/analysis_normal.Rmd | 121 ++++++++++++++---- 11 files changed, 222 insertions(+), 193 deletions(-) rename man/{getBootstrapBands.Rd => getBootstrapQuantiles.Rd} (72%) diff --git a/DESCRIPTION b/DESCRIPTION index 8744241..d3a86c4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,19 +2,22 @@ Package: BayesianMCPMod Title: Bayesian MCPMod Version: 0.1.3-2 Authors@R: c( - person("Sebastian", "Bossert", - role = "aut", - email = "sebastian.bossert@boehringer-ingelheim.com"), + person("Boehringer Ingelheim Pharma GmbH & Co. KG", + role = c("cph", "fnd")), person("Stephan", "Wojciekowski", role = c("aut", "cre"), email = "stephan.wojciekowski@boehringer-ingelheim.com"), person("Lars", "Andersen", role = "aut", email = "lars.andersen@boehringer-ingelheim.com"), - person("Boehringer Ingelheim Pharma GmbH & Co. KG", - role = c("cph", "fnd")) + person("Steven", "Brooks", + role = "aut", + email = "lars.andersen@boehringer-ingelheim.com"), + person("Sebastian", "Bossert", + role = "aut", + email = "sebastian.bossert@boehringer-ingelheim.com") ) -Description: Provides functions for the analysis of dose finding trials with Bayesian MCPMod. +Description: Provides functions for the planning and analysis of dose finding trials with Bayesian MCPMod. License: Apache License (>= 2) Encoding: UTF-8 LazyData: true diff --git a/NAMESPACE b/NAMESPACE index f89ff4d..0c8b017 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,19 +1,18 @@ # Generated by roxygen2: do not edit by hand S3method(plot,modelFits) -S3method(predict,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(getBootstrapQuantiles) +export(getContr) export(getCritProb) export(getModelFits) export(getPosterior) -export(getPriorList) export(performBayesianMCP) export(performBayesianMCPMod) export(simulateData) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index cf46575..4d56189 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -44,7 +44,7 @@ assessDesign <- function ( data = getModelData(data, model_name), prior_list = prior_list) - crit_pval <- getCritProb( + crit_prob_adj <- getCritProb( mods = mods, dose_levels = dose_levels, dose_weights = n_patients, @@ -59,7 +59,7 @@ assessDesign <- function ( b_mcp_mod <- performBayesianMCPMod( posteriors_list = posterior_list, contr_mat = contr_mat_prior, - crit_prob = crit_pval, + crit_prob_adj = crit_prob_adj, simple = simple) }) @@ -91,34 +91,40 @@ getContr <- function ( ) { - if (is.null(prior_list)) { # frequentist + # frequentist & re-estimation + if (!is.null(se_new_trial) & + is.null(dose_weights) & is.null(prior_list) & is.null(sd_posterior)) { - 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 - - } + w <- NULL + S <- diag((se_new_trial)^2) - } else { # Bayes + # frequentist & no re-estimation + } else if (!is.null(dose_weights) & + is.null(se_new_trial) & is.null(prior_list) & is.null(sd_posterior)) { - 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 - - } + w <- dose_weights + S <- NULL + + # Bayesian & re-estimation + } else if (!is.null(sd_posterior) & + is.null(se_new_trial) & is.null(prior_list) & is.null(dose_weights)) { + + w <- NULL + S <- diag((sd_posterior)^2) + + # Bayesian & no re-estimation + } else if (!is.null(dose_weights) & !is.null(prior_list) & + is.null(se_new_trial) & is.null(sd_posterior)) { + + w <- dose_weights + + suppressMessages(round(unlist(lapply(prior_list, RBesT::ess)))) + S <- NULL + + } else { + + stop (paste("Provided combiations of 'se_new_trial',", + "'dose_weights', 'prior_list', 'sd_posterior' not allowed.", + "See ?getContr for allowed combinations.")) } @@ -164,13 +170,13 @@ getCritProb <- function ( doses = dose_levels, w = dose_weights) - crit_pval <- pnorm(DoseFinding:::critVal( + crit_prob <- pnorm(DoseFinding:::critVal( corMat = contr_mat$corMat, alpha = alpha_crit_val, df = 0, alternative = "one.sided")) - return (crit_pval) + return (crit_prob) } @@ -178,7 +184,7 @@ getCritProb <- function ( #' #' @param posteriors_list tbd #' @param contr_mat tbd -#' @param crit_prob tbd +#' @param crit_prob_adj tbd #' @param simple tbd #' #' @export @@ -186,7 +192,7 @@ performBayesianMCPMod <- function ( posteriors_list, contr_mat, - crit_prob, + crit_prob_adj, simple = FALSE ) { @@ -200,7 +206,7 @@ performBayesianMCPMod <- function ( b_mcp <- performBayesianMCP( posteriors_list = posteriors_list, contr_mat = contr_mat, - crit_prob = crit_prob) + crit_prob_adj = crit_prob_adj) model_shapes <- colnames(contr_mat$contMat) dose_levels <- as.numeric(rownames(contr_mat$contMat)) @@ -209,7 +215,7 @@ performBayesianMCPMod <- function ( if (b_mcp[i, 1]) { - sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "crit_prob") + sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "crit_prob_adj") model_fits <- getModelFits( models = model_shapes, @@ -259,14 +265,14 @@ addSignificance <- function ( #' #' @param posteriors_list tbd #' @param contr_mat tbd -#' @param crit_prob tbd +#' @param crit_prob_adj tbd #' #' @export performBayesianMCP <- function( posteriors_list, contr_mat, - crit_prob + crit_prob_adj ) { @@ -276,10 +282,10 @@ performBayesianMCP <- function( } - b_mcp <- t(sapply(posteriors_list, BayesMCPi, contr_mat, crit_prob)) + b_mcp <- t(sapply(posteriors_list, BayesMCPi, contr_mat, crit_prob_adj)) - attr(b_mcp, "crit_prob") <- crit_prob - class(b_mcp) <- "BayesianMCP" + attr(b_mcp, "crit_prob_adj") <- crit_prob_adj + class(b_mcp) <- "BayesianMCP" return (b_mcp) @@ -289,7 +295,7 @@ BayesMCPi <- function ( posterior_i, contr_mat, - crit_prob + crit_prob_adj ) { @@ -319,10 +325,10 @@ BayesMCPi <- function ( 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, - crit_prob = crit_prob) # TODO attr crit_prob?? + res <- c(sign = ifelse(max(post_probs) > crit_prob_adj, 1, 0), + crit_prob_adj = crit_prob_adj, + max_post_prob = max(post_probs), + post_probs = post_probs) return (res) diff --git a/R/bootstrapping.R b/R/bootstrapping.R index 8585b83..ec82139 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -1,20 +1,20 @@ -#' @title getBootstrapBands +#' @title getBootstrapQuantiles #' #' @param model_fits tbd #' @param n_samples tbd #' @param alpha tbd #' @param avg_fit tbd -#' @param dose_seq tbd +#' @param doses tbd #' #' @return tbd #' @export -getBootstrapBands <- function ( +getBootstrapQuantiles <- function ( model_fits, + quantiles, n_samples = 1e3, - alpha = c(0.05, 0.5), - avg_fit = TRUE, - dose_seq = NULL + doses = NULL, + avg_fit = TRUE ) { @@ -24,11 +24,11 @@ getBootstrapBands <- function ( dose_levels <- model_fits[[1]]$dose_levels model_names <- names(model_fits) - quantile_probs <- c(0.5, sort(unique(c(alpha / 2, 1 - alpha / 2)))) + quantile_probs <- sort(unique(quantiles)) - if (is.null(dose_seq)) { + if (is.null(doses)) { - dose_seq <- seq(min(dose_levels), max(dose_levels), length.out = 100L) + doses <- seq(min(dose_levels), max(dose_levels), length.out = 100L) } @@ -45,7 +45,7 @@ getBootstrapBands <- function ( bnds = DoseFinding::defBnds( mD = max(model_fits[[1]]$dose_levels))[[model]]) - preds <- stats::predict(fit, doseSeq = dose_seq, predType = "ls-means") + preds <- stats::predict(fit, doseSeq = doses, predType = "ls-means") attr(preds, "gAIC") <- DoseFinding::gAIC(fit) return (preds) @@ -71,20 +71,20 @@ getBootstrapBands <- function ( } - sort_indx <- order(rep(seq_along(model_names), length(dose_seq))) + sort_indx <- order(rep(seq_along(model_names), length(doses))) 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( + doses = doses, + models = rep( x = factor(model_names, levels = c("linear", "emax", "exponential", "sigEmax", "logistic", "quadratic", "avgFit")), - each = length(dose_seq)), + each = length(doses)), as.data.frame(quant_mat)) return (cr_bounds_data) diff --git a/R/plot.R b/R/plot.R index 9f8278b..a18af68 100644 --- a/R/plot.R +++ b/R/plot.R @@ -35,11 +35,11 @@ plot.modelFits <- function ( post_summary <- summary.postList( object = attr(model_fits, "posterior"), probs = c(alpha_CrI / 2, 0.5, 1 - alpha_CrI / 2)) - doses <- seq(from = min(dose_levels), - to = max(dose_levels), - length.out = plot_res) + dose_seq <- seq(from = min(dose_levels), + to = max(dose_levels), + length.out = plot_res) - preds_models <- sapply(model_fits, predictModelFit, doses = doses) + preds_models <- sapply(model_fits, predictModelFit, doses = dose_seq) model_names <- names(model_fits) if (avg_fit) { @@ -53,9 +53,9 @@ plot.modelFits <- function ( } gg_data <- data.frame( - dose_seqs = rep(doses, length(model_names)), - fits = as.vector(preds_models), - models = rep(factor(model_names, + doses = rep(dose_seq, length(model_names)), + fits = as.vector(preds_models), + models = rep(factor(model_names, levels = c("linear", "emax", "exponential", "sigEmax", "logistic", "quadratic", "avgFit")), @@ -86,12 +86,12 @@ plot.modelFits <- function ( if (cr_bands) { - crB_data <- getBootstrapBands( + crB_data <- getBootstrapQuantiles( model_fits = model_fits, n_samples = n_bs_smpl, - alpha = alpha_CrB, + quantiles = c(0.5, sort(unique(c(alpha_CrB / 2, 1 - alpha_CrB / 2)))), avg_fit = avg_fit, - dose_seq = doses) + doses = dose_seq) getInx <- function (alpha_CrB) { n <- length(alpha_CrB) @@ -129,7 +129,7 @@ plot.modelFits <- function ( loop_txt <- paste0( "ggplot2::geom_ribbon( data = crB_data, - mapping = ggplot2::aes(x = dose_seqs, + mapping = ggplot2::aes(x = doses, ymin = crB_data[, ", inx[1], "], ymax = crB_data[, ", inx[2], "]), fill = acc_color, @@ -144,7 +144,7 @@ plot.modelFits <- function ( plts <- plts + ggplot2::geom_line( data = crB_data, - mapping = ggplot2::aes(dose_seqs, `50%`), + mapping = ggplot2::aes(doses, `50%`), color = acc_color) } @@ -175,7 +175,7 @@ plot.modelFits <- function ( ## Fitted Models ggplot2::geom_line( data = gg_data, - mapping = ggplot2::aes(dose_seqs, fits)) + + mapping = ggplot2::aes(doses, fits)) + ## Faceting ggplot2::facet_wrap(~ models) diff --git a/R/posterior.R b/R/posterior.R index 04a88d4..d8b6430 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -1,65 +1,3 @@ -#' @title getPriorList -#' -#' @param hist_data tbd -#' @param dose_levels tbd -#' @param dose_names prior_list -#' @param robustify_weight tbd -#' -getPriorList <- function ( - - hist_data, - dose_levels, - dose_names = NULL, - robustify_weight = 0.5 - -) { - - 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) - - if (!is.null(robustify_weight)) { - - prior_ctr <- suppressMessages(RBesT::robustify( - priormix = prior_ctr, - weight = robustify_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 - - attr(prior_list, "dose_levels") <- dose_levels - attr(prior_list, "sd_tot") <- sd_tot - - return (prior_list) - -} - #' @title getPosterior #' #' @param data tbd @@ -80,17 +18,18 @@ getPosterior <- function( if (!is.null(mu_hat) && !is.null(se_hat) && is.null(data)) { posterior_list <- getPosteriorI( - data_i = NULL, prior_list = prior_list, mu_hat = mu_hat, se_hat = se_hat) - } else { + } else if (is.null(mu_hat) && is.null(se_hat) && !is.null(data)) { posterior_list <- lapply(split(data, data$simulation), getPosteriorI, - prior_list = prior_list, - mu_hat = mu_hat, - se_hat = se_hat) + prior_list = prior_list) + + } else { + + stop ("Either 'data' or 'mu_hat' and 'se_hat' must not be NULL.") } diff --git a/R/s3methods.R b/R/s3methods.R index a9bb96f..efcc99b 100644 --- a/R/s3methods.R +++ b/R/s3methods.R @@ -50,29 +50,38 @@ print.BayesianMCP <- function ( n_sim <- nrow(x) cat("Bayesian Multiple Comparison Procedure\n") - 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) + if (n_sim == 1L) { + + attr(x, "crit_prob_adj") <- NULL + class(x) <- NULL + + print.default(x, ...) + + } else { + + cat(" Estimated Success Rate: ", power, "\n") + cat(" N Simulations: ", n_sim) + + } } ## ModelFits ---------------------------------------------- #' @export -predict.ModelFits <- function ( +predict.modelFits <- function ( object, doses = NULL, ... ) { + + predictions <- lapply(object, predictModelFit, doses = doses) + attr(predictions, "doses") <- doses - lapply(object, predictModelFit, doses = doses, ...) + return (predictions) } @@ -94,7 +103,8 @@ print.modelFits <- function ( out_table <- data.frame(predictions, mEff = sapply(x, function (y) y$max_effect), - gAIC = sapply(x, function (y) y$gAIC)) + gAIC = sapply(x, function (y) y$gAIC), + w = sapply(x, function (y) y$model_weight)) out_table <- apply(as.matrix(out_table), 2, round, digits = n_digits) if (!is.null(x[[1]]$significant)) { diff --git a/man/getBootstrapBands.Rd b/man/getBootstrapQuantiles.Rd similarity index 72% rename from man/getBootstrapBands.Rd rename to man/getBootstrapQuantiles.Rd index de4bd80..e33d842 100644 --- a/man/getBootstrapBands.Rd +++ b/man/getBootstrapQuantiles.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/bootstrapping.R -\name{getBootstrapBands} -\alias{getBootstrapBands} -\title{getBootstrapBands} +\name{getBootstrapQuantiles} +\alias{getBootstrapQuantiles} +\title{getBootstrapQuantiles} \usage{ -getBootstrapBands( +getBootstrapQuantiles( model_fits, n_samples = 1000, alpha = c(0.05, 0.5), @@ -27,5 +27,5 @@ getBootstrapBands( tbd } \description{ -getBootstrapBands +getBootstrapQuantiles } diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd index 6493c31..c556204 100644 --- a/man/performBayesianMCP.Rd +++ b/man/performBayesianMCP.Rd @@ -4,14 +4,14 @@ \alias{performBayesianMCP} \title{performBayesianMCP} \usage{ -performBayesianMCP(posteriors_list, contr_mat, crit_prob) +performBayesianMCP(posteriors_list, contr_mat, crit_prob_adj) } \arguments{ \item{posteriors_list}{tbd} \item{contr_mat}{tbd} -\item{crit_prob}{tbd} +\item{crit_prob_adj}{tbd} } \description{ performBayesianMCP diff --git a/man/performBayesianMCPMod.Rd b/man/performBayesianMCPMod.Rd index 83a4a06..f1929c5 100644 --- a/man/performBayesianMCPMod.Rd +++ b/man/performBayesianMCPMod.Rd @@ -4,14 +4,19 @@ \alias{performBayesianMCPMod} \title{performBayesianMCPMod} \usage{ -performBayesianMCPMod(posteriors_list, contr_mat, crit_prob, simple = FALSE) +performBayesianMCPMod( + posteriors_list, + contr_mat, + crit_prob_adj, + simple = FALSE +) } \arguments{ \item{posteriors_list}{tbd} \item{contr_mat}{tbd} -\item{crit_prob}{tbd} +\item{crit_prob_adj}{tbd} \item{simple}{tbd} } diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 610fdf0..1b5b24c 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -33,12 +33,15 @@ More information around BRINTELLIX will be added... # Calculation of a MAP prior In a first step a meta analytic prior will be calculated using historical data from 4 trials (with main endpoint CfB in MADRS score after 8 weeks). -Please note that only information from the control group will be integrated (leading to an informative multicomponent prior for the control group), while for the active groups a non-informative prior will be specified. -```{r Calculation of a MAP prior} +Please note that only information from the control group will be integrated (leading to an informative multi-component prior for the control group), while for the active groups a non-informative prior will be specified. +```{r Historical Data for Control Arm} data("metaData") -testdata <- as.data.frame(metaData) -dataset <- filter(testdata, bname == "BRINTELLIX") -histcontrol <- filter(dataset, dose == 0, primtime == 8, indication == "MAJOR DEPRESSIVE DISORDER",protid!=6) +dataset <- filter(as.data.frame(metaData), bname == "BRINTELLIX") +histcontrol <- filter(dataset, + dose == 0, + primtime == 8, + indication == "MAJOR DEPRESSIVE DISORDER", + protid != 6) ##Create MAP Prior hist_data <- data.frame( @@ -47,42 +50,101 @@ hist_data <- data.frame( se = histcontrol$se, sd = histcontrol$sd, n = histcontrol$sampsize) - +``` +Here, we suggest a function to construct a list of prior distributions for the different dose groups. +This function is adapted to the needs of this example. +Other applications may need a different way to construct prior distributions. +```{r Defining MAP prior function} + +getPriorList <- function ( + + hist_data, + dose_levels, + dose_names = NULL, + robustify_weight = 0.5 + +) { + + 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) + + if (!is.null(robustify_weight)) { + + prior_ctr <- suppressMessages(RBesT::robustify( + priormix = prior_ctr, + weight = robustify_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) + +} +``` +With the dose levels to be investigated, the prior distribution can be constructed. +```{r Getting the MAP prior} dose_levels <- c(0, 2.5, 5, 10) prior_list <- getPriorList( hist_data = hist_data, dose_levels = dose_levels) - ``` - # Specifications new trial We will use the trial with ct.gov number NCT00635219 as exemplary new trial. -To be able to apply the bayesian MCPMod approach, candidate models need to be specified. Since there are only 3 active dose levels we will limit the set of candidate models to a linear, an exponential and an emax model. +To be able to apply the bayesian MCPMod approach, candidate models need to be specified. +Since there are only 3 active dose levels we will limit the set of candidate models to a linear, an exponential, and an emax model. +Note that the linear candidate model does not require a guesstimate. ```{r Pre-Specification of candidate models} -#Pre-Specification (B)MCPMod - -## candidate models for MCPMod -# linear function - no guestimates needed -exp <- DoseFinding::guesst(d = 5, - p = c(0.2), - model = "exponential", - Maxd = max(dose_levels)) -emax <- DoseFinding::guesst(d = 2.5, - p = c(0.9), - model = "emax") +exp_guesst <- DoseFinding::guesst(d = 5, + p = c(0.2), + model = "exponential", + Maxd = max(dose_levels)) +emax_guesst <- DoseFinding::guesst(d = 2.5, + p = c(0.9), + model = "emax") mods <- DoseFinding::Mods( linear = NULL, - emax = emax, - exponential = exp, + emax = emax_guesst, + exponential = exp_guesst, doses = dose_levels, maxEff = -1, placEff = -12.8) -new_trial <- filter(dataset, primtime == 8, indication == "MAJOR DEPRESSIVE DISORDER",protid==6) +new_trial <- filter(dataset, + primtime == 8, + indication == "MAJOR DEPRESSIVE DISORDER", + protid == 6) + n_patients <- c(150, 142, 147, 149) ``` @@ -91,13 +153,15 @@ n_patients <- c(150, 142, 147, 149) As outlined in citePaper, in a first step the posterior is calculated combining the prior information with the estimated results of the new trial. +The simulation takes the variability of the historical data into account. + ```{r Trial results} #Simulation of new trial ##Note: This step should not be required, as we provide summary measures directly from the new trial data_emax <- simulateData( n_patients = n_patients, dose_levels = dose_levels, - sd = attr(prior_list, "sd_tot"), + sd = with(hist_data, sum(sd * n) / sum(n)), mods = mods, true_model = "emax") @@ -113,8 +177,8 @@ summary(posterior) ``` # Execution of Bayesian MCPMod Test step -For the execution of the testing step of bayesian MCPMod a critical value (on the probability scale) will be determined for a given alpha level. In addition a contrast matrix is created. Please note that there are different possibilities how to generate contrasts. -This information is then used as input for the bayesian MCP testing function. +For the execution of the testing step of Bayesian MCPMod a critical value (on the probability scale) will be determined for a given alpha level. In addition a contrast matrix is created. Please note that there are different possibilities how to generate contrasts. +This information is then used as input for the Bayesian MCP testing function. ```{r Execution of Bayesian MCPMod Test step} crit_pval <- getCritProb( @@ -174,9 +238,12 @@ fit <- getModelFits( posterior = post_observed, simple = FALSE) -plot(fit,cr_bands=TRUE) +plot(fit, cr_bands = TRUE) plot(fit_simple, cr_bands = TRUE) +predict(fit, doses = c(4, 7)) +getBootstrapQuantiles(fit, quantiles = c(0.025, 0.975), doses = c(4, 7)) + ``` # Additional notes TBD, whether certain wrapper functions should be mentioned. From 893ea851f92a841706c4a44bdaccb8eeb2e6ed40 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Mon, 11 Dec 2023 16:04:38 +0100 Subject: [PATCH 06/39] - some Rd file that was overlooked --- man/getCritProb.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/getCritProb.Rd b/man/getCritProb.Rd index 21d3482..41ee06b 100644 --- a/man/getCritProb.Rd +++ b/man/getCritProb.Rd @@ -9,7 +9,7 @@ getCritProb(mods, dose_levels, dose_weights, alpha_crit_val = 0.025) \arguments{ \item{mods}{An object of class "Mods" as specified in the Dosefinding package.} -\item{dose_levels}{vector containing the different doseage levels.} +\item{dose_levels}{vector containing the different dosage levels.} \item{dose_weights}{Vector specifying weights for the different doses} From 7c04fc65ad193e71d6f68763468390054023cc0d Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Mon, 11 Dec 2023 16:35:30 +0100 Subject: [PATCH 07/39] - soved merge in vignette --- vignettes/analysis_normal.Rmd | 5 ----- 1 file changed, 5 deletions(-) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 3a35216..1aaeaf5 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -55,7 +55,6 @@ Here, we suggest a function to construct a list of prior distributions for the d This function is adapted to the needs of this example. Other applications may need a different way to construct prior distributions. ```{r Defining MAP prior function} - getPriorList <- function ( hist_data, @@ -114,12 +113,8 @@ dose_levels <- c(0, 2.5, 5, 10) prior_list <- getPriorList( hist_data = hist_data, -<<<<<<< HEAD dose_levels = dose_levels) -======= - dose_levels = dose_levels,robustify_weight = 0.5) ->>>>>>> 3df7e3809383180120cd61e9cc4f1428f597fd94 ``` # Specifications new trial We will use the trial with ct.gov number NCT00635219 as exemplary new trial. From d62f03784ea0346adea6dec1c71cf05fce697215 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Tue, 12 Dec 2023 23:23:48 +0100 Subject: [PATCH 08/39] - renaming functions - assessDesign: added contr matrix argument, Bayesian re-estimation, improving runtime, output information, average success probability, sd as input - different model shapes of the same type can be evaluated - getPosterior can return posterior ESS - new function getESS to get ess from post_list - simulation function can take a dose-response vector - moved getPrior to vignettes --- NAMESPACE | 3 +- R/BMCPMod.R | 161 ++++++++++++++++++---------- R/modelling.R | 2 + R/plot.R | 4 +- R/posterior.R | 44 ++++++-- R/s3methods.R | 17 ++- R/simulation.R | 25 ++++- man/assessDesign.Rd | 7 +- man/{getContrMat.Rd => getContr.Rd} | 12 +-- man/getESS.Rd | 14 +++ man/getPosterior.Rd | 8 +- man/performBayesianMCP.Rd | 6 +- man/performBayesianMCPMod.Rd | 11 +- man/simulateData.Rd | 3 +- vignettes/Simulation_Example.Rmd | 96 ++++++++++++++--- vignettes/analysis_normal.Rmd | 19 ++-- 16 files changed, 309 insertions(+), 123 deletions(-) rename man/{getContrMat.Rd => getContr.Rd} (72%) create mode 100644 man/getESS.Rd diff --git a/NAMESPACE b/NAMESPACE index 9ea6f28..dc6a5e4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,8 +9,9 @@ S3method(print,postList) S3method(summary,postList) export(assessDesign) export(getBootstrapQuantiles) -export(getContrMat) +export(getContr) export(getCritProb) +export(getESS) export(getModelFits) export(getPosterior) export(performBayesianMCP) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 3cea986..e9e6f98 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -16,29 +16,45 @@ assessDesign <- function ( mods, prior_list, - sd = NULL, + sd, n_sim = 1e3, alpha_crit_val = 0.05, - simple = TRUE + simple = TRUE, + reestimate = FALSE, -) { + contr = NULL, + dr_means = NULL - 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)) + dose_levels <- attr(mods, "doses") data <- simulateData( n_patients = n_patients, dose_levels = dose_levels, sd = sd, mods = mods, - n_sim = n_sim) + n_sim = n_sim, + dr_means = dr_means) + + model_names <- colnames(data)[-c(1:3)] - model_names <- names(mods) + crit_prob_adj <- getCritProb( + mods = mods, + dose_levels = dose_levels, + dose_weights = n_patients, + alpha_crit_val = alpha_crit_val) + + if (!reestimate & is.null(contr)) { + + contr <- getContr( + mods = mods, + dose_levels = dose_levels, + dose_weights = n_patients, + prior_list = prior_list) + + } eval_design <- lapply(model_names, function (model_name) { @@ -46,28 +62,32 @@ assessDesign <- function ( data = getModelData(data, model_name), prior_list = prior_list) - crit_prob_adj <- getCritProb( - mods = mods, - dose_levels = dose_levels, - dose_weights = n_patients, - alpha_crit_val = alpha_crit_val) - - contr_mat_prior <- getContr( - mods = mods, - dose_levels = dose_levels, - dose_weights = n_patients, - prior_list = prior_list) + if (reestimate & is.null(contr)) { + + post_sds <- sapply(posterior_list, function (post) summary(post)[, 2]) + + contr <- apply(post_sds, 2, function (post_sd) getContr( + mods = mods, + dose_levels = dose_levels, + sd_posterior = post_sd)) + + } b_mcp_mod <- performBayesianMCPMod( - posteriors_list = posterior_list, - contr_mat = contr_mat_prior, - crit_prob_adj = crit_prob_adj, - simple = simple) + posterior_list = posterior_list, + contr = contr, + crit_prob_adj = crit_prob_adj, + simple = simple) }) names(eval_design) <- model_names + attr(eval_design, "placEff") <- attr(mods, "placEff") + attr(eval_design, "maxEff") <- attr(mods, "maxEff") + attr(eval_design, "sampleSize") <- n_patients + attr(eval_design, "priorESS") <- getESS(prior_list) + return (eval_design) } @@ -81,16 +101,17 @@ assessDesign <- function ( #' @param dose_weights Vector specifying weights for the different doses #' @param prior_list a prior_list object #' -#' @return contr_mat Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. +#' @return contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. #' #' @export -getContrMat <- function ( +getContr <- function ( + mods, dose_levels, dose_weights = NULL, prior_list = NULL, - se_new_trial = NULL, - sd_posterior = NULL + sd_posterior = NULL, + se_new_trial = NULL ) { @@ -170,13 +191,13 @@ getCritProb <- function ( ) { - contr_mat <- DoseFinding::optContr( + contr <- DoseFinding::optContr( models = mods, doses = dose_levels, w = dose_weights) crit_prob <- pnorm(DoseFinding:::critVal( - corMat = contr_mat$corMat, + corMat = contr$corMat, alpha = alpha_crit_val, df = 0, alternative = "one.sided")) @@ -189,8 +210,8 @@ getCritProb <- function ( #' #' @description performs bayesian MCP Test step and modelling. #' -#' @param posteriors_list a getPosterior object -#' @param contr_mat a getContrMat object, contrast matrix to be used for the testing step. +#' @param posterior_list a getPosterior object +#' @param contr a getContrMat object, contrast matrix to be used for the testing step. #' @param crit_prob a getCritProb object #' @param simple boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE. #' @@ -199,28 +220,41 @@ getCritProb <- function ( #' @export performBayesianMCPMod <- function ( - posteriors_list, - contr_mat, + posterior_list, + contr, crit_prob_adj, simple = FALSE ) { - if (class(posteriors_list) == "postList") { + if (class(posterior_list) == "postList") { - posteriors_list <- list(posteriors_list) + posterior_list <- list(posterior_list) } - b_mcp <- performBayesianMCP( - posteriors_list = posteriors_list, - contr_mat = contr_mat, - crit_prob_adj = crit_prob_adj) + if (inherits(contr, "optContr")) { + + model_shapes <- colnames(contr$contMat) + dose_levels <- as.numeric(rownames(contr$contMat)) + + } else if (length(contr) == length(posterior_list)) { + + model_shapes <- colnames(contr[[1]]$contMat) + dose_levels <- as.numeric(rownames(contr[[1]]$contMat)) + + } else { + + stop ("Argument 'contr' must be of type 'optContr'") + + } - model_shapes <- colnames(contr_mat$contMat) - dose_levels <- as.numeric(rownames(contr_mat$contMat)) + b_mcp <- performBayesianMCP( + posterior_list = posterior_list, + contr = contr, + crit_prob_adj = crit_prob_adj) - fits_list <- lapply(seq_along(posteriors_list), function (i) { + fits_list <- lapply(seq_along(posterior_list), function (i) { if (b_mcp[i, 1]) { @@ -229,7 +263,7 @@ performBayesianMCPMod <- function ( model_fits <- getModelFits( models = model_shapes, dose_levels = dose_levels, - posterior = posteriors_list[[i]], + posterior = posterior_list[[i]], simple = simple) model_fits <- addSignificance(model_fits, sign_models) @@ -274,8 +308,8 @@ addSignificance <- function ( #' #' @description performs bayesian MCP Test step. #' -#' @param posteriors_list a getPosterior object -#' @param contr_mat a getContrMat object, contrast matrix to be used for the testing step. +#' @param posterior_list a getPosterior object +#' @param contr a getContrMat object, contrast matrix to be used for the testing step. #' @param crit_prob a getCritProb object, specifying the critical value to be used for the testing (on the probability scale) #' #' @return b_mcp test result @@ -283,22 +317,39 @@ addSignificance <- function ( #' @export performBayesianMCP <- function( - posteriors_list, - contr_mat, + posterior_list, + contr, crit_prob_adj ) { - if (class(posteriors_list) == "postList") { + if (class(posterior_list) == "postList") { - posteriors_list <- list(posteriors_list) + posterior_list <- list(posterior_list) } - b_mcp <- t(sapply(posteriors_list, BayesMCPi, contr_mat, crit_prob_adj)) + if (inherits(contr, "optContr")) { + + b_mcp <- t(sapply(posterior_list, BayesMCPi, contr, crit_prob_adj)) + + } else { + + b_mcp <- t(mapply(BayesMCPi, posterior_list, contr, crit_prob_adj)) + + } + class(b_mcp) <- "BayesianMCP" attr(b_mcp, "crit_prob_adj") <- crit_prob_adj - class(b_mcp) <- "BayesianMCP" + attr(b_mcp, "ess_avg") <- ifelse( + test = is.na(attr(posterior_list[[1]], "ess")), + yes = numeric(0), + no = rowMeans(sapply(posterior_list, function (posteriors) { + + attr(posteriors, "ess") + + }))) + return (b_mcp) @@ -307,7 +358,7 @@ performBayesianMCP <- function( BayesMCPi <- function ( posterior_i, - contr_mat, + contr, crit_prob_adj ) { @@ -336,7 +387,7 @@ BayesMCPi <- function ( } post_combs_i <- getPostCombsI(posterior_i) - post_probs <- apply(contr_mat$contMat, 2, getPostProb, post_combs_i) + post_probs <- apply(contr$contMat, 2, getPostProb, post_combs_i) res <- c(sign = ifelse(max(post_probs) > crit_prob_adj, 1, 0), crit_prob_adj = crit_prob_adj, diff --git a/R/modelling.R b/R/modelling.R index 1328d24..ed85a16 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -21,6 +21,8 @@ getModelFits <- function ( ) { + models <- unique(gsub("\\d", "", models)) + getModelFit <- ifelse(simple, getModelFitSimple, getModelFitOpt) model_fits <- lapply(models, getModelFit, dose_levels, posterior) diff --git a/R/plot.R b/R/plot.R index d77bac0..26e627a 100644 --- a/R/plot.R +++ b/R/plot.R @@ -47,9 +47,9 @@ plot.modelFits <- function ( if (avg_fit) { mod_weights <- sapply(model_fits, function (x) x$model_weight) - avg_mod <- preds_models %*% mod_weights + avg_preds <- preds_models %*% mod_weights - preds_models <- cbind(preds_models, avg_mod) + preds_models <- cbind(preds_models, avg_preds) model_names <- c(model_names, "avgFit") } diff --git a/R/posterior.R b/R/posterior.R index 11407c6..2823f6f 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -12,9 +12,10 @@ #' @export getPosterior <- function( prior_list, - data = NULL, - mu_hat = NULL, - se_hat = NULL + data = NULL, + mu_hat = NULL, + se_hat = NULL, + calc_ess = FALSE ) { @@ -23,12 +24,13 @@ getPosterior <- function( posterior_list <- getPosteriorI( prior_list = prior_list, mu_hat = mu_hat, - se_hat = se_hat) + se_hat = se_hat, + calc_ess = calc_ess) } else if (is.null(mu_hat) && is.null(se_hat) && !is.null(data)) { posterior_list <- lapply(split(data, data$simulation), getPosteriorI, - prior_list = prior_list) + prior_list = prior_list, calc_ess = calc_ess) } else { @@ -48,10 +50,11 @@ getPosterior <- function( getPosteriorI <- function( - data_i = NULL, + data_i = NULL, prior_list, - mu_hat = NULL, - se_hat = NULL + mu_hat = NULL, + se_hat = NULL, + calc_ess = FALSE ) { @@ -82,13 +85,34 @@ getPosteriorI <- function( } - names(post_list) <- names(prior_list) - class(post_list) <- "postList" + names(post_list) <- names(prior_list) + class(post_list) <- "postList" + attr(post_list, "ess") <- ifelse( + test = calc_ess, + yes = getESS(post_list), + no = numeric(0)) return (post_list) } +#' @title getESS +#' +#' @description blubber +#' +#' @param post_list blubb +#' +#' @export +getESS <- function ( + + post_list + +) { + + suppressMessages(sapply(post_list, RBesT::ess)) + +} + getPostCombsI <- function ( posterior_i diff --git a/R/s3methods.R b/R/s3methods.R index efcc99b..98433fd 100644 --- a/R/s3methods.R +++ b/R/s3methods.R @@ -8,8 +8,8 @@ print.BayesianMCPMod <- function ( ) { - n_models <- ncol(x$BayesianMCP) - 2L - model_names <- colnames(x$BayesianMCP)[-c(1, 2)] |> + model_names <- colnames(x$BayesianMCP)[ + grepl("post_probs.", colnames(x$BayesianMCP))] |> sub(pattern = "post_probs.", replacement = "", x = _) model_success <- colMeans(do.call(rbind, lapply(x$Mod, function (y) { @@ -20,7 +20,7 @@ print.BayesianMCPMod <- function ( } else { - model_signs <- rep(FALSE, n_models) + model_signs <- rep(FALSE, length(model_names)) names(model_signs) <- model_names return (model_signs) @@ -32,7 +32,14 @@ print.BayesianMCPMod <- function ( print(x$BayesianMCP) cat("\n") cat("Model Significance Frequencies\n") - print(model_success, ...) + print(c(avg = mean(model_success), model_success), ...) + + if (!is.na(attr(x$BayesianMCP, "ess_avg"))) { + + cat("Average Posterior ESS\n") + print(attr(x$BayesianMCP, "ess_avg"), ...) + + } } @@ -138,7 +145,7 @@ print.modelFits <- function ( cat("Dose Levels\n", paste(dose_names, round(dose_levels, n_digits), sep = " = "), "\n") cat("\n") - cat("Predictions, Maximum Effect, gAIC & Significance\n") + cat("Predictions, Maximum Effect, gAIC, Model Weights & Significance\n") print(out_table, ...) } diff --git a/R/simulation.R b/R/simulation.R index 6b512c8..15c9149 100644 --- a/R/simulation.R +++ b/R/simulation.R @@ -21,8 +21,9 @@ simulateData <- function( dose_levels, sd, mods, - n_sim = 1e3, - true_model = NULL + n_sim = 1e3, + true_model = NULL, + dr_means = NULL ) { @@ -41,10 +42,24 @@ simulateData <- function( ptno = rep(seq_len(sum(n_patients)), times = n_sim), dose = rep(rep(dose_levels, times = n_patients), times = n_sim)) - model_responses <- DoseFinding::getResp(mods, sim_info$dose) - random_noise <- stats::rnorm(nrow(sim_info), mean = 0, sd = sd) + if (is.null(dr_means)) { + + model_responses <- DoseFinding::getResp(mods, sim_info$dose) + + } else { + + stopifnot(identical(length(dose_levels), length(dr_means))) + + model_responses <- matrix( + data = rep(unlist(mapply(rep, dr_means, n_patients)), + n_sim * length(mods)), + ncol = length(mods), + dimnames = list(NULL, names(mods))) + + } - sim_data <- cbind(sim_info, model_responses + random_noise) + random_noise <- stats::rnorm(nrow(sim_info), mean = 0, sd = sd) + sim_data <- cbind(sim_info, model_responses + random_noise) if (!is.null(true_model)) { diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd index 7817f74..9b0d9a9 100644 --- a/man/assessDesign.Rd +++ b/man/assessDesign.Rd @@ -9,10 +9,13 @@ assessDesign( n_patients, mods, prior_list, - sd = NULL, + sd, n_sim = 1000, alpha_crit_val = 0.05, - simple = TRUE + simple = TRUE, + reestimate = FALSE, + contr = NULL, + dr_means = NULL ) } \arguments{ diff --git a/man/getContrMat.Rd b/man/getContr.Rd similarity index 72% rename from man/getContrMat.Rd rename to man/getContr.Rd index 30d46e9..b1185c6 100644 --- a/man/getContrMat.Rd +++ b/man/getContr.Rd @@ -1,16 +1,16 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/BMCPMod.R -\name{getContrMat} -\alias{getContrMat} +\name{getContr} +\alias{getContr} \title{getContr} \usage{ -getContrMat( +getContr( mods, dose_levels, dose_weights = NULL, prior_list = NULL, - se_new_trial = NULL, - sd_posterior = NULL + sd_posterior = NULL, + se_new_trial = NULL ) } \arguments{ @@ -23,7 +23,7 @@ getContrMat( \item{prior_list}{a prior_list object} } \value{ -contr_mat Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. +contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. } \description{ This function calculates contrast vectors that are optimal for detecting certain alternatives. More information and link to publication will be added. diff --git a/man/getESS.Rd b/man/getESS.Rd new file mode 100644 index 0000000..d3ee3ca --- /dev/null +++ b/man/getESS.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/posterior.R +\name{getESS} +\alias{getESS} +\title{getESS} +\usage{ +getESS(post_list) +} +\arguments{ +\item{post_list}{blubb} +} +\description{ +blubber +} diff --git a/man/getPosterior.Rd b/man/getPosterior.Rd index 209a846..084d5e8 100644 --- a/man/getPosterior.Rd +++ b/man/getPosterior.Rd @@ -4,7 +4,13 @@ \alias{getPosterior} \title{getPosterior} \usage{ -getPosterior(prior_list, data = NULL, mu_hat = NULL, se_hat = NULL) +getPosterior( + prior_list, + data = NULL, + mu_hat = NULL, + se_hat = NULL, + calc_ess = FALSE +) } \arguments{ \item{prior_list}{prior_list object} diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd index 06b4df6..86ad352 100644 --- a/man/performBayesianMCP.Rd +++ b/man/performBayesianMCP.Rd @@ -4,12 +4,12 @@ \alias{performBayesianMCP} \title{performBayesianMCP} \usage{ -performBayesianMCP(posteriors_list, contr_mat, crit_prob_adj) +performBayesianMCP(posterior_list, contr, crit_prob_adj) } \arguments{ -\item{posteriors_list}{a getPosterior object} +\item{posterior_list}{a getPosterior object} -\item{contr_mat}{a getContrMat object, contrast matrix to be used for the testing step.} +\item{contr}{a getContrMat object, contrast matrix to be used for the testing step.} \item{crit_prob}{a getCritProb object, specifying the critical value to be used for the testing (on the probability scale)} } diff --git a/man/performBayesianMCPMod.Rd b/man/performBayesianMCPMod.Rd index 78706b3..d3ce9a5 100644 --- a/man/performBayesianMCPMod.Rd +++ b/man/performBayesianMCPMod.Rd @@ -4,17 +4,12 @@ \alias{performBayesianMCPMod} \title{performBayesianMCPMod} \usage{ -performBayesianMCPMod( - posteriors_list, - contr_mat, - crit_prob_adj, - simple = FALSE -) +performBayesianMCPMod(posterior_list, contr, crit_prob_adj, simple = FALSE) } \arguments{ -\item{posteriors_list}{a getPosterior object} +\item{posterior_list}{a getPosterior object} -\item{contr_mat}{a getContrMat object, contrast matrix to be used for the testing step.} +\item{contr}{a getContrMat object, contrast matrix to be used for the testing step.} \item{simple}{boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.} diff --git a/man/simulateData.Rd b/man/simulateData.Rd index 9afd8da..ae01d13 100644 --- a/man/simulateData.Rd +++ b/man/simulateData.Rd @@ -10,7 +10,8 @@ simulateData( sd, mods, n_sim = 1000, - true_model = NULL + true_model = NULL, + dr_means = NULL ) } \arguments{ diff --git a/vignettes/Simulation_Example.Rmd b/vignettes/Simulation_Example.Rmd index 94960e2..e374f3d 100644 --- a/vignettes/Simulation_Example.Rmd +++ b/vignettes/Simulation_Example.Rmd @@ -34,7 +34,7 @@ In a first step a meta analytic prior will be calculated using historical data f Please note that only information from the control group will be integrated (leading to an informative multicomponent prior for the control group), while for the active groups a non-informative prior will be specified. -```{r Calculation of a MAP prior} +```{r Historical Data} data("metaData") testdata <- as.data.frame(metaData) dataset <- filter(testdata, bname == "BRINTELLIX") @@ -48,12 +48,67 @@ hist_data <- data.frame( sd = histcontrol$sd, n = histcontrol$sampsize) +sd_tot <- with(hist_data, sum(sd * n) / sum(n)) +``` + +```{r Calculation of a MAP prior} +getPriorList <- function ( + + hist_data, + dose_levels, + dose_names = NULL, + robust_weight = 0.5 + +) { + + 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) + + if (!is.null(robust_weight)) { + + 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) + +} + dose_levels <- c(0, 2.5, 5, 10, 20) prior_list <- getPriorList( hist_data = hist_data, - dose_levels = dose_levels, - robustify_weight = 0.5) + dose_levels = dose_levels) ``` @@ -73,22 +128,30 @@ exp <- DoseFinding::guesst(d = 5, emax <- DoseFinding::guesst(d = 2.5, p = c(0.9), model = "emax") -sigemax<- DoseFinding::guesst(d = c(2.5,5), - p = c(0.1,0.6), +sigemax <- DoseFinding::guesst(d = c(2.5, 5), + p = c(0.1, 0.6), + model = "sigEmax") +sigemax2 <- DoseFinding::guesst(d = c(2, 4), + p = c(0.3, 0.8), model = "sigEmax") -#beta <- DoseFinding::guesst(d=5, p=0.8, model="betaMod", dMax=1, scal=1.2, Maxd=20) +# beta <- DoseFinding::guesst(d = 5, +# p = 0.8, +# model = "betaMod", +# dMax = 1, +# scal = 1.2, +# Maxd = 20) mods <- DoseFinding::Mods( linear = NULL, emax = emax, exponential = exp, - sigEmax = sigemax, + sigEmax = rbind(sigemax, sigemax2), #betaMod = beta, doses = dose_levels, maxEff = -3, placEff = -12.8) -n_patients=c(60,80,80,80,80) +n_patients = c(60, 80, 80, 80, 80) ``` # Calculation of success probabilities @@ -96,18 +159,23 @@ n_patients=c(60,80,80,80,80) To calculate success probabilities for the different assumed dose-response models and the specified trial design we will apply the assessDesign function. ```{r} -success_probabilities<-assessDesign (n_patients=n_patients, - mods=mods, - prior_list=prior_list) +success_probabilities <- assessDesign( + n_patients = n_patients, + mods = mods, + prior_list = prior_list, + sd = sd_tot) + success_probabilities ``` As alternative design we will evaluate a design with the same overall sample size but putting more patients on the highest dose group (and control). ```{r} -success_probabilities_uneq<-assessDesign (n_patients=c(80,60,60,60,120), - mods=mods, - prior_list=prior_list) +success_probabilities_uneq <- assessDesign( + n_patients = c(80, 60, 60, 60, 120), + mods = mods, + prior_list = prior_list, + sd = sd_tot) success_probabilities_uneq ``` diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 1aaeaf5..5e735d6 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -59,8 +59,8 @@ getPriorList <- function ( hist_data, dose_levels, - dose_names = NULL, - robustify_weight = 0.5 + dose_names = NULL, + robust_weight = 0.5 ) { @@ -77,11 +77,11 @@ getPriorList <- function ( prior_ctr <- RBesT::automixfit(gmap) - if (!is.null(robustify_weight)) { + if (!is.null(robust_weight)) { prior_ctr <- suppressMessages(RBesT::robustify( priormix = prior_ctr, - weight = robustify_weight, + weight = robust_weight, sigma = sd_tot)) } @@ -114,7 +114,6 @@ dose_levels <- c(0, 2.5, 5, 10) prior_list <- getPriorList( hist_data = hist_data, dose_levels = dose_levels) - ``` # Specifications new trial We will use the trial with ct.gov number NCT00635219 as exemplary new trial. @@ -200,9 +199,11 @@ contr_mat_prior <- getContr( # crit_prob = crit_pval) BMCP_result <- performBayesianMCP( - posteriors_list = posterior, - contr_mat = contr_mat_prior, - crit_prob = crit_pval) + posterior_list = posterior, + contr = contr_mat_prior, + crit_prob_adj = crit_pval) + +BMCP_result #BMCP_result2 <- performBayesianMCPMod( # posteriors_list = posterior_emax, @@ -210,8 +211,6 @@ BMCP_result <- performBayesianMCP( # crit_prob = crit_pval) #BMCP_result2 - -BMCP_result[1,1:5] ``` The testing step is significant indicating a non-flat dose-response shape. In detail the p-values for the emax model is the most significant one. From 1ea587f0f87df5cf03211371a95e1707009cbe54 Mon Sep 17 00:00:00 2001 From: Stephan Wojciekowski Date: Wed, 13 Dec 2023 00:15:08 +0100 Subject: [PATCH 09/39] - addressing some Notes & Warning from devtools::check() --- DESCRIPTION | 2 +- R/BMCPMod.R | 20 +++++++++++++------- R/bootstrapping.R | 4 ++-- R/posterior.R | 7 ++++--- R/simulation.R | 1 + man/assessDesign.Rd | 8 ++++++++ man/getBootstrapQuantiles.Rd | 4 ++-- man/getContr.Rd | 8 ++++++-- man/getPosterior.Rd | 6 ++++-- man/performBayesianMCP.Rd | 2 +- man/performBayesianMCPMod.Rd | 4 ++-- man/simulateData.Rd | 2 ++ 12 files changed, 46 insertions(+), 22 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 313492c..e15756a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -24,7 +24,7 @@ LazyData: true Roxygen: list(markdown = TRUE) RoxygenNote: 7.2.3 Imports: - DoseFinding, + DoseFinding (>= 1.1-1), ggplot2, stats, RBesT, diff --git a/R/BMCPMod.R b/R/BMCPMod.R index e9e6f98..f3f3de6 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -5,9 +5,13 @@ #' @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 tbd #' @param n_sim number of simulations to be performed #' @param alpha_crit_val critical value to be used for the testing (on the probability scale) #' @param simple boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE. +#' @param reestimate tbd Default FALSE +#' @param contr tbd Default NULL +#' @param dr_means tbd Default NULL #' #' @export assessDesign <- function ( @@ -98,8 +102,10 @@ assessDesign <- function ( #' #' @param mods An object of class "Mods" as specified in the Dosefinding package. #' @param dose_levels vector containing the different doseage levels. -#' @param dose_weights Vector specifying weights for the different doses -#' @param prior_list a prior_list object +#' @param dose_weights Vector specifying weights for the different doses. Default NULL +#' @param prior_list a prior_list object. Default NULL +#' @param sd_posterior tbd. Default NULL +#' @param se_new_trial tbd. Default NULL #' #' @return contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. #' @@ -196,7 +202,7 @@ getCritProb <- function ( doses = dose_levels, w = dose_weights) - crit_prob <- pnorm(DoseFinding:::critVal( + crit_prob <- stats::pnorm(DoseFinding::critVal( corMat = contr$corMat, alpha = alpha_crit_val, df = 0, @@ -212,7 +218,7 @@ getCritProb <- function ( #' #' @param posterior_list a getPosterior object #' @param contr a getContrMat object, contrast matrix to be used for the testing step. -#' @param crit_prob a getCritProb object +#' @param crit_prob_adj a getCritProb object #' @param simple boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE. #' #' @return bmcpmod test result as well as modelling result. @@ -227,7 +233,7 @@ performBayesianMCPMod <- function ( ) { - if (class(posterior_list) == "postList") { + if (inherits(posterior_list, "postList")) { posterior_list <- list(posterior_list) @@ -310,7 +316,7 @@ addSignificance <- function ( #' #' @param posterior_list a getPosterior object #' @param contr a getContrMat object, contrast matrix to be used for the testing step. -#' @param crit_prob a getCritProb object, specifying the critical value to be used for the testing (on the probability scale) +#' @param crit_prob_adj a getCritProb object, specifying the critical value to be used for the testing (on the probability scale) #' #' @return b_mcp test result #' @@ -323,7 +329,7 @@ performBayesianMCP <- function( ) { - if (class(posterior_list) == "postList") { + if (inherits(posterior_list, "postList")) { posterior_list <- list(posterior_list) diff --git a/R/bootstrapping.R b/R/bootstrapping.R index ec82139..0957703 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -1,10 +1,10 @@ #' @title getBootstrapQuantiles #' #' @param model_fits tbd +#' @param quantiles tbd #' @param n_samples tbd -#' @param alpha tbd -#' @param avg_fit tbd #' @param doses tbd +#' @param avg_fit tbd #' #' @return tbd #' @export diff --git a/R/posterior.R b/R/posterior.R index 2823f6f..599be98 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -3,11 +3,12 @@ #' #' @description Either the patient level data or both the mu_hat as well as the sd_hat must to be provided. #' -#' @param data dataframe containing the information of dose and response. -#' Also a simulateData object can be provided. #' @param prior_list prior_list object +#' @param data dataframe containing the information of dose and response. Default NULL +#' Also a simulateData object can be provided. #' @param mu_hat vector of estimated mean values -#' @param sd_hat vector of estimated standard deviations. +#' @param se_hat vector of estimated standard deviations. +#' @param calc_ess tbd. Default NULL #' #' @export getPosterior <- function( diff --git a/R/simulation.R b/R/simulation.R index 15c9149..65330cd 100644 --- a/R/simulation.R +++ b/R/simulation.R @@ -10,6 +10,7 @@ #' @param true_model Default value is NULL. #' Assumed true underlying model. Provided via a String. e.g. "emax". #' In case of NULL, all dose-response models, included in the mods input parameter will be used. +#' @param dr_means tbd. Default NULL. #' #' @return sim_data one list object, containing patient level simulated data for all assumed true models. #' Also providing information about simulation iteration, patient number as well as dosage levels. diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd index 9b0d9a9..91e3392 100644 --- a/man/assessDesign.Rd +++ b/man/assessDesign.Rd @@ -25,11 +25,19 @@ assessDesign( \item{prior_list}{a prior_list object specifying the utilized prior for the different dose groups} +\item{sd}{tbd} + \item{n_sim}{number of simulations to be performed} \item{alpha_crit_val}{critical value to be used for the testing (on the probability scale)} \item{simple}{boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.} + +\item{reestimate}{tbd Default FALSE} + +\item{contr}{tbd Default NULL} + +\item{dr_means}{tbd Default NULL} } \description{ This function performs simulation based trial design evaluations for a set of specified dose-response models diff --git a/man/getBootstrapQuantiles.Rd b/man/getBootstrapQuantiles.Rd index 11f88a2..d7f8ecd 100644 --- a/man/getBootstrapQuantiles.Rd +++ b/man/getBootstrapQuantiles.Rd @@ -15,13 +15,13 @@ getBootstrapQuantiles( \arguments{ \item{model_fits}{tbd} +\item{quantiles}{tbd} + \item{n_samples}{tbd} \item{doses}{tbd} \item{avg_fit}{tbd} - -\item{alpha}{tbd} } \value{ tbd diff --git a/man/getContr.Rd b/man/getContr.Rd index b1185c6..0e64c19 100644 --- a/man/getContr.Rd +++ b/man/getContr.Rd @@ -18,9 +18,13 @@ getContr( \item{dose_levels}{vector containing the different doseage levels.} -\item{dose_weights}{Vector specifying weights for the different doses} +\item{dose_weights}{Vector specifying weights for the different doses. Default NULL} -\item{prior_list}{a prior_list object} +\item{prior_list}{a prior_list object. Default NULL} + +\item{sd_posterior}{tbd. Default NULL} + +\item{se_new_trial}{tbd. Default NULL} } \value{ contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. diff --git a/man/getPosterior.Rd b/man/getPosterior.Rd index 084d5e8..44e7560 100644 --- a/man/getPosterior.Rd +++ b/man/getPosterior.Rd @@ -15,12 +15,14 @@ getPosterior( \arguments{ \item{prior_list}{prior_list object} -\item{data}{dataframe containing the information of dose and response. +\item{data}{dataframe containing the information of dose and response. Default NULL Also a simulateData object can be provided.} \item{mu_hat}{vector of estimated mean values} -\item{sd_hat}{vector of estimated standard deviations.} +\item{se_hat}{vector of estimated standard deviations.} + +\item{calc_ess}{tbd. Default NULL} } \description{ Either the patient level data or both the mu_hat as well as the sd_hat must to be provided. diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd index 86ad352..02dd385 100644 --- a/man/performBayesianMCP.Rd +++ b/man/performBayesianMCP.Rd @@ -11,7 +11,7 @@ performBayesianMCP(posterior_list, contr, crit_prob_adj) \item{contr}{a getContrMat object, contrast matrix to be used for the testing step.} -\item{crit_prob}{a getCritProb object, specifying the critical value to be used for the testing (on the probability scale)} +\item{crit_prob_adj}{a getCritProb object, specifying the critical value to be used for the testing (on the probability scale)} } \value{ b_mcp test result diff --git a/man/performBayesianMCPMod.Rd b/man/performBayesianMCPMod.Rd index d3ce9a5..6359109 100644 --- a/man/performBayesianMCPMod.Rd +++ b/man/performBayesianMCPMod.Rd @@ -11,9 +11,9 @@ performBayesianMCPMod(posterior_list, contr, crit_prob_adj, simple = FALSE) \item{contr}{a getContrMat object, contrast matrix to be used for the testing step.} -\item{simple}{boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.} +\item{crit_prob_adj}{a getCritProb object} -\item{crit_prob}{a getCritProb object} +\item{simple}{boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.} } \value{ bmcpmod test result as well as modelling result. diff --git a/man/simulateData.Rd b/man/simulateData.Rd index ae01d13..1b55e08 100644 --- a/man/simulateData.Rd +++ b/man/simulateData.Rd @@ -30,6 +30,8 @@ Default is 1000} \item{true_model}{Default value is NULL. Assumed true underlying model. Provided via a String. e.g. "emax". In case of NULL, all dose-response models, included in the mods input parameter will be used.} + +\item{dr_means}{tbd. Default NULL.} } \value{ sim_data one list object, containing patient level simulated data for all assumed true models. From 25f8f28541ce5be0c55fecee4cdd47a6c8603237 Mon Sep 17 00:00:00 2001 From: wojcieko Date: Thu, 14 Dec 2023 10:32:11 +0100 Subject: [PATCH 10/39] - fixed average success rate in assessDesign and added as attribute - added successRate attribute to BayesianMCP class objects - fixed a warning in print function of BayesianMCPMod that would occur if more than 1 model shape of the same type was specified --- R/BMCPMod.R | 14 ++++++++++---- R/s3methods.R | 10 ++++++---- 2 files changed, 16 insertions(+), 8 deletions(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index f3f3de6..866c1f8 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -85,12 +85,17 @@ assessDesign <- function ( }) + avg_success_rate <- mean(sapply(eval_design, function (bmcpmod) { + attr(bmcpmod$BayesianMCP, "successRate") + })) + names(eval_design) <- model_names - attr(eval_design, "placEff") <- attr(mods, "placEff") - attr(eval_design, "maxEff") <- attr(mods, "maxEff") - attr(eval_design, "sampleSize") <- n_patients - attr(eval_design, "priorESS") <- getESS(prior_list) + attr(eval_design, "avgSuccessRate") <- avg_success_rate + attr(eval_design, "placEff") <- attr(mods, "placEff") + attr(eval_design, "maxEff") <- attr(mods, "maxEff") + attr(eval_design, "sampleSize") <- n_patients + attr(eval_design, "priorESS") <- getESS(prior_list) return (eval_design) @@ -347,6 +352,7 @@ performBayesianMCP <- function( class(b_mcp) <- "BayesianMCP" attr(b_mcp, "crit_prob_adj") <- crit_prob_adj + attr(b_mcp, "successRate") <- mean(b_mcp[, 1]) attr(b_mcp, "ess_avg") <- ifelse( test = is.na(attr(posterior_list[[1]], "ess")), yes = numeric(0), diff --git a/R/s3methods.R b/R/s3methods.R index 98433fd..f4cc5f0 100644 --- a/R/s3methods.R +++ b/R/s3methods.R @@ -10,7 +10,9 @@ print.BayesianMCPMod <- function ( model_names <- colnames(x$BayesianMCP)[ grepl("post_probs.", colnames(x$BayesianMCP))] |> - sub(pattern = "post_probs.", replacement = "", x = _) + sub(pattern = "post_probs.", replacement = "", x = _) |> + gsub(pattern = "\\d", replacement = "", x = _) |> + unique(x = _) model_success <- colMeans(do.call(rbind, lapply(x$Mod, function (y) { @@ -32,7 +34,7 @@ print.BayesianMCPMod <- function ( print(x$BayesianMCP) cat("\n") cat("Model Significance Frequencies\n") - print(c(avg = mean(model_success), model_success), ...) + print(model_success, ...) if (!is.na(attr(x$BayesianMCP, "ess_avg"))) { @@ -53,7 +55,6 @@ print.BayesianMCP <- function ( ) { - power <- mean(x[, 1]) n_sim <- nrow(x) cat("Bayesian Multiple Comparison Procedure\n") @@ -61,13 +62,14 @@ print.BayesianMCP <- function ( if (n_sim == 1L) { attr(x, "crit_prob_adj") <- NULL + attr(x, "success_rate") <- NULL class(x) <- NULL print.default(x, ...) } else { - cat(" Estimated Success Rate: ", power, "\n") + cat(" Estimated Success Rate: ", attr(x, "successRate"), "\n") cat(" N Simulations: ", n_sim) } From 7ff4008ec5e9cdcc2f703b1973e074cd61d236db Mon Sep 17 00:00:00 2001 From: Bossert Date: Fri, 15 Dec 2023 10:15:10 +0000 Subject: [PATCH 11/39] Minor update in documentation of getContr --- R/BMCPMod.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 866c1f8..b21ca57 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -109,7 +109,7 @@ assessDesign <- function ( #' @param dose_levels vector containing the different doseage levels. #' @param dose_weights Vector specifying weights for the different doses. Default NULL #' @param prior_list a prior_list object. Default NULL -#' @param sd_posterior tbd. Default NULL +#' @param sd_posterior a vector of positive numerics. Default NULL #' @param se_new_trial tbd. Default NULL #' #' @return contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. From d2ae16569a25e3b418c8f0bf1a497705049955ac Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Fri, 15 Dec 2023 10:27:01 +0000 Subject: [PATCH 12/39] additional documentation getContr --- R/BMCPMod.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index b21ca57..76f2529 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -110,7 +110,7 @@ assessDesign <- function ( #' @param dose_weights Vector specifying weights for the different doses. Default NULL #' @param prior_list a prior_list object. Default NULL #' @param sd_posterior a vector of positive numerics. Default NULL -#' @param se_new_trial tbd. Default NULL +#' @param se_new_trial a vector of positive numerics. Default NULL #' #' @return contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. #' From 24ef42d34176757721d1f82d16b1072fe6f9d3c1 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Fri, 15 Dec 2023 12:16:52 +0000 Subject: [PATCH 13/39] Addition of performBayesianMCPMod function within vignette --- vignettes/analysis_normal.Rmd | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 5e735d6..70a70af 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -246,6 +246,13 @@ getBootstrapQuantiles(fit, quantiles = c(0.025, 0.975), doses = c(4, 7)) ``` # Additional notes -TBD, whether certain wrapper functions should be mentioned. +It is also possible to perform the testing and modelling step in combined fashion via the performBayesianMCPMod function +```{r} +performBayesianMCPMod( + posterior_list = posterior, + contr = contr_mat_prior, + crit_prob_adj = crit_pval) +``` + From a289dff78e6257a29a8c9e15a4cac42652b2f97d Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Fri, 15 Dec 2023 12:53:45 +0000 Subject: [PATCH 14/39] First updates around documentation --- R/BMCPMod.R | 21 +++++++++++++-------- R/posterior.R | 16 +++++++++------- man/getContr.Rd | 4 ++-- man/getCritProb.Rd | 10 ++++++++-- man/getESS.Rd | 7 +++++-- man/getPosterior.Rd | 12 ++++++++---- man/performBayesianMCP.Rd | 7 ++++--- 7 files changed, 49 insertions(+), 28 deletions(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 76f2529..446ecc7 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -185,6 +185,8 @@ getContr <- function ( #' @title getCritProb #' +#' @description This function calculates multiplicity adjusted +#' #' @param mods An object of class "Mods" as specified in the Dosefinding package. #' @param dose_levels vector containing the different dosage levels. #' @param dose_weights Vector specifying weights for the different doses @@ -197,15 +199,17 @@ getCritProb <- function ( mods, dose_levels, - dose_weights, + dose_weights =NULL, + se_new_trial = NULL, alpha_crit_val = 0.025 ) { - contr <- DoseFinding::optContr( - models = mods, - doses = dose_levels, - w = dose_weights) + contr <- getContr(mods = mods, + dose_levels = dose_levels , + dose_weights = dose_weights, + se_new_trial = se_new_trial, + alpha_crit_val = alpha_crit_val) crit_prob <- stats::pnorm(DoseFinding::critVal( corMat = contr$corMat, @@ -317,13 +321,14 @@ addSignificance <- function ( #' @title performBayesianMCP #' -#' @description performs bayesian MCP Test step. +#' @description performs bayesian MCP Test step, as described in Fleischer et al. (Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670.) +#' Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate is calculated that the contrast is bigger than 0 given the data observed #' -#' @param posterior_list a getPosterior object +#' @param posterior_list a getPosterior object with information about the (mixture) posterior distribution per dose group #' @param contr a getContrMat object, contrast matrix to be used for the testing step. #' @param crit_prob_adj a getCritProb object, specifying the critical value to be used for the testing (on the probability scale) #' -#' @return b_mcp test result +#' @return b_mcp test result, with information about p-values for the individual dose-response shapes #' #' @export performBayesianMCP <- function( diff --git a/R/posterior.R b/R/posterior.R index 599be98..c29e011 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -1,15 +1,16 @@ #' @title getPosterior #' -#' @description Either the patient level data or both the mu_hat as well as the sd_hat must to be provided. +#' @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. +#' This function calculates the posterior for every dose group independently via the RBesT function postmix. #' #' @param prior_list prior_list object #' @param data dataframe containing the information of dose and response. Default NULL #' Also a simulateData object can be provided. -#' @param mu_hat vector of estimated mean values -#' @param se_hat vector of estimated standard deviations. -#' @param calc_ess tbd. Default NULL -#' +#' @param mu_hat vector of estimated mean values (per dose group). +#' @param se_hat vector of estimated standard deviations (per dose group). +#' @param calc_ess boolean variable, indicating whether effective sample size should be calculated. Default FALSE +#' @return posterior_list, a posterior list object is returned with information about (mixture) posterior distribution per dose group #' @export getPosterior <- function( prior_list, @@ -99,10 +100,11 @@ getPosteriorI <- function( #' @title getESS #' -#' @description blubber +#' @description This function calculates the effective sample size for every dose group via the RBesT function ess. #' -#' @param post_list blubb +#' @param post_list a posterior list object, for which the effective sample size (per dose group) should be calculated #' +#' @return a vector of the effective sample sizes (per dose group) #' @export getESS <- function ( diff --git a/man/getContr.Rd b/man/getContr.Rd index 0e64c19..809e932 100644 --- a/man/getContr.Rd +++ b/man/getContr.Rd @@ -22,9 +22,9 @@ getContr( \item{prior_list}{a prior_list object. Default NULL} -\item{sd_posterior}{tbd. Default NULL} +\item{sd_posterior}{a vector of positive numerics. Default NULL} -\item{se_new_trial}{tbd. Default NULL} +\item{se_new_trial}{a vector of positive numerics. Default NULL} } \value{ contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. diff --git a/man/getCritProb.Rd b/man/getCritProb.Rd index 41ee06b..fa648bd 100644 --- a/man/getCritProb.Rd +++ b/man/getCritProb.Rd @@ -4,7 +4,13 @@ \alias{getCritProb} \title{getCritProb} \usage{ -getCritProb(mods, dose_levels, dose_weights, alpha_crit_val = 0.025) +getCritProb( + mods, + dose_levels, + dose_weights = NULL, + se_new_trial = NULL, + alpha_crit_val = 0.025 +) } \arguments{ \item{mods}{An object of class "Mods" as specified in the Dosefinding package.} @@ -19,5 +25,5 @@ getCritProb(mods, dose_levels, dose_weights, alpha_crit_val = 0.025) crit_pval multiplicity adjusted critical value on the probability scale. } \description{ -getCritProb +This function calculates multiplicity adjusted } diff --git a/man/getESS.Rd b/man/getESS.Rd index d3ee3ca..7f08776 100644 --- a/man/getESS.Rd +++ b/man/getESS.Rd @@ -7,8 +7,11 @@ getESS(post_list) } \arguments{ -\item{post_list}{blubb} +\item{post_list}{a posterior list object, for which the effective sample size (per dose group) should be calculated} +} +\value{ +a vector of the effective sample sizes (per dose group) } \description{ -blubber +This function calculates the effective sample size for every dose group via the RBesT function ess. } diff --git a/man/getPosterior.Rd b/man/getPosterior.Rd index 44e7560..36e0328 100644 --- a/man/getPosterior.Rd +++ b/man/getPosterior.Rd @@ -18,12 +18,16 @@ getPosterior( \item{data}{dataframe containing the information of dose and response. Default NULL Also a simulateData object can be provided.} -\item{mu_hat}{vector of estimated mean values} +\item{mu_hat}{vector of estimated mean values (per dose group).} -\item{se_hat}{vector of estimated standard deviations.} +\item{se_hat}{vector of estimated standard deviations (per dose group).} -\item{calc_ess}{tbd. Default NULL} +\item{calc_ess}{boolean variable, indicating whether effective sample size should be calculated. Default FALSE} +} +\value{ +posterior_list, a posterior list object is returned with information about (mixture) posterior distribution per dose group } \description{ -Either the patient level data or both the mu_hat as well as the sd_hat must to be provided. +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. +This function calculates the posterior for every dose group independently via the RBesT function postmix. } diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd index 02dd385..545bcc3 100644 --- a/man/performBayesianMCP.Rd +++ b/man/performBayesianMCP.Rd @@ -7,15 +7,16 @@ performBayesianMCP(posterior_list, contr, crit_prob_adj) } \arguments{ -\item{posterior_list}{a getPosterior object} +\item{posterior_list}{a getPosterior object with information about the (mixture) posterior distribution per dose group} \item{contr}{a getContrMat object, contrast matrix to be used for the testing step.} \item{crit_prob_adj}{a getCritProb object, specifying the critical value to be used for the testing (on the probability scale)} } \value{ -b_mcp test result +b_mcp test result, with information about p-values for the individual dose-response shapes } \description{ -performs bayesian MCP Test step. +performs bayesian MCP Test step, as described in Fleischer et al. (Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670.) +Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate is calculated that the contrast is bigger than 0 given the data observed } From cb4cf39ef75ca89ecdb587b5ba0b0c863745f760 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Fri, 15 Dec 2023 14:19:09 +0000 Subject: [PATCH 15/39] Additional documentation around getContr and other functions --- R/BMCPMod.R | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 446ecc7..023d3c6 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -103,11 +103,22 @@ assessDesign <- function ( #' @title getContr #' -#' @description This function calculates contrast vectors that are optimal for detecting certain alternatives. More information and link to publication will be added. +#' @description This function calculates contrast vectors that are optimal for detecting certain alternatives via applying the function optContr of the DoseFinding package. +#' Hereby 4 different options can be distinguished that are automatically executed based on the input that is provided +#' i) Bayesian approach: If dose_weights and a prior_list are provided an optimized contrasts for the posterior sample size is calculated. +#' In detail, in a first step the dose_weights (typically the number of patients per dose group) and the prior information is combined by calculating for +#' each dose group a posterior effective sample. Based on this posterior effective sample sizes the allocation ratio is derived, which allows for a calculation on +#' pseudo-optimal contrasts via regular MCPMod.are calculated from the +#' regular MCPMod for these specific weights +#' ii) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the +#' regular MCPMod for these specific weights +#' iii)Bayesian approach + re-estimation: If only a sd_posterior (i.e. variability of the posterior distribution) is provided, pseudo-optimal contrasts based on these posterior weights will be calculated +#' iv) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the +#' regular MCPMod for this specific vector of standard errors. For the actual evaluation this vector of standard errors is translated into a (diagonal) matrix of variances #' #' @param mods An object of class "Mods" as specified in the Dosefinding package. #' @param dose_levels vector containing the different doseage levels. -#' @param dose_weights Vector specifying weights for the different doses. Default NULL +#' @param dose_weights Vector specifying weights for the different doses. Please note that in case this information should be provided Default NULL #' @param prior_list a prior_list object. Default NULL #' @param sd_posterior a vector of positive numerics. Default NULL #' @param se_new_trial a vector of positive numerics. Default NULL @@ -223,11 +234,11 @@ getCritProb <- function ( #' @title performBayesianMCPMod #' -#' @description performs bayesian MCP Test step and modelling. +#' @description performs bayesian MCP Test step and modelling in a combined fashion. See performBayesianMCP function for MCT Test step and getModelFits for the modelling step #' -#' @param posterior_list a getPosterior object +#' @param posterior_list a getPosterior object with information about the (mixture) posterior distribution per dose group #' @param contr a getContrMat object, contrast matrix to be used for the testing step. -#' @param crit_prob_adj a getCritProb object +#' @param crit_prob_adj a getCritProb object, specifying the critical value to be used for the testing (on the probability scale). #' @param simple boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE. #' #' @return bmcpmod test result as well as modelling result. @@ -322,7 +333,8 @@ addSignificance <- function ( #' @title performBayesianMCP #' #' @description performs bayesian MCP Test step, as described in Fleischer et al. (Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670.) -#' Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate is calculated that the contrast is bigger than 0 given the data observed +#' Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate the posterior probability is calculated that the contrast is bigger than 0 (based on the posterior distribution of the dose groups). +#' In order to obtain significant test decision we consider the maximum of the posterior probabilities across the different models. This maximum is compared with a (multiplicity adjusted) critical value (on the probability scale). #' #' @param posterior_list a getPosterior object with information about the (mixture) posterior distribution per dose group #' @param contr a getContrMat object, contrast matrix to be used for the testing step. From 6d52c2a3f54883f9c134c07c4ad694d79a19c64b Mon Sep 17 00:00:00 2001 From: wojcieko Date: Fri, 15 Dec 2023 23:37:49 +0100 Subject: [PATCH 16/39] - fixed bug in getCritProb - fixed a bug in getPosteriorI regarding ESS attribute. fyi: ifelse seems to have a bug so that it only assignes only the first value of a vector - rounded ESS in assessDesign --- R/BMCPMod.R | 15 +++++++-------- R/posterior.R | 18 ++++++++++++------ 2 files changed, 19 insertions(+), 14 deletions(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 023d3c6..efeeede 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -95,7 +95,7 @@ assessDesign <- function ( attr(eval_design, "placEff") <- attr(mods, "placEff") attr(eval_design, "maxEff") <- attr(mods, "maxEff") attr(eval_design, "sampleSize") <- n_patients - attr(eval_design, "priorESS") <- getESS(prior_list) + attr(eval_design, "priorESS") <- round(getESS(prior_list), 1) return (eval_design) @@ -210,17 +210,16 @@ getCritProb <- function ( mods, dose_levels, - dose_weights =NULL, - se_new_trial = NULL, + dose_weights = NULL, + se_new_trial = NULL, alpha_crit_val = 0.025 ) { - contr <- getContr(mods = mods, - dose_levels = dose_levels , - dose_weights = dose_weights, - se_new_trial = se_new_trial, - alpha_crit_val = alpha_crit_val) + contr <- getContr(mods = mods, + dose_levels = dose_levels , + dose_weights = dose_weights, + se_new_trial = se_new_trial) crit_prob <- stats::pnorm(DoseFinding::critVal( corMat = contr$corMat, diff --git a/R/posterior.R b/R/posterior.R index c29e011..fcc67ab 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -87,12 +87,18 @@ getPosteriorI <- function( } - names(post_list) <- names(prior_list) - class(post_list) <- "postList" - attr(post_list, "ess") <- ifelse( - test = calc_ess, - yes = getESS(post_list), - no = numeric(0)) + names(post_list) <- names(prior_list) + class(post_list) <- "postList" + + if (calc_ess) { + + attr(post_list, "ess") <- getESS(post_list) + + } else { + + attr(post_list, "ess") <- numeric(0) + + } return (post_list) From fff41e1571019cccf231922e2d76b66743a88113 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Mon, 18 Dec 2023 07:42:34 +0000 Subject: [PATCH 17/39] Changed the trial number for the vignette example --- vignettes/analysis_normal.Rmd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 70a70af..df6eb97 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -41,7 +41,7 @@ histcontrol <- filter(dataset, dose == 0, primtime == 8, indication == "MAJOR DEPRESSIVE DISORDER", - protid != 6) + protid != 5) ##Create MAP Prior hist_data <- data.frame( @@ -142,7 +142,7 @@ mods <- DoseFinding::Mods( new_trial <- filter(dataset, primtime == 8, indication == "MAJOR DEPRESSIVE DISORDER", - protid == 6) + protid == 5) n_patients <- c(150, 142, 147, 149) From f852df6d77455f7118b1d0e7dd978287a8aed64d Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Mon, 18 Dec 2023 13:35:18 +0000 Subject: [PATCH 18/39] Vignette updated. --- R/modelling.R | 4 +- vignettes/analysis_normal.Rmd | 135 ++++++++++++++++++++++------------ 2 files changed, 90 insertions(+), 49 deletions(-) diff --git a/R/modelling.R b/R/modelling.R index ed85a16..66f9302 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -1,6 +1,6 @@ #' @title getModelFits #' -#' @description Fits dose-response curves for the specified dose-repsonse models, based on the posterior distributions. +#' @description Fits dose-response curves for the specified dose-response models, based on the posterior distributions. #' For the simplified fit, multivariate normal distributions will be approximated and reduced by one-dimensional normal distributions. #' For the default case, the Nelder-Mead algorithm is used. Will be further updated and links to publication as well as references will be added. #' @@ -17,6 +17,7 @@ getModelFits <- function ( models, dose_levels, posterior, + #avg_fit=FALSE, if possible we should add the average fit directly here simple = FALSE ) { @@ -32,6 +33,7 @@ getModelFits <- function ( attr(model_fits, "posterior") <- posterior class(model_fits) <- "modelFits" + return (model_fits) } diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index df6eb97..c597fa2 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -28,12 +28,12 @@ In this vignette we will show the use of the Bayesian MCPMod package for continu Hereby the focus is on the utilization of an informative prior and the BayesianMCPMod evaluation of a single trial. We will use data that is included in the clinDR package. More specifically trial results for BRINTELLIX will be used to illustrate the specification of an informative prior and the usage of such a prior for the bayesian evaluation of a (hypothetical) new trial. -More information around BRINTELLIX will be added... +BRINTELLIX is a medication used to treat major depressive disorder. Various clinical trials with different dose groups (including a control group) were conducted. # Calculation of a MAP prior -In a first step a meta analytic prior will be calculated using historical data from 4 trials (with main endpoint CfB in MADRS score after 8 weeks). -Please note that only information from the control group will be integrated (leading to an informative multi-component prior for the control group), while for the active groups a non-informative prior will be specified. +In a first step a meta analytic prior will be calculated using historical data from 4 trials (with main endpoint Change from baseline in MADRS score after 8 weeks). +Please note that only information from the control group will be integrated (leading to an informative mixture prior for the control group), while for the active groups a non-informative prior will be specified. ```{r Historical Data for Control Arm} data("metaData") dataset <- filter(as.data.frame(metaData), bname == "BRINTELLIX") @@ -113,10 +113,13 @@ dose_levels <- c(0, 2.5, 5, 10) prior_list <- getPriorList( hist_data = hist_data, - dose_levels = dose_levels) + dose_levels = dose_levels, + robust_weight = 0.3) + +getESS(prior_list) ``` # Specifications new trial -We will use the trial with ct.gov number NCT00635219 as exemplary new trial. +We will use the trial with ct.gov number NCT00735709 as exemplary new trial. To be able to apply the bayesian MCPMod approach, candidate models need to be specified. Since there are only 3 active dose levels we will limit the set of candidate models to a linear, an exponential, and an emax model. Note that the linear candidate model does not require a guesstimate. @@ -144,32 +147,21 @@ new_trial <- filter(dataset, indication == "MAJOR DEPRESSIVE DISORDER", protid == 5) -n_patients <- c(150, 142, 147, 149) +n_patients <- c(128, 124, 129, 122) ``` # Combination of prior information and trial results -As outlined in citePaper, in a first step the posterior is calculated combining the prior information with the estimated results of the new trial. +As outlined in Fleischer et al. [Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670.], in a first step the posterior is calculated combining the prior information with the estimated results of the new trial. Via the summary function it is possible to print out the summary information of the posterior distributions. -The simulation takes the variability of the historical data into account. ```{r Trial results} -#Simulation of new trial -##Note: This step should not be required, as we provide summary measures directly from the new trial -data_emax <- simulateData( - n_patients = n_patients, - dose_levels = dose_levels, - sd = with(hist_data, sum(sd * n) / sum(n)), - mods = mods, - true_model = "emax") - -posterior <- getPosterior(prior = prior_list, - data = data_emax) posterior <- getPosterior(prior = prior_list, mu_hat = new_trial$rslt, - se_hat = new_trial$se) + se_hat = new_trial$se, + calc_ess = TRUE) summary(posterior) @@ -177,33 +169,54 @@ summary(posterior) ``` # Execution of Bayesian MCPMod Test step -For the execution of the testing step of Bayesian MCPMod a critical value (on the probability scale) will be determined for a given alpha level. In addition a contrast matrix is created. Please note that there are different possibilities how to generate contrasts. -This information is then used as input for the Bayesian MCP testing function. +For the execution of the testing step of Bayesian MCPMod a critical value (on the probability scale) will be determined for a given alpha level. This critical value is calculated using (the re-estimated) contrasts for the frequentist MCPMod to ensure that when using non-informative priors the actual error level for falsely declaring a significant trial in the Bayesian MCPMod is controlled (by the specified alpha level). + +A pseudo-optimal contrast matrix is generated based on the variability of the posterior distribution (see Fleischer et al. 2022 for more details). Please note that there are different possibilities how to establish the contrasts. + + -```{r Execution of Bayesian MCPMod Test step} +```{r Preparation of input for Bayesian MCPMod Test step} crit_pval <- getCritProb( - mods = mods, - dose_levels = dose_levels, - dose_weights = n_patients, - alpha_crit_val = 0.1) + mods = mods, + dose_levels = dose_levels, + se_new_trial = new_trial$se, + alpha_crit_val = 0.05) + +contr_mat<- getContr( + mods = mods, + dose_levels = dose_levels, + sd_posterior = summary(posterior)[,2]) + +#Alternative ways to come up with a contrast would be to use +#i) the frequentist contrast +#contr_mat_prior <- getContr( +# mods = mods, +# dose_levels = dose_levels, +# dose_weights = n_patients, +# prior_list = prior_list) +# ii) re-estimated frequentist contrasts +# contr_mat_prior <- getContr( +# mods = mods, +# dose_levels = dose_levels, +# se_new_trial = new_trial$se) +# iii) Bayesian approach using number of patients for new trial and prior distribution +#contr_mat_prior <- getContr( +# mods = mods, +# dose_levels = dose_levels, +# dose_weights = n_patients, +# prior_list = prior_list) +#This would be the most reasonable output, but since it is not exported it is currently not usable. +#BMCP_result<-BayesMCPi(posterior_i = posterior, + # contr_mat = contr_mat_prior, + # crit_prob = crit_pval) -contr_mat_prior <- getContr( - mods = mods, - dose_levels = dose_levels, - dose_weights = n_patients, - prior_list = prior_list) #This would be the most reasonable output, but since it is not exported it is currently not usable. #BMCP_result<-BayesMCPi(posterior_i = posterior, # contr_mat = contr_mat_prior, # crit_prob = crit_pval) -BMCP_result <- performBayesianMCP( - posterior_list = posterior, - contr = contr_mat_prior, - crit_prob_adj = crit_pval) -BMCP_result #BMCP_result2 <- performBayesianMCPMod( # posteriors_list = posterior_emax, @@ -212,13 +225,25 @@ BMCP_result #BMCP_result2 ``` -The testing step is significant indicating a non-flat dose-response shape. In detail the p-values for the emax model is the most significant one. +The bayesian MCP testing step is then executed based on the posterior information, the provided contrasts and the (multiplicity adjusted) critical value. +```{r Executionr of Bayesian MCPMod Test step} +BMCP_result <- performBayesianMCP( + posterior_list = posterior, + contr = contr_mat, + crit_prob_adj = crit_pval) + +BMCP_result +``` + +The testing step is significant indicating a non-flat dose-response shape. In detail, all 3 models are significant and the p-value for the emax model is the most significant one. # Model fitting and visualization of results -In the model fitting step the posterior distribution is used as basis. Both simplified and full fitting are performed. -Please note that all models are fitted to the data even though they were not significant in the testing step. -For the plotting bootstrap based credible intervals (for 50% and 95%) are shown as well. +In the model fitting step the posterior distribution is used as basis. Both simplified and full fitting are performed. +For the simplified fit the multivariate normal distribution of the control group is approximated and reduced by a one-dimensional normal distribution.The actual fit (on this approximated posterior distribution) is then performed using generalized least squares criterion. +For the full fitting the (multi-dimensional) non-linear optimization problem is adressed via the Nelder Mead algorithm. + +The output of the fit includes information about the predicted effects for the included dose levels the generalized AIC (and the corresponding weights). For the considered case the simplified and the full fit are very similar. ```{r Model fitting} #Model fit post_observed <- posterior @@ -237,21 +262,35 @@ fit <- getModelFits( dose_levels = dose_levels, posterior = post_observed, simple = FALSE) +``` + +Via the predict function one can also receive estimates for dose levels that were not included in the trial. + +It is possible to plot the fitted dose response models and an AIC based average model (black lines). To assess the uncertainty one can in addition visualize credible intervals (yellow shaded areas, the default is set to 50% and 95%). These credible intervals are calculated as follows. +Samples from the posterior distribution are drawn and for every sample the simplified fitting step and a prediction is performed. These fits are then used to identify and visualize the specified quantiles. +The bootstrap based quantiles can also directly be calculated and displayed via the gotbootstrapQuantiles function. + + +```{r Visualization and interpretation of the results} +predict(fit, doses = c(0, 2.5,4, 5,7, 10)) + plot(fit, cr_bands = TRUE) plot(fit_simple, cr_bands = TRUE) -predict(fit, doses = c(4, 7)) -getBootstrapQuantiles(fit, quantiles = c(0.025, 0.975), doses = c(4, 7)) +getBootstrapQuantiles(fit, quantiles = c(0.025,0.5, 0.975), doses = c(0, 2.5,4, 5,7, 10)) ``` + +Technical note: The median quantile of the bootstrap based procedure is not necessary similar to the main model fit, as they are derived via different procedures (main fit, i.e. black lines in the plot, show the best fit of a certain model based on minimizing the residuals for the posterior distribution while bootstrap based 50% quantile indicates the median fit of the random sampling and fitting procedure). + # Additional notes -It is also possible to perform the testing and modelling step in combined fashion via the performBayesianMCPMod function +It is also possible to perform the testing and modelling step in combined fashion via the performBayesianMCPMod function. ```{r} -performBayesianMCPMod( - posterior_list = posterior, - contr = contr_mat_prior, - crit_prob_adj = crit_pval) +#performBayesianMCPMod( + # posterior_list = posterior, + # contr = contr_mat, + # crit_prob_adj = crit_pval) ``` From 1700b38c72879ca310e8044f2617a7e34f32c5e4 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Mon, 18 Dec 2023 15:01:04 +0000 Subject: [PATCH 19/39] Updates around documentation of functions and vignette. --- R/BMCPMod.R | 23 +++++++++++++++-------- R/bootstrapping.R | 11 +++++++++-- R/modelling.R | 12 +++++++++--- R/s3methods.R | 8 ++++---- man/getContr.Rd | 23 +++++++++++++++++------ man/getCritProb.Rd | 12 ++++++++++-- man/getModelFits.Rd | 13 ++++++++++--- man/performBayesianMCP.Rd | 5 +++-- man/performBayesianMCPMod.Rd | 6 +++--- vignettes/analysis_normal.Rmd | 14 +++++++++----- 10 files changed, 89 insertions(+), 38 deletions(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index efeeede..ff7e72f 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -117,11 +117,11 @@ assessDesign <- function ( #' regular MCPMod for this specific vector of standard errors. For the actual evaluation this vector of standard errors is translated into a (diagonal) matrix of variances #' #' @param mods An object of class "Mods" as specified in the Dosefinding package. -#' @param dose_levels vector containing the different doseage levels. -#' @param dose_weights Vector specifying weights for the different doses. Please note that in case this information should be provided Default NULL -#' @param prior_list a prior_list object. Default NULL -#' @param sd_posterior a vector of positive numerics. Default NULL -#' @param se_new_trial a vector of positive numerics. Default NULL +#' @param dose_levels vector containing the different dosage levels. +#' @param dose_weights Vector specifying weights for the different doses. Please note that in case this information is provided together with a prior (i.e. Option i) is planned) these two inputs should be provided on the same scale (e.g. patient numbers). Default NULL +#' @param prior_list a prior_list object, only required as input for Option i). Default NULL +#' @param sd_posterior a vector of positive values with information about the variability of the posterior distribution, only required for Option iii). Default NULL +#' @param se_new_trial a vector of positive values with information about the observed variability, only required for Option iv). Default NULL #' #' @return contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. #' @@ -196,11 +196,18 @@ getContr <- function ( #' @title getCritProb #' -#' @description This function calculates multiplicity adjusted +#' @description This function calculates multiplicity adjusted critical values. The critical values are calculated in such a way that +#' when using non-informative priors the actual error level for falsely declaring a significant trial in the Bayesian MCPMod is controlled (by the specified alpha level). +#' Hereby optimal contrasts of the frequentist MCPMod are applied and two options can be distinguished +#' i) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the +#' regular MCPMod for these specific weights and the corresponding critical value for this set of contrasts is calculated via the critVal function of the DoseFinding package. +#' ii) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the +#' regular MCPMod for this specific vector of standard errors. Here as well the critical value for this set of contrasts is calculated via the critVal function of the DoseFinding package. #' #' @param mods An object of class "Mods" as specified in the Dosefinding package. #' @param dose_levels vector containing the different dosage levels. -#' @param dose_weights Vector specifying weights for the different doses +#' @param dose_weights Vector specifying weights for the different doses, only required for Option i). Default NULL +#' @param se_new_trial a vector of positive values, only required for Option ii). Default NULL #' @param alpha_crit_val significance level. Default set to 0.025. #' #' @return crit_pval multiplicity adjusted critical value on the probability scale. @@ -339,7 +346,7 @@ addSignificance <- function ( #' @param contr a getContrMat object, contrast matrix to be used for the testing step. #' @param crit_prob_adj a getCritProb object, specifying the critical value to be used for the testing (on the probability scale) #' -#' @return b_mcp test result, with information about p-values for the individual dose-response shapes +#' @return b_mcp test result, with information about p-values for the individual dose-response shapes and overall significance #' #' @export performBayesianMCP <- function( diff --git a/R/bootstrapping.R b/R/bootstrapping.R index 0957703..2829637 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -1,7 +1,14 @@ #' @title getBootstrapQuantiles #' -#' @param model_fits tbd -#' @param quantiles tbd +#' @description A function to Calculate credible intervals to assess the uncertainty for the model fit. one can in addition visualize credible intervals (yellow shaded areas, the default is set to 50% and 95%). These credible intervals are calculated as follows. +#' Samples from the posterior distribution are drawn and for every sample the simplified fitting step and a prediction is performed. These fits are then used to identify and visualize the specified quantiles. +#' The bootstrap based quantiles can also directly be calculated and displayed via the gotbootstrapQuantiles function.dose-response curves for the specified dose-response models, based on the posterior distributions. +#' For the simplified fit, multivariate normal distributions will be approximated and reduced by one-dimensional normal distributions. +#' For the default case, the Nelder-Mead algorithm is used. +#' +#' +#' @param model_fits an object of class modelFits, i.e. information about fitted models & corresponding model coefficients as well as the posterior distribution that was the basis for the model fitting +#' @param quantiles a vector of quantiles that should be evaluated #' @param n_samples tbd #' @param doses tbd #' @param avg_fit tbd diff --git a/R/modelling.R b/R/modelling.R index 66f9302..343206d 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -2,14 +2,20 @@ #' #' @description Fits dose-response curves for the specified dose-response models, based on the posterior distributions. #' For the simplified fit, multivariate normal distributions will be approximated and reduced by one-dimensional normal distributions. -#' For the default case, the Nelder-Mead algorithm is used. Will be further updated and links to publication as well as references will be added. -#' +#' For the default case, the Nelder-Mead algorithm is used. +#' In detail, for both approaches the mean vector and the covariance $\Sigma_{l}^{-1}$ of the mixture posterior distributions and the corresponding posterior weights $\omega_{l}^{*}$ for $l \in {1,...,L}$ is used as basis +#' For the full fit a GLS estimator is used to minimize the following expression for the respective dose-response models $m$ +#' $$ +#' \hat{\theta}_{m}=argmin_{\theta_{m}} \sum_{i=l}^{L^{*}} \omega_{l}^{*}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))'\Sigma_{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m})) +#' $$ +#' Herefore the function nloptr of the nloptr package is utilized. +#' In the simplified case $L=1$ as the dimension of the posterior is reduced to 1 first. #' @param models list of model names for which a fit will be performed. #' @param dose_levels a vector containing the different dosage levels. #' @param posterior a getPosterior object, containing the (multivariate) posterior distribution per dosage level. #' @param simple boolean variable, defining whether simplified fit will be applied. Default FALSE. #' -#' @return model_fits returns a list, containing information about the fitted model coefficients, the prediction per dose group as well as maximum effect and generalized AIC per model. +#' @return model_fits returns a list, containing information about the fitted model coefficients, the prediction per dose group as well as maximum effect and generalized AIC (and corresponding weight) per model. #' #' @export getModelFits <- function ( diff --git a/R/s3methods.R b/R/s3methods.R index f4cc5f0..1f704d2 100644 --- a/R/s3methods.R +++ b/R/s3methods.R @@ -36,12 +36,12 @@ print.BayesianMCPMod <- function ( cat("Model Significance Frequencies\n") print(model_success, ...) - if (!is.na(attr(x$BayesianMCP, "ess_avg"))) { +# if (!is.na(attr(x$BayesianMCP, "ess_avg"))) {#Note SB: I have taken this out, as vignettes didn't work - cat("Average Posterior ESS\n") - print(attr(x$BayesianMCP, "ess_avg"), ...) + # cat("Average Posterior ESS\n") + # print(attr(x$BayesianMCP, "ess_avg"), ...) - } + # } } diff --git a/man/getContr.Rd b/man/getContr.Rd index 809e932..647a75d 100644 --- a/man/getContr.Rd +++ b/man/getContr.Rd @@ -16,19 +16,30 @@ getContr( \arguments{ \item{mods}{An object of class "Mods" as specified in the Dosefinding package.} -\item{dose_levels}{vector containing the different doseage levels.} +\item{dose_levels}{vector containing the different dosage levels.} -\item{dose_weights}{Vector specifying weights for the different doses. Default NULL} +\item{dose_weights}{Vector specifying weights for the different doses. Please note that in case this information is provided together with a prior (i.e. Option i) is planned) these two inputs should be provided on the same scale (e.g. patient numbers). Default NULL} -\item{prior_list}{a prior_list object. Default NULL} +\item{prior_list}{a prior_list object, only required as input for Option i). Default NULL} -\item{sd_posterior}{a vector of positive numerics. Default NULL} +\item{sd_posterior}{a vector of positive values with information about the variability of the posterior distribution, only required for Option iii). Default NULL} -\item{se_new_trial}{a vector of positive numerics. Default NULL} +\item{se_new_trial}{a vector of positive values with information about the observed variability, only required for Option iv). Default NULL} } \value{ contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. } \description{ -This function calculates contrast vectors that are optimal for detecting certain alternatives. More information and link to publication will be added. +This function calculates contrast vectors that are optimal for detecting certain alternatives via applying the function optContr of the DoseFinding package. +Hereby 4 different options can be distinguished that are automatically executed based on the input that is provided +i) Bayesian approach: If dose_weights and a prior_list are provided an optimized contrasts for the posterior sample size is calculated. +In detail, in a first step the dose_weights (typically the number of patients per dose group) and the prior information is combined by calculating for +each dose group a posterior effective sample. Based on this posterior effective sample sizes the allocation ratio is derived, which allows for a calculation on +pseudo-optimal contrasts via regular MCPMod.are calculated from the +regular MCPMod for these specific weights +ii) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the +regular MCPMod for these specific weights +iii)Bayesian approach + re-estimation: If only a sd_posterior (i.e. variability of the posterior distribution) is provided, pseudo-optimal contrasts based on these posterior weights will be calculated +iv) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the +regular MCPMod for this specific vector of standard errors. For the actual evaluation this vector of standard errors is translated into a (diagonal) matrix of variances } diff --git a/man/getCritProb.Rd b/man/getCritProb.Rd index fa648bd..e5f3335 100644 --- a/man/getCritProb.Rd +++ b/man/getCritProb.Rd @@ -17,7 +17,9 @@ getCritProb( \item{dose_levels}{vector containing the different dosage levels.} -\item{dose_weights}{Vector specifying weights for the different doses} +\item{dose_weights}{Vector specifying weights for the different doses, only required for Option i). Default NULL} + +\item{se_new_trial}{a vector of positive values, only required for Option ii). Default NULL} \item{alpha_crit_val}{significance level. Default set to 0.025.} } @@ -25,5 +27,11 @@ getCritProb( crit_pval multiplicity adjusted critical value on the probability scale. } \description{ -This function calculates multiplicity adjusted +This function calculates multiplicity adjusted critical values. The critical values are calculated in such a way that +when using non-informative priors the actual error level for falsely declaring a significant trial in the Bayesian MCPMod is controlled (by the specified alpha level). +Hereby optimal contrasts of the frequentist MCPMod are applied and two options can be distinguished +i) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the +regular MCPMod for these specific weights and the corresponding critical value for this set of contrasts is calculated via the critVal function of the DoseFinding package. +ii) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the +regular MCPMod for this specific vector of standard errors. Here as well the critical value for this set of contrasts is calculated via the critVal function of the DoseFinding package. } diff --git a/man/getModelFits.Rd b/man/getModelFits.Rd index 4a87f02..7292b80 100644 --- a/man/getModelFits.Rd +++ b/man/getModelFits.Rd @@ -16,10 +16,17 @@ getModelFits(models, dose_levels, posterior, simple = FALSE) \item{simple}{boolean variable, defining whether simplified fit will be applied. Default FALSE.} } \value{ -model_fits returns a list, containing information about the fitted model coefficients, the prediction per dose group as well as maximum effect and generalized AIC per model. +model_fits returns a list, containing information about the fitted model coefficients, the prediction per dose group as well as maximum effect and generalized AIC (and corresponding weight) per model. } \description{ -Fits dose-response curves for the specified dose-repsonse models, based on the posterior distributions. +Fits dose-response curves for the specified dose-response models, based on the posterior distributions. For the simplified fit, multivariate normal distributions will be approximated and reduced by one-dimensional normal distributions. -For the default case, the Nelder-Mead algorithm is used. Will be further updated and links to publication as well as references will be added. +For the default case, the Nelder-Mead algorithm is used. +In detail, for both approaches the mean vector and the covariance $\Sigma_{l}^{-1}$ of the mixture posterior distributions and the corresponding posterior weights $\omega_{l}^{\emph{}$ for $l \in {1,...,L}$ is used as basis +For the full fit a GLS estimator is used to minimize the following expression for the respective dose-response models $m$ +$$ +\hat{\theta}\emph{{m}=argmin}{\theta_{m}} \sum_{i=l}^{L^{}}} \omega_{l}^{*}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}\emph{{m}))'\Sigma}{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m})) +$$ +Herefore the function nloptr of the nloptr package is utilized. +In the simplified case $L=1$ as the dimension of the posterior is reduced to 1 first. } diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd index 545bcc3..702a9c5 100644 --- a/man/performBayesianMCP.Rd +++ b/man/performBayesianMCP.Rd @@ -14,9 +14,10 @@ performBayesianMCP(posterior_list, contr, crit_prob_adj) \item{crit_prob_adj}{a getCritProb object, specifying the critical value to be used for the testing (on the probability scale)} } \value{ -b_mcp test result, with information about p-values for the individual dose-response shapes +b_mcp test result, with information about p-values for the individual dose-response shapes and overall significance } \description{ performs bayesian MCP Test step, as described in Fleischer et al. (Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670.) -Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate is calculated that the contrast is bigger than 0 given the data observed +Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate the posterior probability is calculated that the contrast is bigger than 0 (based on the posterior distribution of the dose groups). +In order to obtain significant test decision we consider the maximum of the posterior probabilities across the different models. This maximum is compared with a (multiplicity adjusted) critical value (on the probability scale). } diff --git a/man/performBayesianMCPMod.Rd b/man/performBayesianMCPMod.Rd index 6359109..58d7249 100644 --- a/man/performBayesianMCPMod.Rd +++ b/man/performBayesianMCPMod.Rd @@ -7,11 +7,11 @@ performBayesianMCPMod(posterior_list, contr, crit_prob_adj, simple = FALSE) } \arguments{ -\item{posterior_list}{a getPosterior object} +\item{posterior_list}{a getPosterior object with information about the (mixture) posterior distribution per dose group} \item{contr}{a getContrMat object, contrast matrix to be used for the testing step.} -\item{crit_prob_adj}{a getCritProb object} +\item{crit_prob_adj}{a getCritProb object, specifying the critical value to be used for the testing (on the probability scale).} \item{simple}{boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.} } @@ -19,5 +19,5 @@ performBayesianMCPMod(posterior_list, contr, crit_prob_adj, simple = FALSE) bmcpmod test result as well as modelling result. } \description{ -performs bayesian MCP Test step and modelling. +performs bayesian MCP Test step and modelling in a combined fashion. See performBayesianMCP function for MCT Test step and getModelFits for the modelling step } diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index c597fa2..d6863c5 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -285,12 +285,16 @@ getBootstrapQuantiles(fit, quantiles = c(0.025,0.5, 0.975), doses = c(0, 2.5,4, Technical note: The median quantile of the bootstrap based procedure is not necessary similar to the main model fit, as they are derived via different procedures (main fit, i.e. black lines in the plot, show the best fit of a certain model based on minimizing the residuals for the posterior distribution while bootstrap based 50% quantile indicates the median fit of the random sampling and fitting procedure). # Additional notes -It is also possible to perform the testing and modelling step in combined fashion via the performBayesianMCPMod function. +It is also possible to perform the testing and modelling step in a combined fashion via the performBayesianMCPMod function. + ```{r} -#performBayesianMCPMod( - # posterior_list = posterior, - # contr = contr_mat, - # crit_prob_adj = crit_pval) +overall_result<-performBayesianMCPMod( + posterior_list = posterior, + contr = contr_mat, + crit_prob_adj = crit_pval, + simple = FALSE) + +overall_result$Mod ``` From 44ccf9e128d7558a9c06564e8459bc37880592d6 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Mon, 18 Dec 2023 15:26:36 +0000 Subject: [PATCH 20/39] Additional updates (e.g. to fix bug around print of BayesianMCPMod) --- R/bootstrapping.R | 19 ++++++++++--------- R/plot.R | 2 +- R/s3methods.R | 8 ++++---- 3 files changed, 15 insertions(+), 14 deletions(-) diff --git a/R/bootstrapping.R b/R/bootstrapping.R index 2829637..bc2d7e6 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -1,19 +1,20 @@ #' @title getBootstrapQuantiles #' -#' @description A function to Calculate credible intervals to assess the uncertainty for the model fit. one can in addition visualize credible intervals (yellow shaded areas, the default is set to 50% and 95%). These credible intervals are calculated as follows. -#' Samples from the posterior distribution are drawn and for every sample the simplified fitting step and a prediction is performed. These fits are then used to identify and visualize the specified quantiles. -#' The bootstrap based quantiles can also directly be calculated and displayed via the gotbootstrapQuantiles function.dose-response curves for the specified dose-response models, based on the posterior distributions. -#' For the simplified fit, multivariate normal distributions will be approximated and reduced by one-dimensional normal distributions. -#' For the default case, the Nelder-Mead algorithm is used. +#' @description A function for the calculation of credible intervals to assess the uncertainty for the model fit. +#' Hereby the credible intervals are calculated as follows. +#' Samples from the posterior distribution are drawn (via the RBesT function rmix) and for every sample the simplified fitting step (see getModelFits function) and a prediction is performed. +#' These fits are then used to identify the specified quantiles. +#' This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in O'Quigley, Iasonos, J., and B. Bornkamp. 2017. Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials. Boca Raton: CRC Press. +#' Instead of drawing n bootstrap samples from the sampling distribution of the trial dose-response estimates, here the samples are directly taken from the posterior distribution. #' #' #' @param model_fits an object of class modelFits, i.e. information about fitted models & corresponding model coefficients as well as the posterior distribution that was the basis for the model fitting #' @param quantiles a vector of quantiles that should be evaluated -#' @param n_samples tbd -#' @param doses tbd -#' @param avg_fit tbd +#' @param n_samples number of samples that should be drawn as basis for the +#' @param doses a vector of doses for which a prediction should be performed +#' @param avg_fit boolean variable, defining whether an average fit (based on generalized AIC weights) should be performed in addition to the individual models. Default TRUE. #' -#' @return tbd +#' @return A data frame with entries doses, models, and quantiles #' @export getBootstrapQuantiles <- function ( diff --git a/R/plot.R b/R/plot.R index 26e627a..93a60fc 100644 --- a/R/plot.R +++ b/R/plot.R @@ -1,7 +1,7 @@ #' @title plot.modelFits #' #' @description plot function based on the ggplot2 package. Providing visualizations for each model and a average Fit. -#' More details to be added, as well as references. +#' Black line shows the fitted dose response models and an AIC based average model (black lines). To assess the uncertainty one can in addition visualize credible intervals (yellow shaded areas, the default is set to 50% and 95%). These credible intervals are calculated as follows. #' @param x an object of type getModelFits #' @param gAIC logical value indicating whether gAIC values are shown in the plot. Default TRUE #' @param avg_fit logical value indicating whether average fit is presented in the plot. Default TRUE diff --git a/R/s3methods.R b/R/s3methods.R index 1f704d2..5c94755 100644 --- a/R/s3methods.R +++ b/R/s3methods.R @@ -36,12 +36,12 @@ print.BayesianMCPMod <- function ( cat("Model Significance Frequencies\n") print(model_success, ...) -# if (!is.na(attr(x$BayesianMCP, "ess_avg"))) {#Note SB: I have taken this out, as vignettes didn't work + if (any(!is.na(attr(x$BayesianMCP, "ess_avg")))) { - # cat("Average Posterior ESS\n") - # print(attr(x$BayesianMCP, "ess_avg"), ...) + cat("Average Posterior ESS\n") + print(attr(x$BayesianMCP, "ess_avg"), ...) - # } + } } From c16c9917deb268fb7ecbc159649d2a25ca115e8f Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Mon, 18 Dec 2023 15:47:20 +0000 Subject: [PATCH 21/39] Additional updates around the function documentations. --- R/bootstrapping.R | 1 - R/plot.R | 4 +++- man/getBootstrapQuantiles.Rd | 19 ++++++++++++------- man/plot.modelFits.Rd | 4 +++- vignettes/analysis_normal.Rmd | 4 ++-- 5 files changed, 20 insertions(+), 12 deletions(-) diff --git a/R/bootstrapping.R b/R/bootstrapping.R index bc2d7e6..f9d8ad0 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -7,7 +7,6 @@ #' This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in O'Quigley, Iasonos, J., and B. Bornkamp. 2017. Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials. Boca Raton: CRC Press. #' Instead of drawing n bootstrap samples from the sampling distribution of the trial dose-response estimates, here the samples are directly taken from the posterior distribution. #' -#' #' @param model_fits an object of class modelFits, i.e. information about fitted models & corresponding model coefficients as well as the posterior distribution that was the basis for the model fitting #' @param quantiles a vector of quantiles that should be evaluated #' @param n_samples number of samples that should be drawn as basis for the diff --git a/R/plot.R b/R/plot.R index 93a60fc..d2cd89f 100644 --- a/R/plot.R +++ b/R/plot.R @@ -1,7 +1,9 @@ #' @title plot.modelFits #' #' @description plot function based on the ggplot2 package. Providing visualizations for each model and a average Fit. -#' Black line shows the fitted dose response models and an AIC based average model (black lines). To assess the uncertainty one can in addition visualize credible intervals (yellow shaded areas, the default is set to 50% and 95%). These credible intervals are calculated as follows. +#' Black lines show the fitted dose response models and an AIC based average model. Dots indicate the posterior median and vertical lines show corresponding credible intervals (i.e. the variability of the posterior distribution of the respective dose group). +#' To assess the uncertainty of the model fit one can in addition visualize credible bands (yellow/orange shaded areas). The calculation of these bands is performed via the getBootstrapQuantiles function. +#' The default setting is that these credible bands are not calculated. #' @param x an object of type getModelFits #' @param gAIC logical value indicating whether gAIC values are shown in the plot. Default TRUE #' @param avg_fit logical value indicating whether average fit is presented in the plot. Default TRUE diff --git a/man/getBootstrapQuantiles.Rd b/man/getBootstrapQuantiles.Rd index d7f8ecd..7744bc5 100644 --- a/man/getBootstrapQuantiles.Rd +++ b/man/getBootstrapQuantiles.Rd @@ -13,19 +13,24 @@ getBootstrapQuantiles( ) } \arguments{ -\item{model_fits}{tbd} +\item{model_fits}{an object of class modelFits, i.e. information about fitted models & corresponding model coefficients as well as the posterior distribution that was the basis for the model fitting} -\item{quantiles}{tbd} +\item{quantiles}{a vector of quantiles that should be evaluated} -\item{n_samples}{tbd} +\item{n_samples}{number of samples that should be drawn as basis for the} -\item{doses}{tbd} +\item{doses}{a vector of doses for which a prediction should be performed} -\item{avg_fit}{tbd} +\item{avg_fit}{boolean variable, defining whether an average fit (based on generalized AIC weights) should be performed in addition to the individual models. Default TRUE.} } \value{ -tbd +A data frame with entries doses, models, and quantiles } \description{ -getBootstrapQuantiles +A function for the calculation of credible intervals to assess the uncertainty for the model fit. +Hereby the credible intervals are calculated as follows. +Samples from the posterior distribution are drawn (via the RBesT function rmix) and for every sample the simplified fitting step (see getModelFits function) and a prediction is performed. +These fits are then used to identify the specified quantiles. +This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in O'Quigley, Iasonos, J., and B. Bornkamp. 2017. Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials. Boca Raton: CRC Press. +Instead of drawing n bootstrap samples from the sampling distribution of the trial dose-response estimates, here the samples are directly taken from the posterior distribution. } diff --git a/man/plot.modelFits.Rd b/man/plot.modelFits.Rd index 0cb1c26..be11a8a 100644 --- a/man/plot.modelFits.Rd +++ b/man/plot.modelFits.Rd @@ -43,5 +43,7 @@ plts returns a ggplot2 object } \description{ plot function based on the ggplot2 package. Providing visualizations for each model and a average Fit. -More details to be added, as well as references. +Black lines show the fitted dose response models and an AIC based average model. Dots indicate the posterior median and vertical lines show corresponding credible intervals (i.e. the variability of the posterior distribution of the respective dose group). +To assess the uncertainty of the model fit one can in addition visualize credible bands (yellow/orange shaded areas). The calculation of these bands is performed via the getBootstrapQuantiles function. +The default setting is that these credible bands are not calculated. } diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index d6863c5..5253bf8 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -266,9 +266,9 @@ fit <- getModelFits( Via the predict function one can also receive estimates for dose levels that were not included in the trial. -It is possible to plot the fitted dose response models and an AIC based average model (black lines). To assess the uncertainty one can in addition visualize credible intervals (yellow shaded areas, the default is set to 50% and 95%). These credible intervals are calculated as follows. +It is possible to plot the fitted dose response models and an AIC based average model (black lines). To assess the uncertainty one can in addition visualize credible bands (yellow/orange shaded areas, the default is set to 50% and 95%). These credible bands are calculated as follows. Samples from the posterior distribution are drawn and for every sample the simplified fitting step and a prediction is performed. These fits are then used to identify and visualize the specified quantiles. -The bootstrap based quantiles can also directly be calculated and displayed via the gotbootstrapQuantiles function. +The bootstrap based quantiles can also be directly calculated and displayed via the getbootstrapQuantiles function. ```{r Visualization and interpretation of the results} From 20e0baa2d84a399df33c9bf57d418d069f733508 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Mon, 18 Dec 2023 19:51:27 +0000 Subject: [PATCH 22/39] Additional examples for documentation added --- R/BMCPMod.R | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++- R/posterior.R | 13 ++++++++++ 2 files changed, 78 insertions(+), 1 deletion(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index ff7e72f..237352c 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -123,6 +123,17 @@ assessDesign <- function ( #' @param sd_posterior a vector of positive values with information about the variability of the posterior distribution, only required for Option iii). Default NULL #' @param se_new_trial a vector of positive values with information about the observed variability, only required for Option iv). Default NULL #' +#' @examples +#' # example code +#' models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +#' doses = c(0, 0.5, 2,4, 8)) +#' dose_levels=c(0, 0.5, 2, 4, 8) +#' sd_posterior = c(2.8,3,2.5,3.5,4) +#' contr_mat<- getContr( +#' mods = models, +#' dose_levels = dose_levels, +#' sd_posterior = sd_posterior) +#' #' @return contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. #' #' @export @@ -209,7 +220,17 @@ getContr <- function ( #' @param dose_weights Vector specifying weights for the different doses, only required for Option i). Default NULL #' @param se_new_trial a vector of positive values, only required for Option ii). Default NULL #' @param alpha_crit_val significance level. Default set to 0.025. -#' +#' +#' @examples +#' # example code +#' models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +#' doses = c(0, 0.5, 2,4, 8)) +#' dose_levels=c(0, 0.5, 2, 4, 8) +#' critVal<- getCritProb( +#' mods = models, +#' dose_weights =c(50,50,50,50,50), #reflecting the planned sample size +#' dose_levels = dose_levels, +#' alpha_crit_val = 0.05) #' @return crit_pval multiplicity adjusted critical value on the probability scale. #' #' @export @@ -246,6 +267,27 @@ getCritProb <- function ( #' @param contr a getContrMat object, contrast matrix to be used for the testing step. #' @param crit_prob_adj a getCritProb object, specifying the critical value to be used for the testing (on the probability scale). #' @param simple boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE. +#' @examples +#' # example code +#' models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +#' doses = c(0, 0.5, 2,4, 8)) +#' dose_levels=c(0, 0.5, 2, 4, 8) +#' sd_posterior = c(2.8,3,2.5,3.5,4) +#' contr_mat<- getContr( +#' mods = models, +#' dose_levels = dose_levels, +#' sd_posterior = sd_posterior) +#' critVal<- getCritProb( +#' mods = models, +#' dose_weights =c(50,50,50,50,50), #reflecting the planned sample size +#' dose_levels = dose_levels, +#' alpha_crit_val = 0.05) +#' posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), +#' DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), +#' 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)) +#' performBayesianMCPMod(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal,simple = FALSE) #' #' @return bmcpmod test result as well as modelling result. #' @@ -346,6 +388,28 @@ addSignificance <- function ( #' @param contr a getContrMat object, contrast matrix to be used for the testing step. #' @param crit_prob_adj a getCritProb object, specifying the critical value to be used for the testing (on the probability scale) #' +#' @examples +#' # example code +#' models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +#' doses = c(0, 0.5, 2,4, 8)) +#' dose_levels=c(0, 0.5, 2, 4, 8) +#' sd_posterior = c(2.8,3,2.5,3.5,4) +#' contr_mat<- getContr( +#' mods = models, +#' dose_levels = dose_levels, +#' sd_posterior = sd_posterior) +#' critVal<- getCritProb( +#' mods = models, +#' dose_weights =c(50,50,50,50,50), #reflecting the planned sample size +#' dose_levels = dose_levels, +#' alpha_crit_val = 0.05) +#' posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), +#' DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), +#' 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)) +#' performBayesianMCP(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal) +#' #' @return b_mcp test result, with information about p-values for the individual dose-response shapes and overall significance #' #' @export diff --git a/R/posterior.R b/R/posterior.R index fcc67ab..b79ec08 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -11,6 +11,19 @@ #' @param se_hat vector of estimated standard deviations (per dose group). #' @param calc_ess boolean variable, indicating whether effective sample size should be calculated. Default FALSE #' @return posterior_list, a posterior list object is returned with information about (mixture) posterior distribution per dose group +#' @examples +#' # example code +#' prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), +#' DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), +#' 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)) +#' mu<-c(0,1,1.5,2,2.5) +#' se<-c(5,4,6,7,8) +#' getPosterior <- function( +#' prior_list = prior_list, +#' mu_hat = mu, +#' se_hat = se) #' @export getPosterior <- function( prior_list, From 6440e34c2ddf7b91058b6b51958863f532bda746 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Mon, 18 Dec 2023 20:07:33 +0000 Subject: [PATCH 23/39] Additional documentations added (and minor update in vignette) --- R/bootstrapping.R | 13 +++++++++++++ R/modelling.R | 13 +++++++++++-- R/plot.R | 16 ++++++++++++++-- vignettes/analysis_normal.Rmd | 2 +- 4 files changed, 39 insertions(+), 5 deletions(-) diff --git a/R/bootstrapping.R b/R/bootstrapping.R index f9d8ad0..1cb306d 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -13,6 +13,19 @@ #' @param doses a vector of doses for which a prediction should be performed #' @param avg_fit boolean variable, defining whether an average fit (based on generalized AIC weights) should be performed in addition to the individual models. Default TRUE. #' +#' @examples +#' # example code +#' posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2), +#' DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2), +#' DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) , +#' DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) , +#' DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2)) +#' models=c("emax","exponential","sigEmax","linear") +#' dose_levels=c(0,1,2,4,8) +#' fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels) +#' fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE) +#' getBootstrapQuantiles(fit, quantiles = c(0.025,0.5, 0.975), doses = c(0, 1,2,3,4,6,8)) +#' getBootstrapQuantiles(fit_simple, n_samples=2000, quantiles = c(0.025,0.5, 0.975), doses = c(0, 1,2,3,4,6,8)) #' @return A data frame with entries doses, models, and quantiles #' @export getBootstrapQuantiles <- function ( diff --git a/R/modelling.R b/R/modelling.R index 343206d..bd143ce 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -14,7 +14,17 @@ #' @param dose_levels a vector containing the different dosage levels. #' @param posterior a getPosterior object, containing the (multivariate) posterior distribution per dosage level. #' @param simple boolean variable, defining whether simplified fit will be applied. Default FALSE. -#' +#' @examples +#' # example code +#' posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2), +#' DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2), +#' DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) , +#' DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) , +#' DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2)) +#' models=c("emax","exponential","sigEmax","linear") +#' dose_levels=c(0,1,2,4,8) +#' fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels) +#' fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE) #' @return model_fits returns a list, containing information about the fitted model coefficients, the prediction per dose group as well as maximum effect and generalized AIC (and corresponding weight) per model. #' #' @export @@ -23,7 +33,6 @@ getModelFits <- function ( models, dose_levels, posterior, - #avg_fit=FALSE, if possible we should add the average fit directly here simple = FALSE ) { diff --git a/R/plot.R b/R/plot.R index d2cd89f..b53821c 100644 --- a/R/plot.R +++ b/R/plot.R @@ -2,7 +2,7 @@ #' #' @description plot function based on the ggplot2 package. Providing visualizations for each model and a average Fit. #' Black lines show the fitted dose response models and an AIC based average model. Dots indicate the posterior median and vertical lines show corresponding credible intervals (i.e. the variability of the posterior distribution of the respective dose group). -#' To assess the uncertainty of the model fit one can in addition visualize credible bands (yellow/orange shaded areas). The calculation of these bands is performed via the getBootstrapQuantiles function. +#' To assess the uncertainty of the model fit one can in addition visualize credible bands (default coloring as orange shaded areas). The calculation of these bands is performed via the getBootstrapQuantiles function. #' The default setting is that these credible bands are not calculated. #' @param x an object of type getModelFits #' @param gAIC logical value indicating whether gAIC values are shown in the plot. Default TRUE @@ -14,7 +14,19 @@ #' @param n_bs_smpl number of bootstrap samples being used. Default set to 1000. #' @param acc_color color of the credible bands. Default set to "orange" #' @param ... optional parameter to be passed. -#' +#' @examples +#' # example code +#' posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2), +#' DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2), +#' DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) , +#' DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) , +#' DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2)) +#' models=c("emax","exponential","sigEmax","linear") +#' dose_levels=c(0,1,2,4,8) +#' fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels) +#' fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE) +#' plot(fit, cr_bands = TRUE) +#' plot(fit_simple, cr_bands = TRUE,alpha_CrB =c(0.05,0.1,0.5)) #' @return plts returns a ggplot2 object #' @export plot.modelFits <- function ( diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 5253bf8..0107f6c 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -266,7 +266,7 @@ fit <- getModelFits( Via the predict function one can also receive estimates for dose levels that were not included in the trial. -It is possible to plot the fitted dose response models and an AIC based average model (black lines). To assess the uncertainty one can in addition visualize credible bands (yellow/orange shaded areas, the default is set to 50% and 95%). These credible bands are calculated as follows. +It is possible to plot the fitted dose response models and an AIC based average model (black lines). To assess the uncertainty one can in addition visualize credible bands (orange shaded areas, the default is set to 50% and 95%). These credible bands are calculated as follows. Samples from the posterior distribution are drawn and for every sample the simplified fitting step and a prediction is performed. These fits are then used to identify and visualize the specified quantiles. The bootstrap based quantiles can also be directly calculated and displayed via the getbootstrapQuantiles function. From 81dab25e6ab0e7c2ff9c42231d514dce5f24786b Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Mon, 18 Dec 2023 20:51:43 +0000 Subject: [PATCH 24/39] Additional updates around the documentation. --- R/modelling.R | 19 +++++++++++++------ R/s3methods.R | 17 +++++++++++++++++ man/getBootstrapQuantiles.Rd | 14 ++++++++++++++ man/getContr.Rd | 12 ++++++++++++ man/getCritProb.Rd | 11 +++++++++++ man/getModelFits.Rd | 31 +++++++++++++++++++++++++------ man/getPosterior.Rd | 14 ++++++++++++++ man/performBayesianMCP.Rd | 23 +++++++++++++++++++++++ man/performBayesianMCPMod.Rd | 23 +++++++++++++++++++++++ man/plot.modelFits.Rd | 16 +++++++++++++++- 10 files changed, 167 insertions(+), 13 deletions(-) diff --git a/R/modelling.R b/R/modelling.R index bd143ce..013dec2 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -3,13 +3,20 @@ #' @description Fits dose-response curves for the specified dose-response models, based on the posterior distributions. #' For the simplified fit, multivariate normal distributions will be approximated and reduced by one-dimensional normal distributions. #' For the default case, the Nelder-Mead algorithm is used. -#' In detail, for both approaches the mean vector and the covariance $\Sigma_{l}^{-1}$ of the mixture posterior distributions and the corresponding posterior weights $\omega_{l}^{*}$ for $l \in {1,...,L}$ is used as basis -#' For the full fit a GLS estimator is used to minimize the following expression for the respective dose-response models $m$ -#' $$ -#' \hat{\theta}_{m}=argmin_{\theta_{m}} \sum_{i=l}^{L^{*}} \omega_{l}^{*}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))'\Sigma_{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m})) -#' $$ +#' In detail, for both approaches the mean vector \eqn{\theta^{Y}} and the covariance \eqn{\Sigma} of the (mixture) posterior distributions and the corresponding posterior weights \eqn{\tilde{\omega}_{l}} for $l \in {1,...,L}$ are used as basis +#' For the full fit a GLS estimator is used to minimize the following expression for the respective dose-response models \eqn{m} +#' \deqn{ \hat{\theta}_{m}=argmin_{\theta_{m}} \sum_{i=l}^{L} \tilde{\omega}_{l}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))'\Sigma_{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))} #' Herefore the function nloptr of the nloptr package is utilized. -#' In the simplified case $L=1$ as the dimension of the posterior is reduced to 1 first. +#' In the simplified case $L=1$, as the dimension of the posterior is reduced to 1 first. +#' The generalized AIC values are calculated via the formula +#' \deqn{gAIC_{m} = \sum_{l=1}^{L} \tilde{\omega}_{l} \sum_{i=0}^{K} \frac{1}{\Sigma_{l_{i,i}}} (\theta_{l_i}^Y - f(dose_{i},\hat{\theta}_{m}))^2 + 2p } +#' where \eqn{p} denotes the number of estimated parameters and \eqn{K} the number of active dose levels. +#' Here as well for the simplified case the formula reduces to one summand as $L=1$. +#' Corresponding gAIC based weights for model $M$ are calculated as outlined in Schorning et al. 2016: +#' \deqn{ +#' \Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_(m))} +#' } +#' where \eqn{Q} denotes the number of models included in the averaging procedure. #' @param models list of model names for which a fit will be performed. #' @param dose_levels a vector containing the different dosage levels. #' @param posterior a getPosterior object, containing the (multivariate) posterior distribution per dosage level. diff --git a/R/s3methods.R b/R/s3methods.R index 5c94755..2583ca5 100644 --- a/R/s3methods.R +++ b/R/s3methods.R @@ -77,7 +77,24 @@ print.BayesianMCP <- function ( } ## ModelFits ---------------------------------------------- +#' @description This function performs model predictions based on the provided model and dose specifications +#' +#' @param object a modelFits object containing information about the fitted model coefficients +#' @param doses a vector specifying the doses for which a prediction should be done getContrMat object, contrast matrix to be used for the testing step. +#' @examples +#' # example code +#' posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2), +#' DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2), +#' DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) , +#' DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) , +#' DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2)) +#' models=c("emax","exponential","sigEmax","linear") +#' dose_levels=c(0,1,2,4,8) +#' fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels) +#' predict(fit, doses = c(0, 1,3,4,6,8)) +#' +#' @return a list with the model predictions for the specified models and doses #' @export predict.modelFits <- function ( diff --git a/man/getBootstrapQuantiles.Rd b/man/getBootstrapQuantiles.Rd index 7744bc5..cc0d1a4 100644 --- a/man/getBootstrapQuantiles.Rd +++ b/man/getBootstrapQuantiles.Rd @@ -34,3 +34,17 @@ These fits are then used to identify the specified quantiles. This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in O'Quigley, Iasonos, J., and B. Bornkamp. 2017. Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials. Boca Raton: CRC Press. Instead of drawing n bootstrap samples from the sampling distribution of the trial dose-response estimates, here the samples are directly taken from the posterior distribution. } +\examples{ +# example code +posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2), + DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2), + DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) , + DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) , + DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2)) +models=c("emax","exponential","sigEmax","linear") +dose_levels=c(0,1,2,4,8) +fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels) +fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE) +getBootstrapQuantiles(fit, quantiles = c(0.025,0.5, 0.975), doses = c(0, 1,2,3,4,6,8)) +getBootstrapQuantiles(fit_simple, n_samples=2000, quantiles = c(0.025,0.5, 0.975), doses = c(0, 1,2,3,4,6,8)) +} diff --git a/man/getContr.Rd b/man/getContr.Rd index 647a75d..3abb8cc 100644 --- a/man/getContr.Rd +++ b/man/getContr.Rd @@ -43,3 +43,15 @@ iii)Bayesian approach + re-estimation: If only a sd_posterior (i.e. variability iv) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the regular MCPMod for this specific vector of standard errors. For the actual evaluation this vector of standard errors is translated into a (diagonal) matrix of variances } +\examples{ +# example code +models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +doses = c(0, 0.5, 2,4, 8)) +dose_levels=c(0, 0.5, 2, 4, 8) +sd_posterior = c(2.8,3,2.5,3.5,4) +contr_mat<- getContr( +mods = models, +dose_levels = dose_levels, +sd_posterior = sd_posterior) + +} diff --git a/man/getCritProb.Rd b/man/getCritProb.Rd index e5f3335..e6ff1bb 100644 --- a/man/getCritProb.Rd +++ b/man/getCritProb.Rd @@ -35,3 +35,14 @@ regular MCPMod for these specific weights and the corresponding critical value f ii) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the regular MCPMod for this specific vector of standard errors. Here as well the critical value for this set of contrasts is calculated via the critVal function of the DoseFinding package. } +\examples{ +# example code +models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +doses = c(0, 0.5, 2,4, 8)) +dose_levels=c(0, 0.5, 2, 4, 8) +critVal<- getCritProb( + mods = models, + dose_weights =c(50,50,50,50,50), #reflecting the planned sample size + dose_levels = dose_levels, + alpha_crit_val = 0.05) +} diff --git a/man/getModelFits.Rd b/man/getModelFits.Rd index 7292b80..715feed 100644 --- a/man/getModelFits.Rd +++ b/man/getModelFits.Rd @@ -22,11 +22,30 @@ model_fits returns a list, containing information about the fitted model coeffic Fits dose-response curves for the specified dose-response models, based on the posterior distributions. For the simplified fit, multivariate normal distributions will be approximated and reduced by one-dimensional normal distributions. For the default case, the Nelder-Mead algorithm is used. -In detail, for both approaches the mean vector and the covariance $\Sigma_{l}^{-1}$ of the mixture posterior distributions and the corresponding posterior weights $\omega_{l}^{\emph{}$ for $l \in {1,...,L}$ is used as basis -For the full fit a GLS estimator is used to minimize the following expression for the respective dose-response models $m$ -$$ -\hat{\theta}\emph{{m}=argmin}{\theta_{m}} \sum_{i=l}^{L^{}}} \omega_{l}^{*}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}\emph{{m}))'\Sigma}{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m})) -$$ +In detail, for both approaches the mean vector \eqn{\theta^{Y}} and the covariance \eqn{\Sigma} of the (mixture) posterior distributions and the corresponding posterior weights \eqn{\tilde{\omega}_{l}} for $l \in {1,...,L}$ are used as basis +For the full fit a GLS estimator is used to minimize the following expression for the respective dose-response models \eqn{m} +\deqn{ \hat{\theta}_{m}=argmin_{\theta_{m}} \sum_{i=l}^{L} \tilde{\omega}_{l}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))'\Sigma_{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))} Herefore the function nloptr of the nloptr package is utilized. -In the simplified case $L=1$ as the dimension of the posterior is reduced to 1 first. +In the simplified case $L=1$, as the dimension of the posterior is reduced to 1 first. +The generalized AIC values are calculated via the formula +\deqn{gAIC_{m} = \sum_{l=1}^{L} \tilde{\omega}_{l} \sum_{i=0}^{K} \frac{1}{\Sigma_{l_{i,i}}} (\theta_{l_i}^Y - f(dose_{i},\hat{\theta}_{m}))^2 + 2p } +where \eqn{p} denotes the number of estimated parameters and \eqn{K} the number of active dose levels. +Here as well for the simplified case the formula reduces to one summand as $L=1$. +Corresponding gAIC based weights for model $M$ are calculated as outlined in Schorning et al. 2016: +\deqn{ +\Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_(m))} +} +where \eqn{Q} denotes the number of models included in the averaging procedure. +} +\examples{ +# example code +posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2), + DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2), + DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) , + DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) , + DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2)) +models=c("emax","exponential","sigEmax","linear") +dose_levels=c(0,1,2,4,8) +fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels) +fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE) } diff --git a/man/getPosterior.Rd b/man/getPosterior.Rd index 36e0328..005b6d0 100644 --- a/man/getPosterior.Rd +++ b/man/getPosterior.Rd @@ -31,3 +31,17 @@ posterior_list, a posterior list object is returned with information about (mixt 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. This function calculates the posterior for every dose group independently via the RBesT function postmix. } +\examples{ +# example code +prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), + DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), + 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)) +mu<-c(0,1,1.5,2,2.5) +se<-c(5,4,6,7,8) +getPosterior <- function( + prior_list = prior_list, + mu_hat = mu, + se_hat = se) +} diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd index 702a9c5..fa89d63 100644 --- a/man/performBayesianMCP.Rd +++ b/man/performBayesianMCP.Rd @@ -21,3 +21,26 @@ performs bayesian MCP Test step, as described in Fleischer et al. (Bayesian MCPM Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate the posterior probability is calculated that the contrast is bigger than 0 (based on the posterior distribution of the dose groups). In order to obtain significant test decision we consider the maximum of the posterior probabilities across the different models. This maximum is compared with a (multiplicity adjusted) critical value (on the probability scale). } +\examples{ +# example code +models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +doses = c(0, 0.5, 2,4, 8)) +dose_levels=c(0, 0.5, 2, 4, 8) +sd_posterior = c(2.8,3,2.5,3.5,4) +contr_mat<- getContr( +mods = models, +dose_levels = dose_levels, +sd_posterior = sd_posterior) +critVal<- getCritProb( + mods = models, + dose_weights =c(50,50,50,50,50), #reflecting the planned sample size + dose_levels = dose_levels, + alpha_crit_val = 0.05) +posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), + DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), + 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)) +performBayesianMCP(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal) + +} diff --git a/man/performBayesianMCPMod.Rd b/man/performBayesianMCPMod.Rd index 58d7249..e43ec82 100644 --- a/man/performBayesianMCPMod.Rd +++ b/man/performBayesianMCPMod.Rd @@ -21,3 +21,26 @@ bmcpmod test result as well as modelling result. \description{ performs bayesian MCP Test step and modelling in a combined fashion. See performBayesianMCP function for MCT Test step and getModelFits for the modelling step } +\examples{ +# example code +models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +doses = c(0, 0.5, 2,4, 8)) +dose_levels=c(0, 0.5, 2, 4, 8) +sd_posterior = c(2.8,3,2.5,3.5,4) +contr_mat<- getContr( +mods = models, +dose_levels = dose_levels, +sd_posterior = sd_posterior) +critVal<- getCritProb( + mods = models, + dose_weights =c(50,50,50,50,50), #reflecting the planned sample size + dose_levels = dose_levels, + alpha_crit_val = 0.05) +posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), + DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), + 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)) +performBayesianMCPMod(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal,simple = FALSE) + +} diff --git a/man/plot.modelFits.Rd b/man/plot.modelFits.Rd index be11a8a..8ec1128 100644 --- a/man/plot.modelFits.Rd +++ b/man/plot.modelFits.Rd @@ -44,6 +44,20 @@ plts returns a ggplot2 object \description{ plot function based on the ggplot2 package. Providing visualizations for each model and a average Fit. Black lines show the fitted dose response models and an AIC based average model. Dots indicate the posterior median and vertical lines show corresponding credible intervals (i.e. the variability of the posterior distribution of the respective dose group). -To assess the uncertainty of the model fit one can in addition visualize credible bands (yellow/orange shaded areas). The calculation of these bands is performed via the getBootstrapQuantiles function. +To assess the uncertainty of the model fit one can in addition visualize credible bands (default coloring as orange shaded areas). The calculation of these bands is performed via the getBootstrapQuantiles function. The default setting is that these credible bands are not calculated. } +\examples{ +# example code +posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2), + DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2), + DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) , + DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) , + DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2)) +models=c("emax","exponential","sigEmax","linear") +dose_levels=c(0,1,2,4,8) +fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels) +fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE) +plot(fit, cr_bands = TRUE) +plot(fit_simple, cr_bands = TRUE,alpha_CrB =c(0.05,0.1,0.5)) +} From da6e8524800c220be3267c0c332077918d07bbb3 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Mon, 18 Dec 2023 21:00:37 +0000 Subject: [PATCH 25/39] References added to the documentation. --- R/BMCPMod.R | 4 ++-- R/bootstrapping.R | 4 ++-- R/modelling.R | 3 ++- man/getBootstrapQuantiles.Rd | 5 ++++- man/getModelFits.Rd | 5 ++++- man/performBayesianMCP.Rd | 5 ++++- 6 files changed, 18 insertions(+), 8 deletions(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 237352c..f43a195 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -380,10 +380,10 @@ addSignificance <- function ( #' @title performBayesianMCP #' -#' @description performs bayesian MCP Test step, as described in Fleischer et al. (Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670.) +#' @description performs bayesian MCP Test step, as described in [Fleischer et al.,2022]. #' Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate the posterior probability is calculated that the contrast is bigger than 0 (based on the posterior distribution of the dose groups). #' In order to obtain significant test decision we consider the maximum of the posterior probabilities across the different models. This maximum is compared with a (multiplicity adjusted) critical value (on the probability scale). -#' +#' @references Fleischer F, Bossert S, Deng Q, Loley C, Gierse J. Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670. doi:10.1002/pst.2193 #' @param posterior_list a getPosterior object with information about the (mixture) posterior distribution per dose group #' @param contr a getContrMat object, contrast matrix to be used for the testing step. #' @param crit_prob_adj a getCritProb object, specifying the critical value to be used for the testing (on the probability scale) diff --git a/R/bootstrapping.R b/R/bootstrapping.R index 1cb306d..16764c2 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -4,9 +4,9 @@ #' Hereby the credible intervals are calculated as follows. #' Samples from the posterior distribution are drawn (via the RBesT function rmix) and for every sample the simplified fitting step (see getModelFits function) and a prediction is performed. #' These fits are then used to identify the specified quantiles. -#' This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in O'Quigley, Iasonos, J., and B. Bornkamp. 2017. Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials. Boca Raton: CRC Press. +#' This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in [O'Quigley et al., 2017]. #' Instead of drawing n bootstrap samples from the sampling distribution of the trial dose-response estimates, here the samples are directly taken from the posterior distribution. -#' +#' @references O'Quigley, J., Iasonos, A., & Bornkamp, B. (Eds.). (2017). Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials (1st ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9781315151984 #' @param model_fits an object of class modelFits, i.e. information about fitted models & corresponding model coefficients as well as the posterior distribution that was the basis for the model fitting #' @param quantiles a vector of quantiles that should be evaluated #' @param n_samples number of samples that should be drawn as basis for the diff --git a/R/modelling.R b/R/modelling.R index 013dec2..2022c91 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -12,11 +12,12 @@ #' \deqn{gAIC_{m} = \sum_{l=1}^{L} \tilde{\omega}_{l} \sum_{i=0}^{K} \frac{1}{\Sigma_{l_{i,i}}} (\theta_{l_i}^Y - f(dose_{i},\hat{\theta}_{m}))^2 + 2p } #' where \eqn{p} denotes the number of estimated parameters and \eqn{K} the number of active dose levels. #' Here as well for the simplified case the formula reduces to one summand as $L=1$. -#' Corresponding gAIC based weights for model $M$ are calculated as outlined in Schorning et al. 2016: +#' Corresponding gAIC based weights for model $M$ are calculated as outlined in [Schorning et al., 2016] #' \deqn{ #' \Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_(m))} #' } #' where \eqn{Q} denotes the number of models included in the averaging procedure. +#' @references Schorning K, Bornkamp B, Bretz F, Dette H. 2016. “Model selection versus model averaging in dose finding studies”. Stat Med; 35; 4021-4040. #' @param models list of model names for which a fit will be performed. #' @param dose_levels a vector containing the different dosage levels. #' @param posterior a getPosterior object, containing the (multivariate) posterior distribution per dosage level. diff --git a/man/getBootstrapQuantiles.Rd b/man/getBootstrapQuantiles.Rd index cc0d1a4..764d0de 100644 --- a/man/getBootstrapQuantiles.Rd +++ b/man/getBootstrapQuantiles.Rd @@ -31,7 +31,7 @@ A function for the calculation of credible intervals to assess the uncertainty f Hereby the credible intervals are calculated as follows. Samples from the posterior distribution are drawn (via the RBesT function rmix) and for every sample the simplified fitting step (see getModelFits function) and a prediction is performed. These fits are then used to identify the specified quantiles. -This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in O'Quigley, Iasonos, J., and B. Bornkamp. 2017. Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials. Boca Raton: CRC Press. +This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in \link{O'Quigley et al., 2017}. Instead of drawing n bootstrap samples from the sampling distribution of the trial dose-response estimates, here the samples are directly taken from the posterior distribution. } \examples{ @@ -48,3 +48,6 @@ fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dos getBootstrapQuantiles(fit, quantiles = c(0.025,0.5, 0.975), doses = c(0, 1,2,3,4,6,8)) getBootstrapQuantiles(fit_simple, n_samples=2000, quantiles = c(0.025,0.5, 0.975), doses = c(0, 1,2,3,4,6,8)) } +\references{ +O'Quigley, J., Iasonos, A., & Bornkamp, B. (Eds.). (2017). Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials (1st ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9781315151984 +} diff --git a/man/getModelFits.Rd b/man/getModelFits.Rd index 715feed..a2d1327 100644 --- a/man/getModelFits.Rd +++ b/man/getModelFits.Rd @@ -31,7 +31,7 @@ The generalized AIC values are calculated via the formula \deqn{gAIC_{m} = \sum_{l=1}^{L} \tilde{\omega}_{l} \sum_{i=0}^{K} \frac{1}{\Sigma_{l_{i,i}}} (\theta_{l_i}^Y - f(dose_{i},\hat{\theta}_{m}))^2 + 2p } where \eqn{p} denotes the number of estimated parameters and \eqn{K} the number of active dose levels. Here as well for the simplified case the formula reduces to one summand as $L=1$. -Corresponding gAIC based weights for model $M$ are calculated as outlined in Schorning et al. 2016: +Corresponding gAIC based weights for model $M$ are calculated as outlined in \link{Schorning et al., 2016} \deqn{ \Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_(m))} } @@ -49,3 +49,6 @@ dose_levels=c(0,1,2,4,8) fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels) fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE) } +\references{ +Schorning K, Bornkamp B, Bretz F, Dette H. 2016. “Model selection versus model averaging in dose finding studies”. Stat Med; 35; 4021-4040. +} diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd index fa89d63..4e950e1 100644 --- a/man/performBayesianMCP.Rd +++ b/man/performBayesianMCP.Rd @@ -17,7 +17,7 @@ performBayesianMCP(posterior_list, contr, crit_prob_adj) b_mcp test result, with information about p-values for the individual dose-response shapes and overall significance } \description{ -performs bayesian MCP Test step, as described in Fleischer et al. (Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670.) +performs bayesian MCP Test step, as described in \link{Fleischer et al.,2022}. Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate the posterior probability is calculated that the contrast is bigger than 0 (based on the posterior distribution of the dose groups). In order to obtain significant test decision we consider the maximum of the posterior probabilities across the different models. This maximum is compared with a (multiplicity adjusted) critical value (on the probability scale). } @@ -44,3 +44,6 @@ posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigm performBayesianMCP(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal) } +\references{ +Fleischer F, Bossert S, Deng Q, Loley C, Gierse J. Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670. doi:10.1002/pst.2193 +} From 69d65697c63d4e9038c331cd63027cbed9f4fa15 Mon Sep 17 00:00:00 2001 From: wojcieko Date: Mon, 18 Dec 2023 23:13:40 +0100 Subject: [PATCH 26/39] - beautification --- R/s3methods.R | 30 +++++++++++++++++------------- vignettes/analysis_normal.Rmd | 4 ++-- 2 files changed, 19 insertions(+), 15 deletions(-) diff --git a/R/s3methods.R b/R/s3methods.R index 2583ca5..33390d9 100644 --- a/R/s3methods.R +++ b/R/s3methods.R @@ -77,22 +77,26 @@ print.BayesianMCP <- function ( } ## ModelFits ---------------------------------------------- -#' @description This function performs model predictions based on the provided model and dose specifications +#' @description This function performs model predictions based on the provided +#' model and dose specifications #' -#' @param object a modelFits object containing information about the fitted model coefficients -#' @param doses a vector specifying the doses for which a prediction should be done getContrMat object, contrast matrix to be used for the testing step. - +#' @param object a modelFits object containing information about the fitted +#' model coefficients +#' @param doses a vector specifying the doses for which a prediction should be +#' done getContrMat object, contrast matrix to be used for the testing step. #' @examples #' # example code -#' posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2), -#' DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2), -#' DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) , -#' DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) , -#' DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2)) -#' models=c("emax","exponential","sigEmax","linear") -#' dose_levels=c(0,1,2,4,8) -#' fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels) -#' predict(fit, doses = c(0, 1,3,4,6,8)) +#' posterior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2), +#' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2), +#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) , +#' DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) , +#' DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1), sigma = 2)) +#' models <- c("emax", "exponential", "sigEmax", "linear") +#' dose_levels <- c(0, 1, 2, 4, 8) +#' fit <- getModelFits(models = models, +#' posterior = posterior_list, +#' dose_levels = dose_levels) +#' predict(fit, doses = c(0, 1, 3, 4, 6, 8)) #' #' @return a list with the model predictions for the specified models and doses #' @export diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 0107f6c..6da9c1a 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -289,8 +289,8 @@ It is also possible to perform the testing and modelling step in a combined fash ```{r} overall_result<-performBayesianMCPMod( - posterior_list = posterior, - contr = contr_mat, + posterior_list = posterior, + contr = contr_mat, crit_prob_adj = crit_pval, simple = FALSE) From 9a4ce45d59f7dc4b203faa45d429722583a4b955 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Tue, 19 Dec 2023 09:30:28 +0000 Subject: [PATCH 27/39] Bug fixes around BMCTTest functions --- R/BMCPMod.R | 11 +++++------ R/posterior.R | 8 ++++++-- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index f43a195..81b0eb9 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -287,7 +287,10 @@ getCritProb <- function ( #' 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)) -#' performBayesianMCPMod(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal,simple = FALSE) +#' performBayesianMCPMod(posterior_list=posterior_list, +#' contr=contr_mat, +#' crit_prob_adj=critVal, +#' simple = FALSE) #' #' @return bmcpmod test result as well as modelling result. #' @@ -421,11 +424,7 @@ performBayesianMCP <- function( ) { - if (inherits(posterior_list, "postList")) { - - posterior_list <- list(posterior_list) - - } + posterior_list <- list(posterior_list) if (inherits(contr, "optContr")) { diff --git a/R/posterior.R b/R/posterior.R index b79ec08..6693fbb 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -20,10 +20,13 @@ #' DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13) ,sigma = 2)) #' mu<-c(0,1,1.5,2,2.5) #' se<-c(5,4,6,7,8) -#' getPosterior <- function( +#' posterior_list <- getPosterior( #' prior_list = prior_list, #' mu_hat = mu, #' se_hat = se) +#' +#' summary(posterior_list) +#' #' @export getPosterior <- function( prior_list, @@ -92,7 +95,8 @@ getPosteriorI <- function( } - post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = se_hat) + post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = se_hat, + SIMPLIFY = FALSE) if (is.null(names(prior_list))) { From d49821a77e7fce905b749ab74ef4f560ed7f7636 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Tue, 19 Dec 2023 10:28:25 +0000 Subject: [PATCH 28/39] Further updates around documentation and simulation vignette. --- R/BMCPMod.R | 60 +++++++++++++++++++++++++------- R/simulation.R | 8 +++-- vignettes/Simulation_Example.Rmd | 40 ++++++++++++++------- 3 files changed, 80 insertions(+), 28 deletions(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 81b0eb9..829bdb5 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -5,14 +5,36 @@ #' @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 tbd +#' @param sd a positive value, specification of assumed sd #' @param n_sim number of simulations to be performed -#' @param alpha_crit_val critical value to be used for the testing (on the probability scale) +#' @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 tbd Default FALSE -#' @param contr tbd Default NULL -#' @param dr_means tbd Default NULL +#' @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 +#' +#' @examples +#' # example code +#' models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +#' doses = c(0, 0.5, 2,4, 8),maxEff= 6) +#' sd = 12 +#' prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2), +#' DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), +#' 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, +#' mods = mods, +#' prior_list = prior_list, +#' sd = sd) +#' +#' success_probabilities + #' @export assessDesign <- function ( @@ -282,15 +304,18 @@ getCritProb <- function ( #' dose_weights =c(50,50,50,50,50), #reflecting the planned sample size #' dose_levels = dose_levels, #' alpha_crit_val = 0.05) -#' posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), +#' prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), #' DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), #' 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)) -#' performBayesianMCPMod(posterior_list=posterior_list, -#' contr=contr_mat, -#' crit_prob_adj=critVal, -#' simple = FALSE) +#' mu<-c(0,1,1.5,2,2.5) +#' se<-c(5,4,6,7,8) +#' posterior_list <- getPosterior( +#' prior_list = prior_list, +#' mu_hat = mu, +#' se_hat = se) +#' performBayesianMCPMod(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal,simple = FALSE) #' #' @return bmcpmod test result as well as modelling result. #' @@ -406,11 +431,17 @@ addSignificance <- function ( #' dose_weights =c(50,50,50,50,50), #reflecting the planned sample size #' dose_levels = dose_levels, #' alpha_crit_val = 0.05) -#' posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), +#' prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), #' DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), #' 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)) +#' mu<-c(0,1,1.5,2,2.5) +#' se<-c(5,4,6,7,8) +#' posterior_list <- getPosterior( +#' prior_list = prior_list, +#' mu_hat = mu, +#' se_hat = se) #' performBayesianMCP(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal) #' #' @return b_mcp test result, with information about p-values for the individual dose-response shapes and overall significance @@ -423,8 +454,11 @@ performBayesianMCP <- function( crit_prob_adj ) { - - posterior_list <- list(posterior_list) + if (inherits(posterior_list, "postList")) { + + posterior_list <- list(posterior_list) + + } if (inherits(contr, "optContr")) { diff --git a/R/simulation.R b/R/simulation.R index 65330cd..48561cf 100644 --- a/R/simulation.R +++ b/R/simulation.R @@ -1,8 +1,12 @@ #' @title simulateData #' +#' @description +#' Function to simulate patient level data for a normally distributed endpoint +#' +#' #' @param n_patients vector containing number of patients as a numerical #' value per dose-group. -#' @param dose_levels vector containing the different doseage levels. +#' @param dose_levels vector containing the different dosage levels. #' @param sd standard deviation on patient level. #' @param mods An object of class "Mods" as specified in the Dosefinding package. #' @param n_sim number of simulations to be performed, @@ -10,7 +14,7 @@ #' @param true_model Default value is NULL. #' Assumed true underlying model. Provided via a String. e.g. "emax". #' In case of NULL, all dose-response models, included in the mods input parameter will be used. -#' @param dr_means tbd. Default NULL. +#' @param dr_means a vector, with information about assumed effects per dose group. Default NULL. #' #' @return sim_data one list object, containing patient level simulated data for all assumed true models. #' Also providing information about simulation iteration, patient number as well as dosage levels. diff --git a/vignettes/Simulation_Example.Rmd b/vignettes/Simulation_Example.Rmd index e374f3d..619b7c6 100644 --- a/vignettes/Simulation_Example.Rmd +++ b/vignettes/Simulation_Example.Rmd @@ -26,12 +26,12 @@ 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 to other vignette] we focus on the indication MDD and make use of historical data that is included in the clinDR package. +As in [the analysis example vignette](analysis_normal.Rmd) 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 -In a first step a meta analytic prior will be calculated using historical data from 5 trials (with main endpoint CfB in MADRS score after 8 weeks). -Please note that only information from the control group will be integrated (leading to an informative multicomponent prior for the control group), while for the active groups a non-informative prior will be specified. +In a first step a meta analytic prior will be calculated using historical data from 5 trials (with main endpoint Change from baseline in MADRS score after 8 weeks). +Please note that only information from the control group will be integrated (leading to an informative mixture prior for the control group), while for the active groups a non-informative prior will be specified. ```{r Historical Data} @@ -51,7 +51,9 @@ hist_data <- data.frame( sd_tot <- with(hist_data, sum(sd * n) / sum(n)) ``` -```{r Calculation of a MAP prior} +We will make use of the same getPriorList function as in the [analysis example vignette](analysis_normal.Rmd) to create a MAP prior. + +```{r Calculation of a MAP prior, echo = FALSE} getPriorList <- function ( hist_data, @@ -103,19 +105,24 @@ getPriorList <- function ( return (prior_list) } +``` + +```{r Evaluation of a MAP prior} dose_levels <- c(0, 2.5, 5, 10, 20) prior_list <- getPriorList( hist_data = hist_data, - dose_levels = dose_levels) + dose_levels = dose_levels, + robust_weight = 0.3) +getESS(prior_list) ``` # Specification of new trial design -For the hypothetical new trial we plan with 4 active dose levels $2.5, 5, 10, 20$ and we specify a broad set of potential dose-response relationships, including a linear, an exponential, an emax and a sigEMax model. -Furthermore we assume a maximum effect of -3 on top of control (i.e. assuming that active treatment can reduce the MADRS score after 8 weeks by up to 15.8) and plan a trial with 80 patients for all active groups and 60 patients for control. +For the hypothetical new trial we plan with 4 active dose levels \eqn{2.5, 5, 10, 20} and we specify a broad set of potential dose-response relationships, including a linear, an exponential, an emax and a sigEMax model. +Furthermore we assume a maximum effect of -3 on top of control (i.e. assuming that active treatment can reduce the MADRS score after 8 weeks by up to 15.8) and plan a trial with 80 patients for all active groups and 60 patients for control. ```{r} #Pre-Specification (B)MCPMod @@ -134,12 +141,7 @@ sigemax <- DoseFinding::guesst(d = c(2.5, 5), sigemax2 <- DoseFinding::guesst(d = c(2, 4), p = c(0.3, 0.8), model = "sigEmax") -# beta <- DoseFinding::guesst(d = 5, -# p = 0.8, -# model = "betaMod", -# dMax = 1, -# scal = 1.2, -# Maxd = 20) + mods <- DoseFinding::Mods( linear = NULL, @@ -181,4 +183,16 @@ success_probabilities_uneq For this specific trial setting the adapted allocation ratio leads to increased success probabilities under all assumed dose response relationships. +Instead of specifying the assumed effects via the models it is also possible to directly specify the effects for the individual dose levels via the dr_means input. This allows e.g. also the simulation of scenarios with a prior-data conflict. + +```{r} +success_probabilities <- assessDesign( + n_patients = c(60, 80, 80, 80, 80), + mods = mods, + prior_list = prior_list, + sd = sd_tot, + dr_means = c(-12,-14,-15,-16,-17)) +success_probabilities +``` + From 1f76915d0397775e536d0f4af2ec04edc5f17c96 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Tue, 19 Dec 2023 10:57:46 +0000 Subject: [PATCH 29/39] Further updates (to fix problems with links and documentation). --- R/BMCPMod.R | 2 +- R/bootstrapping.R | 2 +- R/modelling.R | 8 +++---- R/s3methods.R | 1 + man/assessDesign.Rd | 33 +++++++++++++++++++++++----- man/getBootstrapQuantiles.Rd | 2 +- man/getModelFits.Rd | 8 +++---- man/getPosterior.Rd | 5 ++++- man/performBayesianMCP.Rd | 10 +++++++-- man/performBayesianMCPMod.Rd | 8 ++++++- man/predict.modelFits.Rd | 37 ++++++++++++++++++++++++++++++++ man/simulateData.Rd | 6 +++--- vignettes/Simulation_Example.Rmd | 4 ++-- 13 files changed, 101 insertions(+), 25 deletions(-) create mode 100644 man/predict.modelFits.Rd diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 829bdb5..be8ccb4 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -408,7 +408,7 @@ addSignificance <- function ( #' @title performBayesianMCP #' -#' @description performs bayesian MCP Test step, as described in [Fleischer et al.,2022]. +#' @description performs bayesian MCP Test step, as described in Fleischer et al. (2022). #' Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate the posterior probability is calculated that the contrast is bigger than 0 (based on the posterior distribution of the dose groups). #' In order to obtain significant test decision we consider the maximum of the posterior probabilities across the different models. This maximum is compared with a (multiplicity adjusted) critical value (on the probability scale). #' @references Fleischer F, Bossert S, Deng Q, Loley C, Gierse J. Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670. doi:10.1002/pst.2193 diff --git a/R/bootstrapping.R b/R/bootstrapping.R index 16764c2..7570a8e 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -4,7 +4,7 @@ #' Hereby the credible intervals are calculated as follows. #' Samples from the posterior distribution are drawn (via the RBesT function rmix) and for every sample the simplified fitting step (see getModelFits function) and a prediction is performed. #' These fits are then used to identify the specified quantiles. -#' This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in [O'Quigley et al., 2017]. +#' This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in O'Quigley et al. (2017). #' Instead of drawing n bootstrap samples from the sampling distribution of the trial dose-response estimates, here the samples are directly taken from the posterior distribution. #' @references O'Quigley, J., Iasonos, A., & Bornkamp, B. (Eds.). (2017). Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials (1st ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9781315151984 #' @param model_fits an object of class modelFits, i.e. information about fitted models & corresponding model coefficients as well as the posterior distribution that was the basis for the model fitting diff --git a/R/modelling.R b/R/modelling.R index 2022c91..1c42024 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -3,16 +3,16 @@ #' @description Fits dose-response curves for the specified dose-response models, based on the posterior distributions. #' For the simplified fit, multivariate normal distributions will be approximated and reduced by one-dimensional normal distributions. #' For the default case, the Nelder-Mead algorithm is used. -#' In detail, for both approaches the mean vector \eqn{\theta^{Y}} and the covariance \eqn{\Sigma} of the (mixture) posterior distributions and the corresponding posterior weights \eqn{\tilde{\omega}_{l}} for $l \in {1,...,L}$ are used as basis +#' In detail, for both approaches the mean vector \eqn{\theta^{Y}} and the covariance \eqn{\Sigma} of the (mixture) posterior distributions and the corresponding posterior weights \eqn{\tilde{\omega}_{l}} for \eqn{l \in {1,...,L}} are used as basis #' For the full fit a GLS estimator is used to minimize the following expression for the respective dose-response models \eqn{m} #' \deqn{ \hat{\theta}_{m}=argmin_{\theta_{m}} \sum_{i=l}^{L} \tilde{\omega}_{l}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))'\Sigma_{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))} #' Herefore the function nloptr of the nloptr package is utilized. -#' In the simplified case $L=1$, as the dimension of the posterior is reduced to 1 first. +#' In the simplified case \eqn{L=1}, as the dimension of the posterior is reduced to 1 first. #' The generalized AIC values are calculated via the formula #' \deqn{gAIC_{m} = \sum_{l=1}^{L} \tilde{\omega}_{l} \sum_{i=0}^{K} \frac{1}{\Sigma_{l_{i,i}}} (\theta_{l_i}^Y - f(dose_{i},\hat{\theta}_{m}))^2 + 2p } #' where \eqn{p} denotes the number of estimated parameters and \eqn{K} the number of active dose levels. -#' Here as well for the simplified case the formula reduces to one summand as $L=1$. -#' Corresponding gAIC based weights for model $M$ are calculated as outlined in [Schorning et al., 2016] +#' Here as well for the simplified case the formula reduces to one summand as \eqn{L=1}. +#' Corresponding gAIC based weights for model $M$ are calculated as outlined in Schorning et al. (2016) #' \deqn{ #' \Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_(m))} #' } diff --git a/R/s3methods.R b/R/s3methods.R index 33390d9..6fd93ce 100644 --- a/R/s3methods.R +++ b/R/s3methods.R @@ -77,6 +77,7 @@ print.BayesianMCP <- function ( } ## ModelFits ---------------------------------------------- +#' @title predict.modelFits #' @description This function performs model predictions based on the provided #' model and dose specifications #' diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd index 91e3392..61f9476 100644 --- a/man/assessDesign.Rd +++ b/man/assessDesign.Rd @@ -25,20 +25,43 @@ assessDesign( \item{prior_list}{a prior_list object specifying the utilized prior for the different dose groups} -\item{sd}{tbd} +\item{sd}{a positive value, specification of assumed sd} \item{n_sim}{number of simulations to be performed} -\item{alpha_crit_val}{critical value to be used for the testing (on the probability scale)} +\item{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.} \item{simple}{boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.} -\item{reestimate}{tbd Default FALSE} +\item{reestimate}{boolean variable, defining whether critical value should be calculated with re-estimated contrasts (see getCritProb function for more details). Default FALSE} -\item{contr}{tbd Default NULL} +\item{contr}{Allows specification of a fixed contrasts matrix. Default NULL} -\item{dr_means}{tbd Default NULL} +\item{dr_means}{a vector, allows specification of individual (not model based) assumed effects per dose group. Default NULL} +} +\value{ +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 } \description{ This function performs simulation based trial design evaluations for a set of specified dose-response models } +\examples{ +# example code + models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +doses = c(0, 0.5, 2,4, 8),maxEff= 6) +sd = 12 +prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2), + DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), + 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, +mods = mods, +prior_list = prior_list, +sd = sd) + +success_probabilities +} diff --git a/man/getBootstrapQuantiles.Rd b/man/getBootstrapQuantiles.Rd index 764d0de..4e63309 100644 --- a/man/getBootstrapQuantiles.Rd +++ b/man/getBootstrapQuantiles.Rd @@ -31,7 +31,7 @@ A function for the calculation of credible intervals to assess the uncertainty f Hereby the credible intervals are calculated as follows. Samples from the posterior distribution are drawn (via the RBesT function rmix) and for every sample the simplified fitting step (see getModelFits function) and a prediction is performed. These fits are then used to identify the specified quantiles. -This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in \link{O'Quigley et al., 2017}. +This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in O'Quigley et al. (2017). Instead of drawing n bootstrap samples from the sampling distribution of the trial dose-response estimates, here the samples are directly taken from the posterior distribution. } \examples{ diff --git a/man/getModelFits.Rd b/man/getModelFits.Rd index a2d1327..a20a053 100644 --- a/man/getModelFits.Rd +++ b/man/getModelFits.Rd @@ -22,16 +22,16 @@ model_fits returns a list, containing information about the fitted model coeffic Fits dose-response curves for the specified dose-response models, based on the posterior distributions. For the simplified fit, multivariate normal distributions will be approximated and reduced by one-dimensional normal distributions. For the default case, the Nelder-Mead algorithm is used. -In detail, for both approaches the mean vector \eqn{\theta^{Y}} and the covariance \eqn{\Sigma} of the (mixture) posterior distributions and the corresponding posterior weights \eqn{\tilde{\omega}_{l}} for $l \in {1,...,L}$ are used as basis +In detail, for both approaches the mean vector \eqn{\theta^{Y}} and the covariance \eqn{\Sigma} of the (mixture) posterior distributions and the corresponding posterior weights \eqn{\tilde{\omega}_{l}} for \eqn{l \in {1,...,L}} are used as basis For the full fit a GLS estimator is used to minimize the following expression for the respective dose-response models \eqn{m} \deqn{ \hat{\theta}_{m}=argmin_{\theta_{m}} \sum_{i=l}^{L} \tilde{\omega}_{l}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))'\Sigma_{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))} Herefore the function nloptr of the nloptr package is utilized. -In the simplified case $L=1$, as the dimension of the posterior is reduced to 1 first. +In the simplified case \eqn{L=1}, as the dimension of the posterior is reduced to 1 first. The generalized AIC values are calculated via the formula \deqn{gAIC_{m} = \sum_{l=1}^{L} \tilde{\omega}_{l} \sum_{i=0}^{K} \frac{1}{\Sigma_{l_{i,i}}} (\theta_{l_i}^Y - f(dose_{i},\hat{\theta}_{m}))^2 + 2p } where \eqn{p} denotes the number of estimated parameters and \eqn{K} the number of active dose levels. -Here as well for the simplified case the formula reduces to one summand as $L=1$. -Corresponding gAIC based weights for model $M$ are calculated as outlined in \link{Schorning et al., 2016} +Here as well for the simplified case the formula reduces to one summand as \eqn{L=1}. +Corresponding gAIC based weights for model $M$ are calculated as outlined in Schorning et al. (2016) \deqn{ \Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_(m))} } diff --git a/man/getPosterior.Rd b/man/getPosterior.Rd index 005b6d0..e30d5d2 100644 --- a/man/getPosterior.Rd +++ b/man/getPosterior.Rd @@ -40,8 +40,11 @@ prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13) ,sigma = 2)) mu<-c(0,1,1.5,2,2.5) se<-c(5,4,6,7,8) -getPosterior <- function( +posterior_list <- getPosterior( prior_list = prior_list, mu_hat = mu, se_hat = se) + +summary(posterior_list) + } diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd index 4e950e1..3fec02c 100644 --- a/man/performBayesianMCP.Rd +++ b/man/performBayesianMCP.Rd @@ -17,7 +17,7 @@ performBayesianMCP(posterior_list, contr, crit_prob_adj) b_mcp test result, with information about p-values for the individual dose-response shapes and overall significance } \description{ -performs bayesian MCP Test step, as described in \link{Fleischer et al.,2022}. +performs bayesian MCP Test step, as described in Fleischer et al. (2022). Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate the posterior probability is calculated that the contrast is bigger than 0 (based on the posterior distribution of the dose groups). In order to obtain significant test decision we consider the maximum of the posterior probabilities across the different models. This maximum is compared with a (multiplicity adjusted) critical value (on the probability scale). } @@ -36,11 +36,17 @@ critVal<- getCritProb( dose_weights =c(50,50,50,50,50), #reflecting the planned sample size dose_levels = dose_levels, alpha_crit_val = 0.05) -posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), +prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), 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)) +mu<-c(0,1,1.5,2,2.5) +se<-c(5,4,6,7,8) +posterior_list <- getPosterior( + prior_list = prior_list, + mu_hat = mu, + se_hat = se) performBayesianMCP(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal) } diff --git a/man/performBayesianMCPMod.Rd b/man/performBayesianMCPMod.Rd index e43ec82..e02ffd1 100644 --- a/man/performBayesianMCPMod.Rd +++ b/man/performBayesianMCPMod.Rd @@ -36,11 +36,17 @@ critVal<- getCritProb( dose_weights =c(50,50,50,50,50), #reflecting the planned sample size dose_levels = dose_levels, alpha_crit_val = 0.05) -posterior_list = list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), +prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), 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)) +mu<-c(0,1,1.5,2,2.5) +se<-c(5,4,6,7,8) +posterior_list <- getPosterior( + prior_list = prior_list, + mu_hat = mu, + se_hat = se) performBayesianMCPMod(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal,simple = FALSE) } diff --git a/man/predict.modelFits.Rd b/man/predict.modelFits.Rd new file mode 100644 index 0000000..7fecd10 --- /dev/null +++ b/man/predict.modelFits.Rd @@ -0,0 +1,37 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/s3methods.R +\name{predict.modelFits} +\alias{predict.modelFits} +\title{predict.modelFits} +\usage{ +\method{predict}{modelFits}(object, doses = NULL, ...) +} +\arguments{ +\item{object}{a modelFits object containing information about the fitted +model coefficients} + +\item{doses}{a vector specifying the doses for which a prediction should be +done getContrMat object, contrast matrix to be used for the testing step.} +} +\value{ +a list with the model predictions for the specified models and doses +} +\description{ +This function performs model predictions based on the provided +model and dose specifications +} +\examples{ +# example code +posterior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2), + DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2), + DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) , + DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) , + DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1), sigma = 2)) +models <- c("emax", "exponential", "sigEmax", "linear") +dose_levels <- c(0, 1, 2, 4, 8) +fit <- getModelFits(models = models, + posterior = posterior_list, + dose_levels = dose_levels) +predict(fit, doses = c(0, 1, 3, 4, 6, 8)) + +} diff --git a/man/simulateData.Rd b/man/simulateData.Rd index 1b55e08..1367b3a 100644 --- a/man/simulateData.Rd +++ b/man/simulateData.Rd @@ -18,7 +18,7 @@ simulateData( \item{n_patients}{vector containing number of patients as a numerical value per dose-group.} -\item{dose_levels}{vector containing the different doseage levels.} +\item{dose_levels}{vector containing the different dosage levels.} \item{sd}{standard deviation on patient level.} @@ -31,12 +31,12 @@ Default is 1000} Assumed true underlying model. Provided via a String. e.g. "emax". In case of NULL, all dose-response models, included in the mods input parameter will be used.} -\item{dr_means}{tbd. Default NULL.} +\item{dr_means}{a vector, with information about assumed effects per dose group. Default NULL.} } \value{ sim_data one list object, containing patient level simulated data for all assumed true models. Also providing information about simulation iteration, patient number as well as dosage levels. } \description{ -simulateData +Function to simulate patient level data for a normally distributed endpoint } diff --git a/vignettes/Simulation_Example.Rmd b/vignettes/Simulation_Example.Rmd index 619b7c6..8f6c247 100644 --- a/vignettes/Simulation_Example.Rmd +++ b/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 [the analysis example vignette](analysis_normal.Rmd) we focus on the indication MDD and make use of historical data that is included in the clinDR package. +As in [the analysis example vignette](analysis_normal.html) 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 @@ -51,7 +51,7 @@ hist_data <- data.frame( sd_tot <- with(hist_data, sum(sd * n) / sum(n)) ``` -We will make use of the same getPriorList function as in the [analysis example vignette](analysis_normal.Rmd) to create a MAP prior. +We will make use of the same getPriorList function as in the [analysis example vignette](analysis_normal.html) to create a MAP prior. ```{r Calculation of a MAP prior, echo = FALSE} getPriorList <- function ( From fdca49e3d66a58f8e8d3fcf026c1ca7fb1891917 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Tue, 19 Dec 2023 13:00:16 +0000 Subject: [PATCH 30/39] Bug fix for assess design (to have correct display for cases where individual dr vector is defined). --- R/BMCPMod.R | 18 +++++++++++------- R/simulation.R | 7 +++++++ 2 files changed, 18 insertions(+), 7 deletions(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index be8ccb4..37ea2a1 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -20,11 +20,11 @@ #' models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, #' doses = c(0, 0.5, 2,4, 8),maxEff= 6) #' sd = 12 -#' prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2), -#' DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), -#' 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)) +#' prior_list<-list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2), +#' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), +#' 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( @@ -114,8 +114,12 @@ assessDesign <- function ( names(eval_design) <- model_names attr(eval_design, "avgSuccessRate") <- avg_success_rate - attr(eval_design, "placEff") <- attr(mods, "placEff") - attr(eval_design, "maxEff") <- attr(mods, "maxEff") + attr(eval_design, "placEff") <- ifelse(test = is.null(dr_means), + yes = attr(mods, "placEff"), + no = dr_means[1]) + attr(eval_design, "maxEff") <- ifelse(test = is.null(dr_means), + yes = attr(mods, "maxEff"), + no = diff(range(dr_means))) attr(eval_design, "sampleSize") <- n_patients attr(eval_design, "priorESS") <- round(getESS(prior_list), 1) diff --git a/R/simulation.R b/R/simulation.R index 48561cf..09794e9 100644 --- a/R/simulation.R +++ b/R/simulation.R @@ -72,6 +72,13 @@ simulateData <- function( } + if (!is.null(dr_means)) { + + sim_data <- sim_data[1:4] + colnames(sim_data) <- c(colnames(sim_data)[-4], "dose_resp") + + } + return (sim_data) } From 62aceacfaa63d32c4faea957513f1b996ae27309 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Tue, 19 Dec 2023 14:33:54 +0000 Subject: [PATCH 31/39] Additional updates around documentation. --- R/modelling.R | 2 +- man/assessDesign.Rd | 10 +++++----- man/getModelFits.Rd | 2 +- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/R/modelling.R b/R/modelling.R index 1c42024..3ef5daf 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -12,7 +12,7 @@ #' \deqn{gAIC_{m} = \sum_{l=1}^{L} \tilde{\omega}_{l} \sum_{i=0}^{K} \frac{1}{\Sigma_{l_{i,i}}} (\theta_{l_i}^Y - f(dose_{i},\hat{\theta}_{m}))^2 + 2p } #' where \eqn{p} denotes the number of estimated parameters and \eqn{K} the number of active dose levels. #' Here as well for the simplified case the formula reduces to one summand as \eqn{L=1}. -#' Corresponding gAIC based weights for model $M$ are calculated as outlined in Schorning et al. (2016) +#' Corresponding gAIC based weights for model \eqn{M} are calculated as outlined in Schorning et al. (2016) #' \deqn{ #' \Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_(m))} #' } diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd index 61f9476..5b11892 100644 --- a/man/assessDesign.Rd +++ b/man/assessDesign.Rd @@ -50,11 +50,11 @@ This function performs simulation based trial design evaluations for a set of sp models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8),maxEff= 6) sd = 12 -prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2), - DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), - 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)) +prior_list<-list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2), + DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), + 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( diff --git a/man/getModelFits.Rd b/man/getModelFits.Rd index a20a053..fb49d23 100644 --- a/man/getModelFits.Rd +++ b/man/getModelFits.Rd @@ -31,7 +31,7 @@ The generalized AIC values are calculated via the formula \deqn{gAIC_{m} = \sum_{l=1}^{L} \tilde{\omega}_{l} \sum_{i=0}^{K} \frac{1}{\Sigma_{l_{i,i}}} (\theta_{l_i}^Y - f(dose_{i},\hat{\theta}_{m}))^2 + 2p } where \eqn{p} denotes the number of estimated parameters and \eqn{K} the number of active dose levels. Here as well for the simplified case the formula reduces to one summand as \eqn{L=1}. -Corresponding gAIC based weights for model $M$ are calculated as outlined in Schorning et al. (2016) +Corresponding gAIC based weights for model \eqn{M} are calculated as outlined in Schorning et al. (2016) \deqn{ \Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_(m))} } From f17ff92a94f38cb346fe1bc9a7e89e7dcaae1f1c Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Tue, 19 Dec 2023 14:51:14 +0000 Subject: [PATCH 32/39] Example for simulation added --- R/simulation.R | 16 ++++++++++++++++ man/simulateData.Rd | 16 ++++++++++++++++ 2 files changed, 32 insertions(+) diff --git a/R/simulation.R b/R/simulation.R index 09794e9..0828739 100644 --- a/R/simulation.R +++ b/R/simulation.R @@ -16,6 +16,22 @@ #' In case of NULL, all dose-response models, included in the mods input parameter will be used. #' @param dr_means a vector, with information about assumed effects per dose group. Default NULL. #' +#' @examples +#' # example code +#' models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +#' doses = c(0, 0.5, 2,4, 8),maxEff= 6) +#' dose_levels = c(0, 0.5, 2,4, 8) +#' sd = 12 +#' n_patients <- c(40,60,60,60,60) +#' sim_data <- simulateData( +#' n_patients = n_patients, +#' dose_levels = dose_levels, +#' mods = models, +#' n_sim = 100, +#' sd = sd) +#' +#' sim_data + #' @return sim_data one list object, containing patient level simulated data for all assumed true models. #' Also providing information about simulation iteration, patient number as well as dosage levels. #' diff --git a/man/simulateData.Rd b/man/simulateData.Rd index 1367b3a..8135e09 100644 --- a/man/simulateData.Rd +++ b/man/simulateData.Rd @@ -40,3 +40,19 @@ Also providing information about simulation iteration, patient number as well as \description{ Function to simulate patient level data for a normally distributed endpoint } +\examples{ +# example code + models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +doses = c(0, 0.5, 2,4, 8),maxEff= 6) +dose_levels = c(0, 0.5, 2,4, 8) +sd = 12 +n_patients <- c(40,60,60,60,60) +sim_data <- simulateData( +n_patients = n_patients, +dose_levels = dose_levels, +mods = models, +n_sim = 100, +sd = sd) + +sim_data +} From 46f6a7af682e17f8fea64a768ccd4a0206ecdf3b Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Wed, 20 Dec 2023 10:05:06 +0000 Subject: [PATCH 33/39] Minor correction in formulas of getModelfits documentation --- R/modelling.R | 4 ++-- man/getModelFits.Rd | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/R/modelling.R b/R/modelling.R index 3ef5daf..ed50b13 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -5,7 +5,7 @@ #' For the default case, the Nelder-Mead algorithm is used. #' In detail, for both approaches the mean vector \eqn{\theta^{Y}} and the covariance \eqn{\Sigma} of the (mixture) posterior distributions and the corresponding posterior weights \eqn{\tilde{\omega}_{l}} for \eqn{l \in {1,...,L}} are used as basis #' For the full fit a GLS estimator is used to minimize the following expression for the respective dose-response models \eqn{m} -#' \deqn{ \hat{\theta}_{m}=argmin_{\theta_{m}} \sum_{i=l}^{L} \tilde{\omega}_{l}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))'\Sigma_{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))} +#' \deqn{ \hat{\theta}_{m}=argmin_{\theta_{m}} \sum_{l=1}^{L} \tilde{\omega}_{l}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))'\Sigma_{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))} #' Herefore the function nloptr of the nloptr package is utilized. #' In the simplified case \eqn{L=1}, as the dimension of the posterior is reduced to 1 first. #' The generalized AIC values are calculated via the formula @@ -14,7 +14,7 @@ #' Here as well for the simplified case the formula reduces to one summand as \eqn{L=1}. #' Corresponding gAIC based weights for model \eqn{M} are calculated as outlined in Schorning et al. (2016) #' \deqn{ -#' \Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_(m))} +#' \Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_{m})} #' } #' where \eqn{Q} denotes the number of models included in the averaging procedure. #' @references Schorning K, Bornkamp B, Bretz F, Dette H. 2016. “Model selection versus model averaging in dose finding studies”. Stat Med; 35; 4021-4040. diff --git a/man/getModelFits.Rd b/man/getModelFits.Rd index fb49d23..f1b0953 100644 --- a/man/getModelFits.Rd +++ b/man/getModelFits.Rd @@ -24,7 +24,7 @@ For the simplified fit, multivariate normal distributions will be approximated a For the default case, the Nelder-Mead algorithm is used. In detail, for both approaches the mean vector \eqn{\theta^{Y}} and the covariance \eqn{\Sigma} of the (mixture) posterior distributions and the corresponding posterior weights \eqn{\tilde{\omega}_{l}} for \eqn{l \in {1,...,L}} are used as basis For the full fit a GLS estimator is used to minimize the following expression for the respective dose-response models \eqn{m} -\deqn{ \hat{\theta}_{m}=argmin_{\theta_{m}} \sum_{i=l}^{L} \tilde{\omega}_{l}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))'\Sigma_{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))} +\deqn{ \hat{\theta}_{m}=argmin_{\theta_{m}} \sum_{l=1}^{L} \tilde{\omega}_{l}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))'\Sigma_{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))} Herefore the function nloptr of the nloptr package is utilized. In the simplified case \eqn{L=1}, as the dimension of the posterior is reduced to 1 first. The generalized AIC values are calculated via the formula @@ -33,7 +33,7 @@ where \eqn{p} denotes the number of estimated parameters and \eqn{K} the number Here as well for the simplified case the formula reduces to one summand as \eqn{L=1}. Corresponding gAIC based weights for model \eqn{M} are calculated as outlined in Schorning et al. (2016) \deqn{ -\Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_(m))} +\Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_{m})} } where \eqn{Q} denotes the number of models included in the averaging procedure. } From 6772ec785c82b092e7f382a10adf2b2585073ab8 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Wed, 20 Dec 2023 10:27:53 +0000 Subject: [PATCH 34/39] Deletion of unnecessary code chunks in analysis example. --- vignettes/analysis_normal.Rmd | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 6da9c1a..870c640 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -205,25 +205,6 @@ contr_mat<- getContr( # dose_levels = dose_levels, # dose_weights = n_patients, # prior_list = prior_list) -#This would be the most reasonable output, but since it is not exported it is currently not usable. -#BMCP_result<-BayesMCPi(posterior_i = posterior, - # contr_mat = contr_mat_prior, - # crit_prob = crit_pval) - - -#This would be the most reasonable output, but since it is not exported it is currently not usable. -#BMCP_result<-BayesMCPi(posterior_i = posterior, - # contr_mat = contr_mat_prior, - # crit_prob = crit_pval) - - - -#BMCP_result2 <- performBayesianMCPMod( -# posteriors_list = posterior_emax, -# contr_mat = contr_mat_prior, -# crit_prob = crit_pval) - -#BMCP_result2 ``` The bayesian MCP testing step is then executed based on the posterior information, the provided contrasts and the (multiplicity adjusted) critical value. ```{r Executionr of Bayesian MCPMod Test step} From 3390fd7ba86aad605cc174afede09162098b79d6 Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Wed, 20 Dec 2023 10:47:08 +0000 Subject: [PATCH 35/39] Correction of typo (DoseFinding) and deletion of tests (so that package can be build). --- R/BMCPMod.R | 8 +- R/simulation.R | 2 +- man/assessDesign.Rd | 2 +- man/getContr.Rd | 4 +- man/getCritProb.Rd | 2 +- man/simulateData.Rd | 2 +- tests/testthat/test-bootstrapping.R | 152 ------------------------- tests/testthat/test-modelling.R | 27 ----- tests/testthat/test-plot.R | 171 ---------------------------- tests/testthat/test-posterior.R | 43 ------- 10 files changed, 10 insertions(+), 403 deletions(-) delete mode 100644 tests/testthat/test-bootstrapping.R delete mode 100644 tests/testthat/test-modelling.R delete mode 100644 tests/testthat/test-plot.R delete mode 100644 tests/testthat/test-posterior.R diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 37ea2a1..30fa6b0 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -3,7 +3,7 @@ #' @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 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 @@ -142,7 +142,7 @@ assessDesign <- function ( #' iv) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the #' regular MCPMod for this specific vector of standard errors. For the actual evaluation this vector of standard errors is translated into a (diagonal) matrix of variances #' -#' @param mods An object of class "Mods" as specified in the Dosefinding package. +#' @param mods An object of class "Mods" as specified in the DoseFinding package. #' @param dose_levels vector containing the different dosage levels. #' @param dose_weights Vector specifying weights for the different doses. Please note that in case this information is provided together with a prior (i.e. Option i) is planned) these two inputs should be provided on the same scale (e.g. patient numbers). Default NULL #' @param prior_list a prior_list object, only required as input for Option i). Default NULL @@ -160,7 +160,7 @@ assessDesign <- function ( #' dose_levels = dose_levels, #' sd_posterior = sd_posterior) #' -#' @return contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. +#' @return contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the DoseFinding package. #' #' @export getContr <- function ( @@ -241,7 +241,7 @@ getContr <- function ( #' ii) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the #' regular MCPMod for this specific vector of standard errors. Here as well the critical value for this set of contrasts is calculated via the critVal function of the DoseFinding package. #' -#' @param mods An object of class "Mods" as specified in the Dosefinding package. +#' @param mods An object of class "Mods" as specified in the DoseFinding package. #' @param dose_levels vector containing the different dosage levels. #' @param dose_weights Vector specifying weights for the different doses, only required for Option i). Default NULL #' @param se_new_trial a vector of positive values, only required for Option ii). Default NULL diff --git a/R/simulation.R b/R/simulation.R index 0828739..c4bd480 100644 --- a/R/simulation.R +++ b/R/simulation.R @@ -8,7 +8,7 @@ #' value per dose-group. #' @param dose_levels vector containing the different dosage levels. #' @param sd standard deviation on patient level. -#' @param mods An object of class "Mods" as specified in the Dosefinding package. +#' @param mods An object of class "Mods" as specified in the Doseinding package. #' @param n_sim number of simulations to be performed, #' Default is 1000 #' @param true_model Default value is NULL. diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd index 5b11892..7aaa902 100644 --- a/man/assessDesign.Rd +++ b/man/assessDesign.Rd @@ -21,7 +21,7 @@ assessDesign( \arguments{ \item{n_patients}{Vector specifying the planned number of patients per dose group} -\item{mods}{An object of class "Mods" as specified in the Dosefinding package.} +\item{mods}{An object of class "Mods" as specified in the DoseFinding package.} \item{prior_list}{a prior_list object specifying the utilized prior for the different dose groups} diff --git a/man/getContr.Rd b/man/getContr.Rd index 3abb8cc..8399bd2 100644 --- a/man/getContr.Rd +++ b/man/getContr.Rd @@ -14,7 +14,7 @@ getContr( ) } \arguments{ -\item{mods}{An object of class "Mods" as specified in the Dosefinding package.} +\item{mods}{An object of class "Mods" as specified in the DoseFinding package.} \item{dose_levels}{vector containing the different dosage levels.} @@ -27,7 +27,7 @@ getContr( \item{se_new_trial}{a vector of positive values with information about the observed variability, only required for Option iv). Default NULL} } \value{ -contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package. +contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the DoseFinding package. } \description{ This function calculates contrast vectors that are optimal for detecting certain alternatives via applying the function optContr of the DoseFinding package. diff --git a/man/getCritProb.Rd b/man/getCritProb.Rd index e6ff1bb..2e4f400 100644 --- a/man/getCritProb.Rd +++ b/man/getCritProb.Rd @@ -13,7 +13,7 @@ getCritProb( ) } \arguments{ -\item{mods}{An object of class "Mods" as specified in the Dosefinding package.} +\item{mods}{An object of class "Mods" as specified in the DoseFinding package.} \item{dose_levels}{vector containing the different dosage levels.} diff --git a/man/simulateData.Rd b/man/simulateData.Rd index 8135e09..86db979 100644 --- a/man/simulateData.Rd +++ b/man/simulateData.Rd @@ -22,7 +22,7 @@ value per dose-group.} \item{sd}{standard deviation on patient level.} -\item{mods}{An object of class "Mods" as specified in the Dosefinding package.} +\item{mods}{An object of class "Mods" as specified in the Doseinding package.} \item{n_sim}{number of simulations to be performed, Default is 1000} diff --git a/tests/testthat/test-bootstrapping.R b/tests/testthat/test-bootstrapping.R deleted file mode 100644 index c4caf87..0000000 --- a/tests/testthat/test-bootstrapping.R +++ /dev/null @@ -1,152 +0,0 @@ -# Test cases -test_that("test getBootstrapBands", { - library(clinDR) - library(dplyr) - - data("metaData") - testdata <- as.data.frame(metaData) - dataset <- filter(testdata, bname == "VIAGRA") - histcontrol <- filter(dataset, dose == 0, primtime == 12, indication == "ERECTILE DYSFUNCTION") - - ##Create MAP Prior - hist_data <- data.frame( - trial = c(1, 2, 3, 4), - est = histcontrol$rslt, - se = histcontrol$se, - sd = histcontrol$sd, - n = histcontrol$sampsize) - - 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::robustify( - priormix = RBesT::automixfit(gmap), - weight = 0.5, - mean = 1.4, - sigma = sd_tot) - - #RBesT::ess(prior_ctr) - - ## derive prior for treatment - ## weak informative using same parameters as for robustify component - prior_trt <- RBesT::mixnorm( - comp1 = c(w = 1, m = 1.4, n = 1), - sigma = sd_tot, - param = "mn") - dose_levels <- c(0, 50, 100, 200) - ## combine priors in list - prior_list <- c(list(prior_ctr), rep(list(prior_trt), times = length(dose_levels[-1]))) - - #Pre-Specification (B)MCPMod - - ## candidate models for MCPMod - # linear function - no guestimates needed - exp <- DoseFinding::guesst(d = 50, - p = c(0.2), - model = "exponential", - Maxd = max(dose_levels)) - emax <- DoseFinding::guesst(d = 100, - p = c(0.9), - model = "emax") - - - mods <- DoseFinding::Mods( - linear = NULL, - emax = emax, - exponential = exp, - doses = dose_levels, - maxEff = 10, - placEff = 1.4) - - #Simulation of new trial - ##Note: This part will be simplified and direct results from one trial will be used - mods_sim <- DoseFinding::Mods( - emax = emax, - doses = dose_levels, - maxEff = 12, - placEff = 1.4) - - n_patients <- c(50, 50, 50, 50) - data <- simulateData( - n_patients = n_patients, - dose_levels = dose_levels, - sd = sd_tot, - mods = mods_sim, - n_sim = 1) - - data_emax <- data[, c("simulation", "dose", "emax")] - names(data_emax)[3] <- "response" - - posterior_emax <- getPosterior( - data = data_emax, - prior_list = prior_list) - - #Evaluation of Bayesian MCPMod - - contr_mat <- DoseFinding::optContr( - models = mods, - doses = dose_levels, - w = n_patients) - ##Calculation of critical value can be done with critVal function - crit_val_equal <- DoseFinding:::critVal(contr_mat$corMat, alpha = 0.05, df = 0, alternative = "one.sided") - crit_pval <- pnorm(crit_val_equal) - - ess_prior <- round(unlist(lapply(prior_list, RBesT::ess))) - - ### Evaluation of Bayesian MCPMod - contr_mat_prior <- DoseFinding::optContr( - models = mods, - doses = dose_levels, - w = n_patients + ess_prior) - - BMCP_result <- performBayesianMCP( - posteriors_list = list(posterior_emax), - contr_mat = contr_mat_prior, - crit_prob = crit_pval) - - BMCP_result - - #Model fit - #This part is currently not working - post_observed <- posterior_emax - model_shapes <- c("linear", "emax", "exponential") - - # Option a) Simplified approach by using frequentist idea - fit_simple <- getModelFits( - models = model_shapes, - dose_levels = dose_levels, -posterior = post_observed, - simple = TRUE) - - # Option b) Making use of the complete posterior distribution - fit <- getModelFits( - models = model_shapes, - dose_levels = dose_levels, - posterior = post_observed, - simple = FALSE) - - result_simple <- getBootstrapBands(fit_simple) - result <- getBootstrapBands(fit) - expect_type(result_simple, "list") - expect_type(result, "list") - - result_2_simple <- getBootstrapBands(fit_simple, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3)) - result_2 <- getBootstrapBands(fit, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3)) - expect_type(result_2_simple, "list") - expect_type(result_2, "list") - - result_3_simple <- getBootstrapBands(fit_simple, dose_seq = NULL) - result_3 <- getBootstrapBands(fit, dose_seq = NULL) - expect_type(result_3_simple, "list") - expect_type(result_3, "list") -}) - diff --git a/tests/testthat/test-modelling.R b/tests/testthat/test-modelling.R deleted file mode 100644 index e398d28..0000000 --- a/tests/testthat/test-modelling.R +++ /dev/null @@ -1,27 +0,0 @@ -# Test data -test_data <- data.frame( - simulation = rep(1, 6), - dose = c(0, 1, 2, 3, 4, 5), - response = c(0, 1, 2, 3, 4, 5) -) - -# Mock getPosterior function -getPosterior <- function(data, prior_list, mu_hat, sd_hat) { - list( - means = c(0, 1, 2, 3, 4, 5), - vars = c(1, 1, 1, 1, 1, 1), - weights = c(1, 1, 1, 1, 1, 1) - ) -} - -# Test predictModelFit function -test_that("predictModelFit works correctly", { - model_fit <- list( - model = "emax", - coeffs = c(e0 = 0, eMax = 1, ed50 = 2), - dose_levels = c(0, 1, 2, 3, 4, 5) - ) - - pred_vals <- predictModelFit(model_fit) - expect_type(pred_vals, "double") -}) diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R deleted file mode 100644 index 732580d..0000000 --- a/tests/testthat/test-plot.R +++ /dev/null @@ -1,171 +0,0 @@ -test_that("Test plot.modelFits with different model_fits input", { - library(BayesianMCPMod) - library(clinDR) - library(dplyr) - data("metaData") - testdata <- as.data.frame(metaData) - dataset <- filter(testdata, bname == "VIAGRA") - histcontrol <- filter(dataset, dose == 0, primtime == 12, indication == "ERECTILE DYSFUNCTION") - - ## Create MAP Prior - hist_data <- data.frame( - trial = c(1, 2, 3, 4), - est = histcontrol$rslt, - se = histcontrol$se, - sd = histcontrol$sd, - n = histcontrol$sampsize - ) - - 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::robustify( - priormix = RBesT::automixfit(gmap), - weight = 0.5, - mean = 1.4, - sigma = sd_tot - ) - - # RBesT::ess(prior_ctr) - - ## derive prior for treatment - ## weak informative using same parameters as for robustify component - prior_trt <- RBesT::mixnorm( - comp1 = c(w = 1, m = 1.4, n = 1), - sigma = sd_tot, - param = "mn" - ) - dose_levels <- c(0, 50, 100, 200) - ## combine priors in list - prior_list <- c(list(prior_ctr), rep(list(prior_trt), times = length(dose_levels[-1]))) - - # Pre-Specification (B)MCPMod - - ## candidate models for MCPMod - # linear function - no guestimates needed - exp <- DoseFinding::guesst( - d = 50, - p = c(0.2), - model = "exponential", - Maxd = max(dose_levels) - ) - emax <- DoseFinding::guesst( - d = 100, - p = c(0.9), - model = "emax" - ) - - - mods <- DoseFinding::Mods( - linear = NULL, - emax = emax, - exponential = exp, - doses = dose_levels, - maxEff = 10, - placEff = 1.4 - ) - - # Simulation of new trial - ## Note: This part will be simplified and direct results from one trial will be used - mods_sim <- DoseFinding::Mods( - emax = emax, - doses = dose_levels, - maxEff = 12, - placEff = 1.4 - ) - - n_patients <- c(50, 50, 50, 50) - data <- simulateData( - n_patients = n_patients, - dose_levels = dose_levels, - sd = sd_tot, - mods = mods_sim, - n_sim = 1 - ) - - data_emax <- data[, c("simulation", "dose", "emax")] - names(data_emax)[3] <- "response" - - posterior_emax <- getPosterior( - data = data_emax, - prior_list = prior_list - ) - - # Evaluation of Bayesian MCPMod - - contr_mat <- DoseFinding::optContr( - models = mods, - doses = dose_levels, - w = n_patients - ) - ## Calculation of critical value can be done with critVal function - crit_val_equal <- DoseFinding:::critVal(contr_mat$corMat, alpha = 0.05, df = 0, alternative = "one.sided") - crit_pval <- stats::pnorm(crit_val_equal) - - ess_prior <- round(unlist(lapply(prior_list, RBesT::ess))) - - ### Evaluation of Bayesian MCPMod - contr_mat_prior <- DoseFinding::optContr( - models = mods, - doses = dose_levels, - w = n_patients + ess_prior - ) - - BMCP_result <- performBayesianMCP( - posteriors_list = list(posterior_emax), - contr_mat = contr_mat_prior, - crit_prob = crit_pval - ) - - # Model fit - # This part is currently not working - post_observed <- posterior_emax - model_shapes <- c("linear", "emax", "exponential") - - - # Option a) Simplified approach by using frequentist idea - fit_simple <- getModelFits( - models = model_shapes, - dose_levels = dose_levels, - posterior = post_observed, - simple = TRUE - ) - - # Option b) Making use of the complete posterior distribution - fit <- getModelFits( - models = model_shapes, - dose_levels = dose_levels, - posterior = post_observed, - simple = FALSE - ) - - # Test with default parameters and more models - plot1 <- plot.modelFits(fit) - expect_s3_class(plot1, "ggplot") - - # Test with cr_intv = TRUE and more models - plot2 <- plot.modelFits(fit, cr_intv = TRUE) - expect_s3_class(plot2, "ggplot") - - # Test with gAIC = FALSE and more models - plot3 <- plot.modelFits(fit, gAIC = FALSE) - expect_s3_class(plot3, "ggplot") - - # Test with avg_fit = FALSE and more models - plot4 <- plot.modelFits(fit, avg_fit = FALSE) - expect_s3_class(plot4, "ggplot") - - # Test with all non-default parameters and more models - plot5 <- plot.modelFits(fit, cr_intv = TRUE, gAIC = FALSE, avg_fit = FALSE) - expect_s3_class(plot5, "ggplot") -}) \ No newline at end of file diff --git a/tests/testthat/test-posterior.R b/tests/testthat/test-posterior.R deleted file mode 100644 index 47eeb69..0000000 --- a/tests/testthat/test-posterior.R +++ /dev/null @@ -1,43 +0,0 @@ -test_that("getPosterior works correctly", { - # Prepare test data and parameters - data <- data.frame(simulation = rep(1, 4), - dose = c(0, 1, 2, 3), - response = c(10, 20, 30, 40)) - prior_list <- list(1, 2, 3, 4) - mu_hat <- c(10, 20, 30, 40) - sd_hat <- matrix(c(1, 2, 3, 4), nrow = 4, ncol = 1) - - # Test getPosterior function - posterior_list <- getPosterior(data, prior_list, mu_hat, sd_hat) - expect_type(posterior_list, "character") - expect_s3_class(posterior_list, "postList") -}) - -test_that("getPosteriorI works correctly", { - # Prepare test data and parameters - data_i <- data.frame(dose = c(0, 1, 2, 3), - response = c(10, 20, 30, 40)) - prior_list <- list(1, 2, 3, 4) - mu_hat <- c(10, 20, 30, 40) - sd_hat <- matrix(c(1, 2, 3, 4), nrow = 4, ncol = 1) - - # Test getPosteriorI function - post_list <- getPosteriorI(data_i, prior_list, mu_hat, sd_hat) - expect_type(post_list, "character") - expect_s3_class(post_list, "postList") -}) - -test_that("summary.postList works correctly", { - # Prepare test data - post_list <- list( - Ctr = matrix(c(0.25, 10, 1), nrow = 3, ncol = 1), - DG_1 = matrix(c(0.25, 20, 2), nrow = 3, ncol = 1), - DG_2 = matrix(c(0.25, 30, 3), nrow = 3, ncol = 1), - DG_3 = matrix(c(0.25, 40, 4), nrow = 3, ncol = 1) - ) - class(post_list) <- "postList" - - # Test summary.postList function - summary_tab <- summary.postList(post_list) - expect_type(summary_tab, "character") -}) From 5c4cfaf5ffeb7faa3749ec23758db35599ab81da Mon Sep 17 00:00:00 2001 From: sebastianbossert Date: Wed, 20 Dec 2023 11:02:35 +0000 Subject: [PATCH 36/39] Update of error in example code in assessdesign --- R/BMCPMod.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 30fa6b0..a43b9ba 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -17,7 +17,7 @@ #' #' @examples #' # example code -#' models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, +#' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, #' doses = c(0, 0.5, 2,4, 8),maxEff= 6) #' sd = 12 #' prior_list<-list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2), From 9b9857d76675bfba114e783a92d35408619f342a Mon Sep 17 00:00:00 2001 From: Andersen Date: Wed, 20 Dec 2023 13:05:56 +0100 Subject: [PATCH 37/39] minor fixes in examples --- R/BMCPMod.R | 15 +++++---------- man/assessDesign.Rd | 3 +-- man/getContr.Rd | 3 +-- man/getCritProb.Rd | 3 +-- man/performBayesianMCP.Rd | 3 +-- man/performBayesianMCPMod.Rd | 3 +-- 6 files changed, 10 insertions(+), 20 deletions(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index a43b9ba..10e1b30 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -17,8 +17,7 @@ #' #' @examples #' # example code -#' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, -#' doses = c(0, 0.5, 2,4, 8),maxEff= 6) +#' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8), maxEff= 6) #' sd = 12 #' prior_list<-list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2), #' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), @@ -151,8 +150,7 @@ assessDesign <- function ( #' #' @examples #' # example code -#' models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, -#' doses = c(0, 0.5, 2,4, 8)) +#' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8), maxEff= 6) #' dose_levels=c(0, 0.5, 2, 4, 8) #' sd_posterior = c(2.8,3,2.5,3.5,4) #' contr_mat<- getContr( @@ -249,8 +247,7 @@ getContr <- function ( #' #' @examples #' # example code -#' models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, -#' doses = c(0, 0.5, 2,4, 8)) +#' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8)) #' dose_levels=c(0, 0.5, 2, 4, 8) #' critVal<- getCritProb( #' mods = models, @@ -295,8 +292,7 @@ getCritProb <- function ( #' @param simple boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE. #' @examples #' # example code -#' models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, -#' doses = c(0, 0.5, 2,4, 8)) +#' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8)) #' dose_levels=c(0, 0.5, 2, 4, 8) #' sd_posterior = c(2.8,3,2.5,3.5,4) #' contr_mat<- getContr( @@ -422,8 +418,7 @@ addSignificance <- function ( #' #' @examples #' # example code -#' models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, -#' doses = c(0, 0.5, 2,4, 8)) +#' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8)) #' dose_levels=c(0, 0.5, 2, 4, 8) #' sd_posterior = c(2.8,3,2.5,3.5,4) #' contr_mat<- getContr( diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd index 7aaa902..1f0a539 100644 --- a/man/assessDesign.Rd +++ b/man/assessDesign.Rd @@ -47,8 +47,7 @@ This function performs simulation based trial design evaluations for a set of sp } \examples{ # example code - models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, -doses = c(0, 0.5, 2,4, 8),maxEff= 6) +mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8), maxEff= 6) sd = 12 prior_list<-list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2), DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), diff --git a/man/getContr.Rd b/man/getContr.Rd index 8399bd2..f3813a4 100644 --- a/man/getContr.Rd +++ b/man/getContr.Rd @@ -45,8 +45,7 @@ regular MCPMod for this specific vector of standard errors. For the actual evalu } \examples{ # example code -models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, -doses = c(0, 0.5, 2,4, 8)) +mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8), maxEff= 6) dose_levels=c(0, 0.5, 2, 4, 8) sd_posterior = c(2.8,3,2.5,3.5,4) contr_mat<- getContr( diff --git a/man/getCritProb.Rd b/man/getCritProb.Rd index 2e4f400..bc17bb9 100644 --- a/man/getCritProb.Rd +++ b/man/getCritProb.Rd @@ -37,8 +37,7 @@ regular MCPMod for this specific vector of standard errors. Here as well the cri } \examples{ # example code -models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, -doses = c(0, 0.5, 2,4, 8)) +mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8)) dose_levels=c(0, 0.5, 2, 4, 8) critVal<- getCritProb( mods = models, diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd index 3fec02c..0e3140e 100644 --- a/man/performBayesianMCP.Rd +++ b/man/performBayesianMCP.Rd @@ -23,8 +23,7 @@ In order to obtain significant test decision we consider the maximum of the post } \examples{ # example code -models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, -doses = c(0, 0.5, 2,4, 8)) +mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8)) dose_levels=c(0, 0.5, 2, 4, 8) sd_posterior = c(2.8,3,2.5,3.5,4) contr_mat<- getContr( diff --git a/man/performBayesianMCPMod.Rd b/man/performBayesianMCPMod.Rd index e02ffd1..24a08a9 100644 --- a/man/performBayesianMCPMod.Rd +++ b/man/performBayesianMCPMod.Rd @@ -23,8 +23,7 @@ performs bayesian MCP Test step and modelling in a combined fashion. See perform } \examples{ # example code -models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, -doses = c(0, 0.5, 2,4, 8)) +mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8)) dose_levels=c(0, 0.5, 2, 4, 8) sd_posterior = c(2.8,3,2.5,3.5,4) contr_mat<- getContr( From 2de0df6b80244a277cc2cc541685715896b57d9e Mon Sep 17 00:00:00 2001 From: Andersen Date: Wed, 20 Dec 2023 13:51:17 +0100 Subject: [PATCH 38/39] minor fix --- R/BMCPMod.R | 32 ++++++++++++++++---------------- man/getContr.Rd | 2 +- man/getCritProb.Rd | 2 +- man/performBayesianMCP.Rd | 4 ++-- man/performBayesianMCPMod.Rd | 4 ++-- 5 files changed, 22 insertions(+), 22 deletions(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 10e1b30..4dafc4e 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -92,9 +92,9 @@ assessDesign <- function ( post_sds <- sapply(posterior_list, function (post) summary(post)[, 2]) contr <- apply(post_sds, 2, function (post_sd) getContr( - mods = mods, - dose_levels = dose_levels, - sd_posterior = post_sd)) + mods = mods, + dose_levels = dose_levels, + sd_posterior = post_sd)) } @@ -154,7 +154,7 @@ assessDesign <- function ( #' dose_levels=c(0, 0.5, 2, 4, 8) #' sd_posterior = c(2.8,3,2.5,3.5,4) #' contr_mat<- getContr( -#' mods = models, +#' mods = mods, #' dose_levels = dose_levels, #' sd_posterior = sd_posterior) #' @@ -179,21 +179,21 @@ getContr <- function ( w <- NULL S <- diag((se_new_trial)^2) - # frequentist & no re-estimation + # frequentist & no re-estimation } else if (!is.null(dose_weights) & is.null(se_new_trial) & is.null(prior_list) & is.null(sd_posterior)) { w <- dose_weights S <- NULL - # Bayesian & re-estimation + # Bayesian & re-estimation } else if (!is.null(sd_posterior) & is.null(se_new_trial) & is.null(prior_list) & is.null(dose_weights)) { w <- NULL S <- diag((sd_posterior)^2) - # Bayesian & no re-estimation + # Bayesian & no re-estimation } else if (!is.null(dose_weights) & !is.null(prior_list) & is.null(se_new_trial) & is.null(sd_posterior)) { @@ -204,8 +204,8 @@ getContr <- function ( } else { stop (paste("Provided combiations of 'se_new_trial',", - "'dose_weights', 'prior_list', 'sd_posterior' not allowed.", - "See ?getContr for allowed combinations.")) + "'dose_weights', 'prior_list', 'sd_posterior' not allowed.", + "See ?getContr for allowed combinations.")) } @@ -250,7 +250,7 @@ getContr <- function ( #' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8)) #' dose_levels=c(0, 0.5, 2, 4, 8) #' critVal<- getCritProb( -#' mods = models, +#' mods = mods, #' dose_weights =c(50,50,50,50,50), #reflecting the planned sample size #' dose_levels = dose_levels, #' alpha_crit_val = 0.05) @@ -296,11 +296,11 @@ getCritProb <- function ( #' dose_levels=c(0, 0.5, 2, 4, 8) #' sd_posterior = c(2.8,3,2.5,3.5,4) #' contr_mat<- getContr( -#' mods = models, +#' mods = mods, #' dose_levels = dose_levels, #' sd_posterior = sd_posterior) #' critVal<- getCritProb( -#' mods = models, +#' mods = mods, #' dose_weights =c(50,50,50,50,50), #reflecting the planned sample size #' dose_levels = dose_levels, #' alpha_crit_val = 0.05) @@ -389,7 +389,7 @@ addSignificance <- function ( model_fits, sign_models - + ) { names(sign_models) <- NULL @@ -422,11 +422,11 @@ addSignificance <- function ( #' dose_levels=c(0, 0.5, 2, 4, 8) #' sd_posterior = c(2.8,3,2.5,3.5,4) #' contr_mat<- getContr( -#' mods = models, +#' mods = mods, #' dose_levels = dose_levels, #' sd_posterior = sd_posterior) #' critVal<- getCritProb( -#' mods = models, +#' mods = mods, #' dose_weights =c(50,50,50,50,50), #reflecting the planned sample size #' dose_levels = dose_levels, #' alpha_crit_val = 0.05) @@ -527,4 +527,4 @@ BayesMCPi <- function ( return (res) -} +} \ No newline at end of file diff --git a/man/getContr.Rd b/man/getContr.Rd index f3813a4..c1e5f20 100644 --- a/man/getContr.Rd +++ b/man/getContr.Rd @@ -49,7 +49,7 @@ mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), expo dose_levels=c(0, 0.5, 2, 4, 8) sd_posterior = c(2.8,3,2.5,3.5,4) contr_mat<- getContr( -mods = models, +mods = mods, dose_levels = dose_levels, sd_posterior = sd_posterior) diff --git a/man/getCritProb.Rd b/man/getCritProb.Rd index bc17bb9..39593f6 100644 --- a/man/getCritProb.Rd +++ b/man/getCritProb.Rd @@ -40,7 +40,7 @@ regular MCPMod for this specific vector of standard errors. Here as well the cri mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8)) dose_levels=c(0, 0.5, 2, 4, 8) critVal<- getCritProb( - mods = models, + mods = mods, dose_weights =c(50,50,50,50,50), #reflecting the planned sample size dose_levels = dose_levels, alpha_crit_val = 0.05) diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd index 0e3140e..cc0eaa6 100644 --- a/man/performBayesianMCP.Rd +++ b/man/performBayesianMCP.Rd @@ -27,11 +27,11 @@ mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), expo dose_levels=c(0, 0.5, 2, 4, 8) sd_posterior = c(2.8,3,2.5,3.5,4) contr_mat<- getContr( -mods = models, +mods = mods, dose_levels = dose_levels, sd_posterior = sd_posterior) critVal<- getCritProb( - mods = models, + mods = mods, dose_weights =c(50,50,50,50,50), #reflecting the planned sample size dose_levels = dose_levels, alpha_crit_val = 0.05) diff --git a/man/performBayesianMCPMod.Rd b/man/performBayesianMCPMod.Rd index 24a08a9..b8745ef 100644 --- a/man/performBayesianMCPMod.Rd +++ b/man/performBayesianMCPMod.Rd @@ -27,11 +27,11 @@ mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), expo dose_levels=c(0, 0.5, 2, 4, 8) sd_posterior = c(2.8,3,2.5,3.5,4) contr_mat<- getContr( -mods = models, +mods = mods, dose_levels = dose_levels, sd_posterior = sd_posterior) critVal<- getCritProb( - mods = models, + mods = mods, dose_weights =c(50,50,50,50,50), #reflecting the planned sample size dose_levels = dose_levels, alpha_crit_val = 0.05) From b241cf678ae8365ae68229ac2b5ffe3fe38eb9f7 Mon Sep 17 00:00:00 2001 From: Andersen Date: Wed, 20 Dec 2023 14:37:31 +0100 Subject: [PATCH 39/39] added optional param for predict.modelFits func --- R/s3methods.R | 1 + man/predict.modelFits.Rd | 2 ++ 2 files changed, 3 insertions(+) diff --git a/R/s3methods.R b/R/s3methods.R index 6fd93ce..a507df9 100644 --- a/R/s3methods.R +++ b/R/s3methods.R @@ -85,6 +85,7 @@ print.BayesianMCP <- function ( #' model coefficients #' @param doses a vector specifying the doses for which a prediction should be #' done getContrMat object, contrast matrix to be used for the testing step. +#' @param ... other parameters #' @examples #' # example code #' posterior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2), diff --git a/man/predict.modelFits.Rd b/man/predict.modelFits.Rd index 7fecd10..7563063 100644 --- a/man/predict.modelFits.Rd +++ b/man/predict.modelFits.Rd @@ -12,6 +12,8 @@ model coefficients} \item{doses}{a vector specifying the doses for which a prediction should be done getContrMat object, contrast matrix to be used for the testing step.} + +\item{...}{other parameters} } \value{ a list with the model predictions for the specified models and doses