diff --git a/.Rbuildignore b/.Rbuildignore index 9eaae49..bda7f19 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -11,3 +11,4 @@ docs ^cran-comments\.md$ ^cran_submission_script\.R$ ^CRAN-SUBMISSION$ +^vignettes/*_files$ diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml index 8163e3c..8c2bc42 100644 --- a/.github/workflows/linux.yml +++ b/.github/workflows/linux.yml @@ -24,6 +24,7 @@ jobs: R_KEEP_PKG_SOURCE: yes steps: + - uses: quarto-dev/quarto-actions/setup@v2 - uses: actions/checkout@v2 - uses: r-lib/actions/setup-pandoc@v2 diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml index 2b4c2c9..a4a3309 100644 --- a/.github/workflows/macos.yml +++ b/.github/workflows/macos.yml @@ -26,6 +26,7 @@ jobs: R_KEEP_PKG_SOURCE: yes steps: + - uses: quarto-dev/quarto-actions/setup@v2 - uses: actions/checkout@v2 - uses: r-lib/actions/setup-pandoc@v2 diff --git a/.github/workflows/pkgdown.yml b/.github/workflows/pkgdown.yml index e7d01f1..1fa6370 100644 --- a/.github/workflows/pkgdown.yml +++ b/.github/workflows/pkgdown.yml @@ -20,6 +20,7 @@ jobs: env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: + - uses: quarto-dev/quarto-actions/setup@v2 - uses: actions/checkout@v2 - uses: r-lib/actions/setup-pandoc@v2 @@ -27,7 +28,7 @@ jobs: - uses: r-lib/actions/setup-r@v2 with: use-public-rspm: true - + - name: Install JAGS run: | sudo apt-get update -y diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index 68ff44d..fe07329 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -25,6 +25,7 @@ jobs: R_KEEP_PKG_SOURCE: yes steps: + - uses: quarto-dev/quarto-actions/setup@v2 - uses: actions/checkout@v2 - uses: r-lib/actions/setup-pandoc@v2 @@ -40,6 +41,10 @@ jobs: extra-packages: any::rcmdcheck needs: check + - name: Install BayesianMCPMod + shell: bash + run: R CMD INSTALL --preclean . + - uses: r-lib/actions/check-r-package@v2 with: upload-snapshots: true diff --git a/DESCRIPTION b/DESCRIPTION index 8a05aeb..c9ddc6b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -6,8 +6,8 @@ Authors@R: c( person("Boehringer Ingelheim Pharma GmbH & Co. KG", role = c("cph", "fnd")), person("Stephan", "Wojciekowski", , "stephan.wojciekowski@boehringer-ingelheim.com", role = c("aut", "cre")), person("Lars", "Andersen", , "lars.andersen@boehringer-ingelheim.com", role = "aut"), - person("Steven", "Brooks", , "steven.brooks@boehringer-ingelheim.com", role = "ctb"), - person("Sebastian", "Bossert", , "sebastian.bossert@boehringer-ingelheim.com", role = "aut") + person("Sebastian", "Bossert", , "sebastian.bossert@boehringer-ingelheim.com", role = "aut"), + person("Steven", "Brooks", , "steven.brooks@boehringer-ingelheim.com", role = "ctb") ) Description: Bayesian MCPMod (Fleischer et al. (2022) ) is an innovative method that improves the @@ -38,17 +38,19 @@ Imports: RBesT, stats Suggests: + reactable, + tibble, + quarto, clinDR, dplyr, knitr, rmarkdown, spelling, testthat (>= 3.0.0) -VignetteBuilder: - knitr +VignetteBuilder: quarto Config/testthat/edition: 3 Encoding: UTF-8 Language: en-US LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.1 diff --git a/NAMESPACE b/NAMESPACE index 499ea56..dc6a5e4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,7 +9,6 @@ S3method(print,postList) S3method(summary,postList) export(assessDesign) export(getBootstrapQuantiles) -export(getBootstrapSamples) export(getContr) export(getCritProb) export(getESS) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 10c5ecc..41f00e1 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -361,7 +361,7 @@ getCritProb <- function ( #' crit_prob_adj = critVal, #' simple = FALSE) #' -#' @return Bayesian MCP test result as well as modelling result. +#' @return Bayesian MCP test result as well as modeling result. #' #' @export performBayesianMCPMod <- function ( @@ -399,7 +399,7 @@ performBayesianMCPMod <- function ( stop ("Argument 'contr' must be of type 'optContr'") } - + b_mcp <- performBayesianMCP( posterior_list = posterior_list, contr = contr, @@ -408,8 +408,8 @@ performBayesianMCPMod <- function ( fits_list <- lapply(seq_along(posterior_list), function (i) { if (b_mcp[i, 1]) { - - sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "critProbAdj") + + sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "crit_prob_adj") model_fits <- getModelFits( models = model_shapes, diff --git a/R/bootstrapping.R b/R/bootstrapping.R index df1734f..4a28b7e 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -1,12 +1,14 @@ -#' @title getBootstrapSamples +#' @title getBootstrapQuantiles #' -#' @description A function for the calculation of bootstrapped model predictions. -#' 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. +#' @description A function for the calculation of bootstrapped model predictions. +#' 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). #' 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. 2017. Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials (1st ed.). Chapman and Hall/CRC. doi: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 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 bootstrapped quantiles #' @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. @@ -14,7 +16,7 @@ #' @examples #' 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_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("exponential", "linear") @@ -23,40 +25,43 @@ #' posterior = posterior_list, #' dose_levels = dose_levels, #' simple = TRUE) -#' -#' getBootstrapSamples(model_fits = fit, +#' +#' getBootstrapQuantiles(model_fits = fit, +#' quantiles = c(0.025, 0.5, 0.975), #' n_samples = 10, # speeding up example run time #' doses = c(0, 6, 8)) -#' -#' @return A data frame with columns for model, dose, and bootstrapped samples -#' +#' +#' @return A data frame with columns for model, dose, and bootstrapped samples +#' #' @export -getBootstrapSamples <- function ( - +getBootstrapQuantiles <- function ( + model_fits, + quantiles, n_samples = 1e3, doses = NULL, avg_fit = TRUE - + ) { - + mu_hat_samples <- sapply(attr(model_fits, "posterior"), RBesT::rmix, n = n_samples) sd_hat <- summary.postList(attr(model_fits, "posterior"))[, 2] - + dose_levels <- model_fits[[1]]$dose_levels model_names <- names(model_fits) - + quantile_probs <- sort(unique(quantiles)) + if (is.null(doses)) { - + doses <- seq(min(dose_levels), max(dose_levels), length.out = 100L) - + } - + preds <- apply(mu_hat_samples, 1, function (mu_hat) { - + preds_mu_hat <- sapply(model_names, function (model) { - + fit <- DoseFinding::fitMod( dose = model_fits[[1]]$dose_levels, resp = mu_hat, @@ -65,101 +70,50 @@ getBootstrapSamples <- function ( type = "general", bnds = DoseFinding::defBnds( mD = max(model_fits[[1]]$dose_levels))[[model]]) - + preds <- stats::predict(fit, doseSeq = doses, predType = "ls-means") attr(preds, "gAIC") <- DoseFinding::gAIC(fit) - + return (preds) - + }, simplify = FALSE) - + preds_mu_mat <- do.call(rbind, preds_mu_hat) - + if (avg_fit) { - + avg_fit_indx <- which.min(sapply(preds_mu_hat, attr, "gAIC")) preds_mu_mat <- rbind(preds_mu_mat, avgFit = preds_mu_mat[avg_fit_indx, ]) - + } - + return (preds_mu_mat) - + }) - + if (avg_fit) { - + model_names <- c(model_names, "avgFit") - + } - - bs_samples <- as.data.frame(preds[order(rep(seq_along(model_names), length(doses))), ]) - colnames(bs_samples) <- paste0("preds_", seq_len(n_samples)) - - bs_samples_data <- cbind( - model = rep( - x = factor(model_names, - levels = c("linear", "emax", "exponential", - "sigEmax", "logistic", "quadratic", - "avgFit")), - each = length(doses)), - dose = doses, - bs_samples) - - # tidyr::pivot_longer(bs_samples_data, cols = contains("preds")) - - return (bs_samples_data) - -} -#' @title getBootstrapQuantiles -#' -#' @description Calculates quantiles from bootstrapped dose predictions. -#' Can be used to derive credible intervals to assess the uncertainty for the model fit. -#' @param bs_samples An object of class bootstrappedSample as created by getBootstrapSamples -#' @param quantiles A vector of quantiles that should be evaluated -#' -#' @examples -#' 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("exponential", "linear") -#' dose_levels <- c(0, 1, 2, 4, 8) -#' fit <- getModelFits(models = models, -#' posterior = posterior_list, -#' dose_levels = dose_levels, -#' simple = TRUE) -#' -#' bs_samples <- getBootstrapSamples(model_fits = fit, -#' n_samples = 10, # speeding up example run time -#' doses = c(0, 6, 8)) -#' -#' getBootstrapQuantiles(bs_samples = bs_samples, -#' quantiles = c(0.025, 0.5, 0.975)) -#' @return A data frame with entries doses, models, and quantiles -#' -#' @export -getBootstrapQuantiles <- function ( + 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)) - bs_samples, - quantiles -) { - - quantile_probs <- sort(unique(quantiles)) - - bs_quantiles <- t(apply(X = bs_samples[, -c(1, 2)], - MARGIN = 1, - FUN = stats::quantile, - probs = quantile_probs)) - - bs_quantiles_data <- cbind( - bs_samples[, c(1, 2)], - as.data.frame(bs_quantiles)) - - return (bs_quantiles_data) + cr_bounds_data <- cbind( + doses = doses, + models = rep( + x = factor(model_names, levels = model_names), + each = length(doses)), + as.data.frame(quant_mat)) + + return (cr_bounds_data) } + diff --git a/R/modelling.R b/R/modeling.R similarity index 69% rename from R/modelling.R rename to R/modeling.R index d52da7f..1d52eca 100644 --- a/R/modelling.R +++ b/R/modeling.R @@ -1,12 +1,12 @@ #' @title getModelFits -#' +#' #' @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. +#' 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 \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_{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}))} -#' Therefore the function nloptr of the nloptr package is utilized. +#' Therefore 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 #' \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 } @@ -20,17 +20,17 @@ #' @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. +#' @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 #' 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_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) @@ -38,49 +38,52 @@ #' posterior = posterior_list, #' dose_levels = dose_levels, #' simple = TRUE) -#' +#' #' @return An object of class modelFits. 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 ( - + models, dose_levels, posterior, simple = FALSE - + ) { - + if (inherits(models, "character")) { + models <- stats::setNames(as.list(models), models) + } checkmate::check_list(models, any.missing = FALSE) checkmate::check_double(dose_levels, lower = 0, any.missing = FALSE, len = length(models)) checkmate::check_class(posterior, "postList") checkmate::check_logical(simple) - - models <- unique(gsub("\\d", "", models)) - + + model_names <- unique(gsub("\\d", "", names(models))) + getModelFit <- ifelse(simple, getModelFitSimple, getModelFitOpt) - model_fits <- lapply(models, getModelFit, dose_levels, posterior) - + + model_fits <- lapply(model_names, getModelFit, dose_levels, posterior, list("scal" = attr(models, "scal"))) model_fits <- addModelWeights(model_fits) - - names(model_fits) <- models + + names(model_fits) <- model_names attr(model_fits, "posterior") <- posterior class(model_fits) <- "modelFits" - - + + return (model_fits) - + } ## ignoring mixture posterior getModelFitSimple <- function ( - + model, dose_levels, - posterior - + posterior, + addArgs = NULL + ) { - + fit <- DoseFinding::fitMod( dose = dose_levels, resp = summary.postList(posterior)[, 1], @@ -88,11 +91,11 @@ getModelFitSimple <- function ( model = model, type = "general", bnds = DoseFinding::defBnds(mD = max(dose_levels))[[model]]) - + pred_vals <- stats::predict(fit, predType = "ls-means") max_effect <- max(pred_vals) - min(pred_vals) gAIC <- DoseFinding::gAIC(fit) - + model_fit <- list( model = model, coeffs = fit$coefs, @@ -100,93 +103,109 @@ getModelFitSimple <- function ( pred_values = pred_vals, max_effect = max_effect, gAIC = gAIC) - + return (model_fit) - + } ## taking mixture posterior into account getModelFitOpt <- function ( - + model, dose_levels, - posterior - + posterior, + addArgs = NULL + ) { - + getOptParams <- function ( - + model, dose_levels, - posterior - + posterior, + addArgs = NULL + ) { - - switch (model, - "emax" = { - lb <- c(-Inf, -Inf, 0.001 * max(dose_levels)) - ub <- c(Inf, Inf, 1.5 * max(dose_levels)) - expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + (theta[2] * dose_levels) / (theta[3] + dose_levels)))^2 / (post_combs$vars[i, ])))}, - "sigEmax" = { - lb <- c(-Inf, -Inf, 0.001 * max(dose_levels), 0.5) - ub <- c(Inf, Inf, 1.5 * max(dose_levels), 10) - expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + (theta[2] * dose_levels^theta[4]) / (theta[3]^theta[4] + dose_levels^theta[4])))^2 / (post_combs$vars[i, ])))}, - "exponential" = { - lb <- c(-Inf, -Inf, 0.1 * max(dose_levels)) - ub <- c(Inf, Inf, 2 * max(dose_levels)) - expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] * (exp(dose_levels / theta[3]) - 1)))^2 / (post_combs$vars[i, ])))}, - "quadratic" = { - lb <- NULL - ub <- c(Inf, Inf, 0) - expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] * dose_levels + theta[3] * dose_levels^2))^2 / (post_combs$vars[i, ])))}, - "linear" = { - lb <- NULL - ub <- NULL - expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] * dose_levels))^2 / (post_combs$vars[i, ])))}, - "logistic" = { - lb <- c(-Inf, -Inf, 0.001 * max(dose_levels), 0.01 * max(dose_levels)) - ub <- c(Inf, Inf, 1.5 * max(dose_levels), 0.5 * max(dose_levels)) - expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] / (1 + exp((theta[3] - dose_levels) / theta[4]))))^2 / (post_combs$vars[i, ])))}, - { - stop ("error")}) - + switch ( + model, + "emax" = { + lb <- c(-Inf, -Inf, 0.001 * max(dose_levels)) + ub <- c(Inf, Inf, 1.5 * max(dose_levels)) + expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + (theta[2] * dose_levels) / (theta[3] + dose_levels)))^2 / (post_combs$vars[i, ]))) + }, + "sigEmax" = { + lb <- c(-Inf, -Inf, 0.001 * max(dose_levels), 0.5) + ub <- c(Inf, Inf, 1.5 * max(dose_levels), 10) + expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + (theta[2] * dose_levels^theta[4]) / (theta[3]^theta[4] + dose_levels^theta[4])))^2 / (post_combs$vars[i, ]))) + }, + "exponential" = { + lb <- c(-Inf, -Inf, 0.1 * max(dose_levels)) + ub <- c(Inf, Inf, 2 * max(dose_levels)) + expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] * (exp(dose_levels / theta[3]) - 1)))^2 / (post_combs$vars[i, ]))) + }, + "quadratic" = { + lb <- NULL + ub <- c(Inf, Inf, Inf) + expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] * dose_levels + theta[3] * dose_levels^2))^2 / (post_combs$vars[i, ]))) + }, + "linear" = { + lb <- NULL + ub <- NULL + expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] * dose_levels))^2 / (post_combs$vars[i, ]))) + }, + "logistic" = { + lb <- c(-Inf, -Inf, 0.001 * max(dose_levels), 0.01 * max(dose_levels)) + ub <- c(Inf, Inf, 1.5 * max(dose_levels), 0.5 * max(dose_levels)) + expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] / (1 + exp((theta[3] - dose_levels) / theta[4]))))^2 / (post_combs$vars[i, ]))) + }, + "betaMod" = { + lb <- c(-Inf, -Inf, 0.05, 0.05) + ub <- c(Inf, Inf, 4, 4) + scal <- ifelse(is.null(addArgs), 1.2 * max(dose_levels), addArgs[["scal"]]) #for betaMod shape + expr_i <- substitute( + sum( + ((post_combs$means[i, ] - (theta[1] + theta[2] * (((theta[3] + theta[4])^(theta[3] + theta[4])) / ((theta[3] ^ theta[3]) * (theta[4]^theta[4]))) * (dose_levels / scal)^(theta[3]) * ((1 - dose_levels / scal)^(theta[4]))))^2 / (post_combs$vars[i, ]))) + , + list(scal = scal) + ) + }, + stop(paste("error:", model, "shape not yet implemented")) + ) + simple_fit <- getModelFitSimple( model = model, dose_levels = dose_levels, posterior = posterior) - + param_list <- list( expr_i = expr_i, - opts = list("algorithm" = "NLOPT_LN_NELDERMEAD", maxeval = 1e3), + opts = list("algorithm" = "NLOPT_LN_NELDERMEAD", maxeval = 1e3), #,stopval=0 params = list( x0 = simple_fit$coeffs, lb = lb, ub = ub), c_names = names(simple_fit$coeffs)) - + return (param_list) - + } - + optFun <- function ( - theta, - dose_levels, post_combs, expr_i - ) { - + comps <- sapply(seq_along(post_combs$weights), function (i) eval(expr_i)) - + return (sum(post_combs$weights * comps)) - + } - + param_list <- getOptParams(model, dose_levels, posterior) post_combs <- getPostCombsI(posterior) - + fit <- nloptr::nloptr( eval_f = optFun, x0 = param_list$params$x0, @@ -195,36 +214,34 @@ getModelFitOpt <- function ( opts = param_list$opts, dose_levels = dose_levels, post_combs = post_combs, - expr_i = param_list$expr_i) - + expr_i = param_list$expr_i + ) + names(fit$solution) <- param_list$c_names - + model_fit <- list( model = model, coeffs = fit$solution, dose_levels = dose_levels) - - model_fit$pred_values <- predictModelFit(model_fit) + + model_fit$pred_values <- predictModelFit(model_fit = model_fit, doses = model_fit$dose_levels, addArgs = addArgs) model_fit$max_effect <- max(model_fit$pred_values) - min(model_fit$pred_values) model_fit$gAIC <- getGenAIC(model_fit, post_combs) - + return (model_fit) - + } predictModelFit <- function ( - + model_fit, - doses = NULL - + doses = NULL, + addArgs = NULL + ) { - if (is.null(doses)) { - doses <- model_fit$dose_levels - } - pred_vals <- switch ( model_fit$model, "emax" = {DoseFinding::emax( @@ -258,46 +275,52 @@ predictModelFit <- function ( model_fit$coeffs["eMax"], model_fit$coeffs["ed50"], model_fit$coeffs["delta"])}, - {stop("error")}) - + "betaMod" = {DoseFinding::betaMod( + doses, + model_fit$coeffs["e0"], + model_fit$coeffs["eMax"], + model_fit$coeffs["delta1"], + model_fit$coeffs["delta2"], + ifelse(is.null(addArgs) | is.null(addArgs[["scal"]]), 1.2 * max(doses), addArgs[["scal"]]))}, + {stop(paste("error:", model_fit$model, "shape not yet implemented"))}) + return (pred_vals) - + } getGenAIC <- function ( - + model_fit, post_combs - + ) { - + expr_i <- quote(sum(1 / post_combs$vars[i, ] * (post_combs$means[i, ] - model_fit$pred_values)^2) + 2 * length(model_fit$coeffs)) comb_aic <- sapply(seq_along(post_combs$weights), function (i) eval(expr_i)) - + g_aic <- stats::weighted.mean(comb_aic, w = post_combs$weights) - + return (g_aic) - + } addModelWeights <- function ( - + model_fits - + ) { - g_aics <- sapply(model_fits, function (model_fit) model_fit$gAIC) exp_values <- exp(-0.5 * g_aics) model_weights <- exp_values / sum(exp_values) - + names(model_weights) <- NULL - + model_fits_out <- lapply(seq_along(model_fits), function (i) { - + c(model_fits[[i]], model_weight = model_weights[i]) - + }) - + return (model_fits_out) - -} \ No newline at end of file + +} diff --git a/R/plot.R b/R/plot.R index 0486d97..bd2d9d8 100644 --- a/R/plot.R +++ b/R/plot.R @@ -1,11 +1,11 @@ #' @title plot.modelFits -#' +#' #' @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). +#' 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 (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. +#' The default setting is that these credible bands are not calculated. #' @param x An object of type modelFits #' @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 @@ -19,7 +19,7 @@ #' @examples #' 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_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("exponential", "linear") @@ -28,12 +28,12 @@ #' posterior = posterior_list, #' dose_levels = dose_levels, #' simple = TRUE) -#' -#' plot(fit) +#' +#' plot(fit) #' @return A ggplot2 object #' @export plot.modelFits <- function ( - + x, gAIC = TRUE, avg_fit = TRUE, @@ -44,12 +44,11 @@ plot.modelFits <- function ( n_bs_smpl = 1e3, acc_color = "orange", ... - + ) { - ## R CMD --as-cran appeasement .data <- NULL - + checkmate::check_class(x, "modelFits") checkmate::check_logical(gAIC) checkmate::check_logical(avg_fit) @@ -59,10 +58,10 @@ plot.modelFits <- function ( checkmate::check_double(alpha_CrB, lower = 0, upper = 1, len = 2) checkmate::check_integer(n_bs_smpl, lower = 1, upper = Inf) checkmate::check_string(acc_color, na.ok = TRUE) - + plot_res <- 1e2 model_fits <- x - + dose_levels <- model_fits[[1]]$dose_levels post_summary <- summary.postList( object = attr(model_fits, "posterior"), @@ -70,82 +69,77 @@ plot.modelFits <- function ( dose_seq <- seq(from = min(dose_levels), to = max(dose_levels), length.out = plot_res) - + preds_models <- sapply(model_fits, predictModelFit, doses = dose_seq) model_names <- names(model_fits) - + if (avg_fit) { - + mod_weights <- sapply(model_fits, function (x) x$model_weight) avg_preds <- preds_models %*% mod_weights - + preds_models <- cbind(preds_models, avg_preds) model_names <- c(model_names, "avgFit") - + } - + gg_data <- data.frame( 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")), - each = plot_res)) - + models = rep(factor(model_names, levels = model_names), + each = plot_res)) + if (gAIC) { - + g_AICs <- sapply(model_fits, function (x) x$gAIC) label_gAUC <- paste("AIC:", round(g_AICs, digits = 1)) - + if (avg_fit) { - + mod_weights <- sort(mod_weights, decreasing = TRUE) paste_names <- names(mod_weights) |> gsub("exponential", "exp", x = _) |> gsub("quadratic", "quad", x = _) |> gsub("linear", "lin", x = _) |> gsub("logistic", "log", x = _) |> - gsub("sigEmax", "sigE", x = _) - + gsub("sigEmax", "sigE", x = _) |> + gsub("betaMod", "betaM", x = _)|> + gsub("quadratic", "quad", x = _) + label_avg <- paste0(paste_names, "=", round(mod_weights, 1), collapse = ", ") label_gAUC <- c(label_gAUC, label_avg) - + } - + } - + if (cr_bands) { - - bs_samples <- getBootstrapSamples( + + crB_data <- getBootstrapQuantiles( model_fits = model_fits, n_samples = n_bs_smpl, + quantiles = c(0.5, sort(unique(c(alpha_CrB / 2, 1 - alpha_CrB / 2)))), avg_fit = avg_fit, doses = dose_seq) - - crB_data <- getBootstrapQuantiles( - bs_samples = bs_samples, - quantiles = c(0.5, sort(unique(c(alpha_CrB / 2, 1 - alpha_CrB / 2))))) - colnames(crB_data) <- c("models", colnames(crB_data)[-1]) - + getInx <- function (alpha_CrB) { n <- length(alpha_CrB) - inx_list <- lapply(seq_len(n), function (i) c(i, 2 * n - i + 2) + 2) + inx_list <- lapply(seq_len(n), function (i) c(i, 2 * n - i + 1) + 3) return (inx_list)} - + } - - plts <- ggplot2::ggplot() + + + plts <- ggplot2::ggplot() + ## Layout etc. - ggplot2::theme_bw() + + ggplot2::theme_bw() + ggplot2::labs(x = "Dose", y = "Model Fits") + ggplot2::theme(panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank()) + ## gAIC {if (gAIC) { - + ggplot2::geom_text( data = data.frame( models = unique(gg_data$models), @@ -153,42 +147,42 @@ plot.modelFits <- function ( mapping = ggplot2::aes(label = label_gAUC), x = -Inf, y = Inf, hjust = "inward", vjust = "inward", size = 3) - - } + } - + } + if (cr_bands) { - + ## Bootstrapped Credible Bands for (inx in getInx(alpha_CrB)) { - + loop_txt <- paste0( "ggplot2::geom_ribbon( data = crB_data, - mapping = ggplot2::aes(x = .data$dose, + mapping = ggplot2::aes(x = doses, ymin = crB_data[, ", inx[1], "], ymax = crB_data[, ", inx[2], "]), - fill = acc_color, + fill = acc_color, alpha = 0.25)") - + plts <- plts + eval(parse(text = loop_txt)) - + } rm(inx) - + ## Bootstrapped Fit plts <- plts + ggplot2::geom_line( data = crB_data, - mapping = ggplot2::aes(.data$dose, .data$`50%`), + mapping = ggplot2::aes(.data$doses, .data$`50%`), color = acc_color) - + } - + plts <- plts + ## Posterior Credible Intervals {if (cr_intv) { - + ggplot2::geom_errorbar( data = data.frame(x = dose_levels, ymin = post_summary[, 3], @@ -198,9 +192,9 @@ plot.modelFits <- function ( ymax = .data$ymax), width = 0, alpha = 0.5) - - } - } + + + } + } + ## Posterior Medians ggplot2::geom_point( data = data.frame( @@ -211,10 +205,10 @@ plot.modelFits <- function ( ## Fitted Models ggplot2::geom_line( data = gg_data, - mapping = ggplot2::aes(.data$doses, .data$fits)) + + mapping = ggplot2::aes(.data$doses, .data$fits)) + ## Faceting ggplot2::facet_wrap(~ models) - + return (plts) - -} \ No newline at end of file + +} diff --git a/R/s3methods.R b/R/s3methods.R index e9dca86..c714a0d 100644 --- a/R/s3methods.R +++ b/R/s3methods.R @@ -2,75 +2,113 @@ #' @export print.BayesianMCPMod <- function ( - + x, ... - + ) { - + print(x$BayesianMCP) - + } ## BayesianMCP -------------------------------------------- #' @export -print.BayesianMCP <- function ( - - x, - ... - -) { - +print.BayesianMCP <- function(x, ...) { + cat("Bayesian Multiple Comparison Procedure\n\n") + + cat("Summary:\n") + cat(" Sign:", x[1, "sign"], "\n") + cat(" Critical Probability:", x[1, "crit_prob_adj"], "\n") + cat(" Maximum Posterior Probability:", x[1, "max_post_prob"], "\n\n") n_sim <- nrow(x) - + cat("Effective Sample Size (ESS) per Dose Group:\n") + print(attr(x, "ess_avg"), row.names = FALSE) + cat("\n") + + cat("Posterior Probabilities for Model Shapes:\n") + model_probs <- x[1, grep("^post_probs\\.", colnames(x))] + model_names <- gsub("post_probs\\.", "", names(model_probs)) + model_df <- data.frame(Model = model_names, Probability = unlist(model_probs)) + print(model_df, row.names = FALSE) + cat("Bayesian Multiple Comparison Procedure\n") - + if (n_sim == 1L) { - + attr(x, "critProbAdj") <- NULL attr(x, "successRate") <- NULL class(x) <- NULL - + print.default(x, ...) - + # if (any(!is.na(attr(x, "essAvg")))) { - # + # # cat("Average Posterior ESS\n") # print(attr(x, "essAvg"), ...) - # + # # } - + } else { - + model_successes <- getModelSuccesses(x) - + cat(" Estimated Success Rate: ", attr(x, "successRate"), "\n") cat(" N Simulations: ", n_sim) - + cat("\n") cat("Model Significance Frequencies\n") print(model_successes, ...) - + } - + + invisible(x) } +# print.BayesianMCP <- function ( + # +# x, +# ... +# +# ) { +# +# n_sim <- nrow(x) +# +# cat("Bayesian Multiple Comparison Procedure\n") +# +# 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: ", attr(x, "successRate"), "\n") +# cat(" N Simulations: ", n_sim) +# +# } +# +# } + ## ModelFits ---------------------------------------------- #' @title predict.modelFits #' @description This function performs model predictions based on the provided -#' model and dose specifications -#' +#' 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 +#' done #' @param ... Currently without function #' @examples #' 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_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") @@ -78,56 +116,56 @@ print.BayesianMCP <- function ( #' 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 ( - + object, doses = NULL, ... - + ) { predictions <- lapply(object, predictModelFit, doses = doses) attr(predictions, "doses") <- doses - + return (predictions) - + } #' @export print.modelFits <- function ( - + x, ... - + ) { - + n_digits = 1 - + dose_levels <- x[[1]]$dose_levels dose_names <- names(attr(x, "posterior")) - + predictions <- t(sapply(x, function (y) y$pred_values)) colnames(predictions) <- dose_names - + out_table <- data.frame(predictions, mEff = sapply(x, function (y) y$max_effect), 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)) { - + out_table <- as.data.frame(out_table) out_table$Sign <- sapply(x, function (y) y$significant) - + } - + model_names <- names(x) |> gsub("exponential", "exponential", x = _) |> gsub("quadratic", "quadratic ", x = _) |> @@ -135,18 +173,18 @@ print.modelFits <- function ( gsub("logistic", "logistic ", x = _) |> gsub("emax", "emax ", x = _) |> gsub("sigEmax", "sigEmax ", x = _) - + # cat("Bayesian MCP Model Fits\n\n") cat("Model Coefficients\n") for (i in seq_along(model_names)) { - + coeff_values <- x[[i]]$coeff coeff_names <- names(coeff_values) - + cat(model_names[i], paste(coeff_names, round(coeff_values, n_digits), sep = " = "), "\n", sep = " ") - + } cat("\n") cat("Dose Levels\n", @@ -154,7 +192,7 @@ print.modelFits <- function ( cat("\n") cat("Predictions, Maximum Effect, gAIC, Model Weights & Significance\n") print(out_table, ...) - + } ## plot.ModelFits() @@ -164,57 +202,57 @@ print.modelFits <- function ( #' @export summary.postList <- function ( - + object, ... - + ) { - + summary_list <- lapply(object, summary, ...) names(summary_list) <- names(object) summary_tab <- do.call(rbind, summary_list) - + return (summary_tab) - + } #' @export print.postList <- function ( - + x, ... - + ) { - + getMaxDiff <- function ( - + medians - + ) { - + diffs <- medians - medians[1] - + max_diff <- max(diffs) max_diff_level <- which.max(diffs) - 1 - + out <- c(max_diff, max_diff_level) names(out) <- c("max_diff", "DG") - + return (out) - + } - + summary_tab <- summary.postList(x) - + names(x) <- rownames(summary_tab) class(x) <- NULL - + list_out <- list(summary_tab, getMaxDiff(summary_tab[, 4]), x) names(list_out) <- c("Summary of Posterior Distributions", "Maximum Difference to Control and Dose Group", "Posterior Distributions") - + print(list_out, ...) - + } diff --git a/_pkgdown.yml b/_pkgdown.yml index 1b0d774..934cedf 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -12,4 +12,5 @@ articles: navbar: ~ contents: - analysis_normal + - analysis_normal_noninformative - Simulation_Example diff --git a/inst/WORDLIST b/inst/WORDLIST index 2bded23..af9ba59 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,5 +1,59 @@ +Bjoern +Bornkamp +Bretz +CRC +Codecov +Dette +DoseFinding +GLS +Gierse +Iasonos +Kunzmann +LinLog +Loley +MCP MCPMod +Macos +NCT +Nelder +O'Quigley Pinheiro +RBesT Schmidli +Schorning +Thoma +Thomann +Un +Yao al -et \ No newline at end of file +assessDesign +betaMod +biom +critVal +doi +emax +ess +et +gAIC +getBootstrapQuantiles +getContr +getCritProb +getESS +getModelFits +getPosterior +ggplot +mixnorm +modelFits +nloptr +normMix +optContr +performBayesianMCP +performBayesianMCPMod +postList +postmix +pst +rmix +sd +se +simulateData +summand diff --git a/man/getBootstrapQuantiles.Rd b/man/getBootstrapQuantiles.Rd index 82d5889..4e22634 100644 --- a/man/getBootstrapQuantiles.Rd +++ b/man/getBootstrapQuantiles.Rd @@ -4,24 +4,39 @@ \alias{getBootstrapQuantiles} \title{getBootstrapQuantiles} \usage{ -getBootstrapQuantiles(bs_samples, quantiles) +getBootstrapQuantiles( + model_fits, + quantiles, + n_samples = 1000, + doses = NULL, + avg_fit = TRUE +) } \arguments{ -\item{bs_samples}{An object of class bootstrappedSample as created by getBootstrapSamples} +\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}{A vector of quantiles that should be evaluated} + +\item{n_samples}{Number of samples that should be drawn as basis for the bootstrapped quantiles} + +\item{doses}{A vector of doses for which a prediction should be performed} + +\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{ -A data frame with entries doses, models, and quantiles +A data frame with columns for model, dose, and bootstrapped samples } \description{ -Calculates quantiles from bootstrapped dose predictions. -Can be used to derive credible intervals to assess the uncertainty for the model fit. +A function for the calculation of bootstrapped model predictions. +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). +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{ 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_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("exponential", "linear") @@ -30,11 +45,13 @@ fit <- getModelFits(models = models, posterior = posterior_list, dose_levels = dose_levels, simple = TRUE) - -bs_samples <- getBootstrapSamples(model_fits = fit, - n_samples = 10, # speeding up example run time - doses = c(0, 6, 8)) - -getBootstrapQuantiles(bs_samples = bs_samples, - quantiles = c(0.025, 0.5, 0.975)) + +getBootstrapQuantiles(model_fits = fit, + quantiles = c(0.025, 0.5, 0.975), + n_samples = 10, # speeding up example run time + doses = c(0, 6, 8)) + +} +\references{ +O'Quigley J, Iasonos A, Bornkamp B. 2017. Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials (1st ed.). Chapman and Hall/CRC. doi:10.1201/9781315151984 } diff --git a/man/getBootstrapSamples.Rd b/man/getBootstrapSamples.Rd deleted file mode 100644 index e81068e..0000000 --- a/man/getBootstrapSamples.Rd +++ /dev/null @@ -1,48 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/bootstrapping.R -\name{getBootstrapSamples} -\alias{getBootstrapSamples} -\title{getBootstrapSamples} -\usage{ -getBootstrapSamples(model_fits, n_samples = 1000, doses = NULL, avg_fit = TRUE) -} -\arguments{ -\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{n_samples}{Number of samples that should be drawn as basis for the bootstrapped quantiles} - -\item{doses}{A vector of doses for which a prediction should be performed} - -\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{ -A data frame with columns for model, dose, and bootstrapped samples -} -\description{ -A function for the calculation of bootstrapped model predictions. -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). -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{ -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("exponential", "linear") -dose_levels <- c(0, 1, 2, 4, 8) -fit <- getModelFits(models = models, - posterior = posterior_list, - dose_levels = dose_levels, - simple = TRUE) - -getBootstrapSamples(model_fits = fit, - n_samples = 10, # speeding up example run time - doses = c(0, 6, 8)) - -} -\references{ -O'Quigley J, Iasonos A, Bornkamp B. 2017. Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials (1st ed.). Chapman and Hall/CRC. doi:10.1201/9781315151984 -} diff --git a/man/getModelFits.Rd b/man/getModelFits.Rd index 3c7eab2..0f37483 100644 --- a/man/getModelFits.Rd +++ b/man/getModelFits.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/modelling.R +% Please edit documentation in R/modeling.R \name{getModelFits} \alias{getModelFits} \title{getModelFits} @@ -40,7 +40,7 @@ where \eqn{Q} denotes the number of models included in the averaging procedure. \examples{ 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_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") @@ -53,7 +53,7 @@ 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/performBayesianMCPMod.Rd b/man/performBayesianMCPMod.Rd index bac1a7f..d2fb698 100644 --- a/man/performBayesianMCPMod.Rd +++ b/man/performBayesianMCPMod.Rd @@ -16,7 +16,7 @@ performBayesianMCPMod(posterior_list, contr, crit_prob_adj, simple = FALSE) \item{simple}{Boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits() function. Default FALSE.} } \value{ -Bayesian MCP test result as well as modelling result. +Bayesian MCP test result as well as modeling result. } \description{ Performs Bayesian MCP Test step and modeling in a combined fashion. See performBayesianMCP() function for MCP Test step and getModelFits() for the modeling step diff --git a/man/plot.modelFits.Rd b/man/plot.modelFits.Rd index cbd372a..53e9256 100644 --- a/man/plot.modelFits.Rd +++ b/man/plot.modelFits.Rd @@ -52,7 +52,7 @@ The default setting is that these credible bands are not calculated. \examples{ 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_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("exponential", "linear") @@ -61,6 +61,6 @@ fit <- getModelFits(models = models, posterior = posterior_list, dose_levels = dose_levels, simple = TRUE) - -plot(fit) + +plot(fit) } diff --git a/man/predict.modelFits.Rd b/man/predict.modelFits.Rd index b819293..fc3c835 100644 --- a/man/predict.modelFits.Rd +++ b/man/predict.modelFits.Rd @@ -25,7 +25,7 @@ model and dose specifications \examples{ 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_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") @@ -33,7 +33,7 @@ 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/tests/testthat/test-bootstrapping.R b/tests/testthat/test-bootstrapping.R index 7ad8ae8..a6a2a7e 100644 --- a/tests/testthat/test-bootstrapping.R +++ b/tests/testthat/test-bootstrapping.R @@ -1,37 +1,36 @@ # Test cases test_that("test getBootstrapBands", { - + BMCP_result <- performBayesianMCP( posterior_list = posterior_list, contr = contr_mat, crit_prob_adj = crit_pval) model_shapes <- c("linear", "emax", "exponential") - + # Option b) Making use of the complete posterior distribution fit <- getModelFits( models = model_shapes, dose_levels = dose_levels, posterior = posterior_list, simple = FALSE) - - bs_samples <- getBootstrapSamples(model_fits = fit, - doses = c(0, 1,2,3,4,6,8)) - - result <- getBootstrapQuantiles(bs_samples, quantiles = c(0.025,0.5, 0.975)) + + result <- bootstrap_quantiles <- BayesianMCPMod::getBootstrapQuantiles( + model_fits = fit, + quantiles = c(0.025, 0.5, 0.975), + doses = c(0, 2.5, 4, 5, 7, 10), + n_samples = 6 + ) + expect_type(result, "list") - - bs_samples_2 <- getBootstrapSamples(model_fits = fit, - n_samples = 1e2, - avg_fit = FALSE, - doses = c(1, 2, 3)) - - result_2 <- getBootstrapQuantiles(bs_samples_2, quantiles = c(0.1, 0.9)) - expect_type(result_2, "list") - - bs_samples_3 <- getBootstrapSamples(model_fits = fit, - doses = NULL) - result_3 <- getBootstrapQuantiles(bs_samples_3, quantiles = c(0.025,0.5, 0.975)) - expect_type(result_3, "list") -}) \ No newline at end of file + result_2 <- bootstrap_quantiles <- BayesianMCPMod::getBootstrapQuantiles( + model_fits = fit, + quantiles = c(0.025, 0.5, 0.975), + avg_fit = FALSE, + doses = c(0, 2.5, 4, 5, 7, 10), + n_samples = 6 + ) + + expect_type(result_2, "list") +}) diff --git a/vignettes/.gitignore b/vignettes/.gitignore index 097b241..5d2fa41 100644 --- a/vignettes/.gitignore +++ b/vignettes/.gitignore @@ -1,2 +1,6 @@ *.html *.R +*_cache/* +*_files/* +*archive +*_files diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index af6159f..d5d2323 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -264,14 +264,12 @@ plot(fit, cr_bands = TRUE) The bootstrap based quantiles can also be directly calculated via the getBootstrapQuantiles() function. For this example, only 6 quantiles are bootstrapped for each model fit. ```{r Bootstrap} -bs_sample <- getBootstrapSamples( +getBootstrapQuantiles( model_fits = fit, + quantiles = c(0.025, 0.5, 0.975), doses = c(0, 2.5, 4, 5, 7, 10), - n_samples = 6) - -getBootstrapQuantiles( - bs_sample = bs_sample, - quantiles = c(0.025, 0.5, 0.975)) + n_samples = 6 +) ``` 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. The main fit, i.e. the black lines in the plot, show the best fit of a certain model based on minimizing the residuals for the posterior distribution, while the bootstrap based 50% quantile shows the median fit of the random sampling and fitting procedure. diff --git a/vignettes/analysis_normal_noninformative.qmd b/vignettes/analysis_normal_noninformative.qmd new file mode 100644 index 0000000..f1b5659 --- /dev/null +++ b/vignettes/analysis_normal_noninformative.qmd @@ -0,0 +1,413 @@ +--- +title: "Analysis Example of Bayesian MCPMod for Continuous Data" +format: + html: + fig-height: 3.5 + self-contained: true + toc: true + number-sections: true + bibliography: references.bib +vignette: > + %\VignetteIndexEntry{Analysis Example of Bayesian MCPMod for Continuous Data} + %\VignetteEncoding{UTF-8} + %\VignetteEngine{quarto::html} +--- + +```{r setup} +library(BayesianMCPMod) +library(RBesT) +library(clinDR) +library(dplyr) +library(tibble) +library(reactable) + +set.seed(7015) +``` + +```{r} +#| code-summary: setup +#| code-fold: true +#| message: false +#| warning: false + +#' Display Parameters Table +#' +#' This function generates a markdown table displaying the names and values of parameters +#' from a named list. +#' +#' @param named_list A named list where each name represents a parameter name and the list +#' element represents the parameter value. Date values in the list are automatically +#' converted to character strings for display purposes. +#' +#' @return Prints a markdown table with two columns: "Parameter Name" and "Parameter Values". +#' The function does not return a value but displays the table directly to the output. +#' +#' @importFrom knitr kable +#' @examples +#' params <- list("Start Date" = as.Date("2020-01-01"), +#' "End Date" = as.Date("2020-12-31"), +#' "Threshold" = 10) +#' display_params_table(params) +#' +#' @export +display_params_table <- function(named_list) { + display_table <- data.frame() + value_names <- data.frame() + for (i in 1:length(named_list)) { + # dates will display as numeric by default, so convert to char first + if (class(named_list[[i]]) == "Date") { + named_list[[i]] = as.character(named_list[[i]]) + } + if (!is.null(names(named_list[[i]]))) { + value_names <- rbind(value_names, paste(names(named_list[[i]]), collapse = ', ')) + } + values <- data.frame(I(list(named_list[[i]]))) + display_table <- rbind(display_table, values) + } + + round_numeric <- function(x, digits = 3) { + if (is.numeric(x)) { + return(round(x, digits)) + } else { + return(x) + } + } + + display_table[1] <- lapply(display_table[1], function(sublist) { + lapply(sublist, round_numeric) + }) + + class(display_table[[1]]) <- "list" + + if (nrow(value_names) == 0) { + knitr::kable( + cbind(names(named_list), display_table), + col.names = c("Name", "Value") + ) + } else { + knitr::kable( + cbind(names(named_list), value_names, display_table), + col.names = c("Name", "Value Labels", "Value") + ) + } +} +``` + +# Introduction + +This vignette demonstrates the application of the {BayesianMCPMod} package for +analyzing a phase 2 dose-finding trial using the Bayesian MCPMod approach. + +# Prior Specification + +Ideally, priors are grounded in historical data. This approach allows for the synthesis of +prior knowledge with current data, enhancing the accuracy of trial evaluations. + +The focus of this vignette is more generic, however. We specify weakly-informative +priors across all dose groups to allow the trial data to have a stronger influence on the analysis. + +```{r} +dose_levels <- c(0, 2.5, 5, 10) + +prior_list <- lapply(dose_levels, function(dose_group) { + RBesT::mixnorm(weak = c(w = 1, m = 0, s = 200), sigma = 10) +}) + +names(prior_list) <- c("Ctr", paste0("DG_", dose_levels[-1])) +``` + +# Dose-Response Model Shapes + +Candidate models are specified using the {DoseFinding} package. Models can be +parameterized using guesstimates or by directly providing distribution parameters. +Note that the linear candidate model does not require parameterization. + +[**Note:** The LinLog model is rarely used and not currently supported by `{BayesianMCPMod}`.]{.aside} + +In the code below, the models are "guesstimated" using the `DoseFinding::guesst` function. +The `d` option usually takes a single value (a dose level), and the corresponding `p` +for the maximum effect achieved at `d`. + + +```{r} +# Guesstimate estimation +exp_guesst <- DoseFinding::guesst( + model = "exponential", + d = 5, p = 0.2, Maxd = max(dose_levels) +) +emax_guesst <- DoseFinding::guesst( + model = "emax", + d = 2.5, p = 0.9 +) +sigEmax_guesst <- DoseFinding::guesst( + model = "sigEmax", + d = c(2.5, 5), p = c(0.5, 0.95) +) +logistic_guesst <- DoseFinding::guesst( + model = "logistic", + d = c(5, 10), p = c(0.1, 0.85) +) +``` + +In some cases, you need to provide more information. For instance, `sigEmax` +requires a pair of `d` and `p` values, and `exponential` requires the specification of +the maximum dose for the trial (`Maxd`). + +[See the help files for model specifications by typing `?DoseFinding::guesst` in your console]{.aside} + + + +Of course, you can also specify the models directly on the parameter scale (without using `DoseFinding::guesst`). + +For example, you can get a betaMod model by specifying `delta1` and `delta2` +parameters (`scale` is assumed to be `1.2` of the maximum dose), or a quadratic model with the `delta2` parameter. + +```{r} +betaMod_params <- c(delta1 = 1, delta2 = 1) +quadratic_params <- c(delta2 = -0.1) +``` + +Now, we can go ahead and create a `Mods` object, which will be used in the remainder +of the vignette. + +```{r} +mods <- DoseFinding::Mods( + linear = NULL, + # guesstimate scale + exponential = exp_guesst, + emax = emax_guesst, + sigEmax = sigEmax_guesst, + logistic = logistic_guesst, + # parameter scale + betaMod = betaMod_params, + quadratic = quadratic_params, + # Options for all models + doses = dose_levels, + maxEff = -1, + placEff = -12.8 +) + +plot(mods) +``` + +The `mods` object we just created above contains the full model parameters, which can be helpful for +understanding how the guesstimates are translated onto the parameter scale. + +```{r} +display_params_table(mods) +``` + +And we can see the assumed treatment effects for the specified dose groups below: + +```{r} +knitr::kable(DoseFinding::getResp(mods, doses = dose_levels)) +``` + +## Trial Data + +We will use the trial with ct.gov number NCT00735709 as our phase 2 trial data, +available in the `{clinDR}` package [@nct00735709_2024a]. + +```{r} +data("metaData") + +trial_data <- dplyr::filter( + dplyr::filter(tibble::tibble(metaData), bname == "BRINTELLIX"), + primtime == 8, + indication == "MAJOR DEPRESSIVE DISORDER", + protid == 5 +) + +n_patients <- c(128, 124, 129, 122) +``` + +# Posterior Calculation + +In the first step of Bayesian MCPMod, the posterior is calculated by combining +the prior information with the estimated results of the trial [@fleischer_2022]. + +```{r} +posterior <- getPosterior( + prior_list = prior_list, + mu_hat = trial_data$rslt, + se_hat = trial_data$se, + calc_ess = TRUE +) + +knitr::kable(summary(posterior)) +``` + +# Bayesian MCPMod Test Step + +The testing step of Bayesian MCPMod is executed using a critical value on the +probability scale and a pseudo-optimal contrast matrix. + +The critical value is calculated using (re-estimated) contrasts for frequentist +MCPMod to ensure error control when using weakly-informative priors. + +A pseudo-optimal contrast matrix is generated based on the variability of the +posterior distribution (see [@fleischer_2022] for more details). + +```{r} +crit_pval <- getCritProb( + mods = mods, + dose_levels = dose_levels, + se_new_trial = trial_data$se, + alpha_crit_val = 0.05 +) + +contr_mat <- getContr( + mods = mods, + dose_levels = dose_levels, + sd_posterior = summary(posterior)[, 2] +) +``` + +Please note that there are different ways to derive the contrasts. +The following code shows the implementation of some of these ways but it is not +executed and the contrast specification above is used. + +```{r} +#| eval: false + +# 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 = trial_data$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) +``` + +The Bayesian MCP testing step is then executed: + +```{r} +BMCP_result <- performBayesianMCP( + posterior_list = posterior, + contr = contr_mat, + crit_prob_adj = crit_pval) +``` + +Summary information: + +```{r} +BMCP_result +``` + +The testing step is significant, indicating a non-flat dose-response shape. +All models are significant, with the `emax` model indicating the greatest deviation +from the null hypothesis. + +# Model Fitting and Visualization + +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. In contrast, for the full fit, the +non-linear optimization problem is addressed via the Nelder-Mead algorithm +[@neldermead_2024a] implemented by the `{nloptr}` package. + +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, so we +present the full fit. + +```{r} +# If simple = TRUE, uses approx posterior +# Here we use complete posterior distribution +fit <- getModelFits( + models = mods, + dose_levels = dose_levels, + posterior = posterior, + simple = FALSE) +``` + +Estimates for dose levels not included in the trial: + +```{r} +display_params_table(stats::predict(fit, doses = c(0, 2.5, 4, 5, 7, 10))) +``` + +Plots of fitted dose-response models and an AIC-based average model: + +```{r} +plot(fit) +``` + +To assess the uncertainty, one can additionally visualize credible bands +(orange shaded areas, default levels are 50% and 95%). + +These credible bands are calculated with a bootstrap method as follows: + +- Samples from the posterior distribution are drawn and for every sample the +simplified fitting step and a prediction is performed. + +- These predictions are then used to identify and visualize the specified quantiles. + +```{r} +plot(fit, cr_bands = TRUE) +``` + +The bootstrap based quantiles can also be directly calculated via the +`getBootstrapQuantiles()` function. + +For this example, only 6 quantiles are bootstrapped for each model fit. + +```{r} +bootstrap_quantiles <- getBootstrapQuantiles( + model_fits = fit, + quantiles = c(0.025, 0.5, 0.975), + doses = c(0, 2.5, 4, 5, 7, 10), + n_samples = 6 +) +``` + +```{r} +#| code-fold: true +reactable::reactable( + data = bootstrap_quantiles, + groupBy = "models", + columns = list( + doses = colDef(aggregate = "count", format = list(aggregated = colFormat(suffix = " doses"))), + "2.5%" = colDef(aggregate = "mean", format = list(aggregated = colFormat(prefix = "mean = ", digits = 2), cell = colFormat(digits = 4))), + "50%" = colDef(aggregate = "mean", format = list(aggregated = colFormat(prefix = "mean = ", digits = 2), cell = colFormat(digits = 4))), + "97.5%" = colDef(aggregate = "mean", format = list(aggregated = colFormat(prefix = "mean = ", digits = 2), cell = colFormat(digits = 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. + +The main fit (black line) minimizes residuals for the posterior distribution, +while the bootstrap median is the median fit of random sampling. + +# Additional note + +Testing and modeling can also be combined via `performBayesianMCPMod()`, +but this is not run here. + +```{r} +#| eval: false +performBayesianMCPMod( + posterior_list = posterior, + contr = contr_mat, + crit_prob_adj = crit_pval, + simple = FALSE) +``` diff --git a/vignettes/references.bib b/vignettes/references.bib new file mode 100644 index 0000000..d6195b7 --- /dev/null +++ b/vignettes/references.bib @@ -0,0 +1,25 @@ +@article{fleischer_2022, + title={Bayesian MCPMod}, + author={Fleischer F, Bossert S, Deng Q, Loley C, Gierse J}, + journal={Pharmaceutical Statistics}, + volume={21}, + number={3}, + pages={654--670}, + year={2022} +} + +@misc {neldermead_2024a, + author = "Wikipedia", + title = "Nelder-Mead method", + year = "2024", + url = "https://en.wikipedia.org/wiki/Nelder-Mead_method", + note = " [Online; accessed 25-February-2024]" +} + +@misc {nct00735709_2024a, + author = "ClinicalTrials.gov", + title = "NCT00735709", + year = "2024", + url = "https://clinicaltrials.gov/study/NCT00735709?term=NCT00735709&rank=1", + note = " [Online; accessed 25-February-2024]" +}