diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml new file mode 100644 index 0000000..fab4029 --- /dev/null +++ b/.github/workflows/test-coverage.yaml @@ -0,0 +1,38 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/master/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +on: + push: + branches: [main] + pull_request: + branches: [main] + workflow_dispatch: + +name: test-coverage + +jobs: + test-coverage: + runs-on: ubuntu-latest + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + RENV_CONFIG_SANDBOX_ENABLED: false + + steps: + - uses: actions/checkout@v2 + + - uses: r-lib/actions/setup-pandoc@v2 + + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: covr + + - name: Install BayesianMCPMod + shell: bash + run: R CMD INSTALL --preclean . + + - name: Test coverage + run: covr::codecov() + shell: Rscript {0} \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 034bfc9..8744241 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BayesianMCPMod Title: Bayesian MCPMod -Version: 0.1.3 +Version: 0.1.3-2 Authors@R: c( person("Sebastian", "Bossert", role = "aut", @@ -22,11 +22,19 @@ Roxygen: list(markdown = TRUE) RoxygenNote: 7.2.3 Imports: DoseFinding, + ggplot2, stats, RBesT, - nloptr + nloptr, + clinDR, + knitr, + rmarkdown, + dplyr Suggests: testthat (>= 3.0.0) -Config/testthat/edition: 3 Depends: R (>= 4.1) +VignetteBuilder: knitr +Config/testthat/edition: 3 +URL: https://github.com/Boehringer-Ingelheim/BayesianMCPMod +BugReports: https://github.com/Boehringer-Ingelheim/BayesianMCPMod/issues diff --git a/NAMESPACE b/NAMESPACE index 2e76e4c..f89ff4d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,8 +1,19 @@ # Generated by roxygen2: do not edit by hand -export(BMCPMod) -export(doFit) -export(estimateModels) -export(getGenAICs) -export(postShape) +S3method(plot,modelFits) +S3method(predict,ModelFits) +S3method(print,BayesianMCP) +S3method(print,BayesianMCPMod) +S3method(print,modelFits) +S3method(print,postList) +S3method(summary,postList) +export(assessDesign) +export(getBootstrapBands) +export(getContrMat) +export(getCritProb) +export(getModelFits) +export(getPosterior) +export(getPriorList) +export(performBayesianMCP) +export(performBayesianMCPMod) export(simulateData) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index c9215cf..0a3f56c 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -1,26 +1,277 @@ -#' @title BMCPMod -#' @param ancova1 tbd -#' -#' @param cont_Mat1 tbd +#' @title assessDesign +#' +#' @param n_patients tbd +#' @param mods tbd +#' @param prior_list tbd +#' @param n_sim tbd +#' @param alpha_crit_val tbd +#' @param simple tbd +#' +#' @export +assessDesign <- function ( + + n_patients, + mods, + prior_list, + + n_sim = 1e3, + alpha_crit_val = 0.05, + simple = TRUE + +) { + + dose_levels <- attr(prior_list, "dose_levels") + + data <- simulateData( + n_patients = n_patients, + dose_levels = dose_levels, + sd = attr(prior_list, "sd_tot"), + mods = mods, + n_sim = n_sim) + + model_names <- names(mods) + + eval_design <- lapply(model_names, function (model_name) { + + posterior_list <- getPosterior( + data = getModelData(data, model_name), + prior_list = prior_list) + + crit_pval <- getCritProb( + mods = mods, + dose_levels = dose_levels, + dose_weights = n_patients, + alpha_crit_val = alpha_crit_val) + + contr_mat_prior <- getContrMat( + mods = mods, + dose_levels = dose_levels, + dose_weights = n_patients, + prior_list = prior_list) + + b_mcp_mod <- performBayesianMCPMod( + posteriors_list = posterior_list, + contr_mat = contr_mat_prior, + crit_prob = crit_pval, + simple = simple) + + }) + + names(eval_design) <- model_names + + return (eval_design) + +} + +#' @title getContrMat +#' +#' @param mods tbd +#' @param dose_levels tbd +#' @param dose_weights tbd +#' @param prior_list tbd +#' +#' @export +getContrMat <- function ( + + mods, + dose_levels, + dose_weights, + prior_list + +) { + + ess_prior <- suppressMessages(round(unlist(lapply(prior_list, RBesT::ess)))) + + contr_mat <- DoseFinding::optContr( + models = mods, + doses = dose_levels, + w = dose_weights + ess_prior) + + return (contr_mat) + +} + +#' @title getCritProb +#' +#' @param mods tbd +#' @param dose_levels tbd +#' @param dose_weights tbd +#' @param alpha_crit_val tbd +#' +#' @export +getCritProb <- function ( + + mods, + dose_levels, + dose_weights, + alpha_crit_val + +) { + + contr_mat <- DoseFinding::optContr( + models = mods, + doses = dose_levels, + w = dose_weights) + + crit_pval <- pnorm(DoseFinding:::critVal( + corMat = contr_mat$corMat, + alpha = alpha_crit_val, + df = 0, + alternative = "one.sided")) + + return (crit_pval) + +} + +#' @title performBayesianMCPMod +#' +#' @param posteriors_list tbd +#' @param contr_mat tbd #' @param crit_prob tbd -#' @param n_simulations tbd -#' +#' @param simple tbd +#' #' @export -BMCPMod <- function( - ancova1, - cont_Mat1, - crit_prob, - n_simulations) { - adj1_p <- list() +performBayesianMCPMod <- function ( + + posteriors_list, + contr_mat, + crit_prob, + simple = FALSE + +) { + + if (class(posteriors_list) == "postList") { + + posteriors_list <- list(posteriors_list) + + } + + b_mcp <- performBayesianMCP( + posteriors_list = posteriors_list, + contr_mat = contr_mat, + crit_prob = crit_prob) + + model_shapes <- colnames(contr_mat$contMat) + dose_levels <- as.numeric(rownames(contr_mat$contMat)) + + fits_list <- lapply(seq_along(posteriors_list), function (i) { + + if (b_mcp[i, 1]) { + + sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "crit_prob") + + model_fits <- getModelFits( + models = model_shapes, + dose_levels = dose_levels, + posterior = posteriors_list[[i]], + simple = simple) + + model_fits <- addSignificance(model_fits, sign_models) + + } else { + + NULL + + } + + }) + + bmcpmod <- list(BayesianMCP = b_mcp, Mod = fits_list) + class(bmcpmod) <- "BayesianMCPMod" + + return (bmcpmod) + +} + +addSignificance <- function ( + + model_fits, + sign_models + +) { + + names(sign_models) <- NULL + + model_fits_out <- lapply(seq_along(model_fits), function (i) { + + c(model_fits[[i]], significant = sign_models[i]) + + }) + + attributes(model_fits_out) <- attributes(model_fits) + + return (model_fits_out) + +} - for (i in 1:n_simulations) { - ancova <- ancova1[[i]] - adj1_p[[i]] <- BayesMCPMod( - ancova, - cont_Mat1, - crit_prob - ) +#' @title BayesianMCP +#' +#' @param posteriors_list tbd +#' @param contr_mat tbd +#' @param crit_prob tbd +#' +#' @export +performBayesianMCP <- function( + + posteriors_list, + contr_mat, + crit_prob + +) { + + if (class(posteriors_list) == "postList") { + + posteriors_list <- list(posteriors_list) + } + + b_mcp <- t(sapply(posteriors_list, BayesMCPi, contr_mat, crit_prob)) + + attr(b_mcp, "crit_prob") <- crit_prob + class(b_mcp) <- "BayesianMCP" + + return (b_mcp) + +} - return(adj1_p) +BayesMCPi <- function ( + + posterior_i, + contr_mat, + crit_prob + +) { + + getPostProb <- function ( + + contr_j, # j: dose level + post_combs_i # i: simulation outcome + + ) { + + ## Test statistic = sum over all components of + ## posterior weight * normal probability distribution of + ## critical values for doses * estimated mean / sqrt(product of critical values for doses) + + ## Calculation for each component of the posterior + contr_theta <- apply(post_combs_i$means, 1, `%*%`, contr_j) + contr_var <- apply(post_combs_i$vars, 1, `%*%`, contr_j^2) + contr_weights <- post_combs_i$weights + + ## P(c_m * theta > 0 | Y = y) for a shape m (and dose j) + post_probs <- sum(contr_weights * stats::pnorm(contr_theta / sqrt(contr_var))) + + return (post_probs) + + } + + post_combs_i <- getPostCombsI(posterior_i) + post_probs <- apply(contr_mat$contMat, 2, getPostProb, post_combs_i) + + res <- c(sign = ifelse(max(post_probs) > crit_prob, 1, 0), + p_val = max(post_probs), + post_probs = post_probs) + + return (res) + } diff --git a/R/bootstrapping.R b/R/bootstrapping.R index e69de29..8585b83 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -0,0 +1,92 @@ +#' @title getBootstrapBands +#' +#' @param model_fits tbd +#' @param n_samples tbd +#' @param alpha tbd +#' @param avg_fit tbd +#' @param dose_seq tbd +#' +#' @return tbd +#' @export +getBootstrapBands <- function ( + + model_fits, + n_samples = 1e3, + alpha = c(0.05, 0.5), + avg_fit = TRUE, + dose_seq = NULL + +) { + + mu_hat_samples <- sapply(attr(model_fits, "posterior"), + RBesT::rmix, n = n_samples) + sd_hat <- summary.postList(attr(model_fits, "posterior"))[, 2] + + dose_levels <- model_fits[[1]]$dose_levels + model_names <- names(model_fits) + quantile_probs <- c(0.5, sort(unique(c(alpha / 2, 1 - alpha / 2)))) + + if (is.null(dose_seq)) { + + dose_seq <- seq(min(dose_levels), max(dose_levels), length.out = 100L) + + } + + preds <- apply(mu_hat_samples, 1, function (mu_hat) { + + preds_mu_hat <- sapply(model_names, function (model) { + + fit <- DoseFinding::fitMod( + dose = model_fits[[1]]$dose_levels, + resp = mu_hat, + S = diag(sd_hat^2), + model = model, + type = "general", + bnds = DoseFinding::defBnds( + mD = max(model_fits[[1]]$dose_levels))[[model]]) + + preds <- stats::predict(fit, doseSeq = dose_seq, predType = "ls-means") + attr(preds, "gAIC") <- DoseFinding::gAIC(fit) + + return (preds) + + }, simplify = FALSE) + + preds_mu_mat <- do.call(rbind, preds_mu_hat) + + if (avg_fit) { + + avg_fit_indx <- which.min(sapply(preds_mu_hat, attr, "gAIC")) + preds_mu_mat <- rbind(preds_mu_mat, avgFit = preds_mu_mat[avg_fit_indx, ]) + + } + + return (preds_mu_mat) + + }) + + if (avg_fit) { + + model_names <- c(model_names, "avgFit") + + } + + sort_indx <- order(rep(seq_along(model_names), length(dose_seq))) + quant_mat <- t(apply(X = preds[sort_indx, ], + MARGIN = 1, + FUN = stats::quantile, + probs = quantile_probs)) + + cr_bounds_data <- cbind( + dose_seqs = dose_seq, + models = rep( + x = factor(model_names, + levels = c("linear", "emax", "exponential", + "sigEmax", "logistic", "quadratic", + "avgFit")), + each = length(dose_seq)), + as.data.frame(quant_mat)) + + return (cr_bounds_data) + +} diff --git a/R/dataSimulation.R b/R/dataSimulation.R deleted file mode 100644 index f601578..0000000 --- a/R/dataSimulation.R +++ /dev/null @@ -1,43 +0,0 @@ -#' @title simulateData -#' @param n_patients tbd -#' -#' @param dose_levels tbd -#' @param pcov tbd -#' @param SD tbd -#' @param placebo_effect tbd -#' @param max_effect tbd -#' @param mod tbd -#' @param n_simulations tbd -#' -#' @export -simulateData <- function( - n_patients, - dose_levels, - pcov, - SD, - placebo_effect, - max_effect, - mod, - n_simulations) { - ntrt <- length(dose_levels) - sum_patients <- sum(n_patients) - - simulno <- rep(1:n_simulations, - each = sum_patients - ) - ptno <- rep(1:sum_patients, - times = n_simulations - ) - dose <- rep(rep(dose_levels, times = n_patients), - times = n_simulations - ) - - df <- data.frame(simulation = simulno, ptno, dose) - df <- cbind( - df, - DoseFinding::getResp(mod, dose) + - rep(stats::rnorm(length(dose), mean = 0, sd = SD), times = length(mod)) - ) - - return(df) -} diff --git a/R/helper.R b/R/helper.R deleted file mode 100644 index d5aeab1..0000000 --- a/R/helper.R +++ /dev/null @@ -1,116 +0,0 @@ -#' @title BayesMCPMod -#' -#' @param ancova tbd -#' @param contMat1 tbd -#' @param crit_prob tbd -#' -#' @return tbd - -BayesMCPMod <- function( - ancova, - contMat1, - crit_prob) { - out <- list() - contMat <- contMat1[[1]] - - m <- length(contMat[1, ]) - probs <- rep(0, m) - adjpval <- rep(0, m) - sign <- 0 - - for (j in 1:m) { - conttheta <- sapply(seq_len(nrow(ancova[["est."]])), function(x) t(contMat[, j]) %*% ancova[["est."]][x, ]) - contvar <- sapply(seq_len(nrow(ancova[["est."]])), function(x) t(contMat[, j]^2) %*% ancova[["cov."]][x, ]^2) - p_ij <- sapply(seq_len(nrow(ancova[["est."]])), function(x) stats::pnorm(conttheta[x] / sqrt(contvar[x]))) - - - probs[j] <- sum(ancova[["obs."]] * p_ij) - } - - max_prob <- max(probs) - - if (max_prob > crit_prob) { - sign <- 1 - } - - out[["sign."]] <- sign - out[["adjpval."]] <- max_prob - out[["posteriorProb"]] <- probs - - return(out) -} - -#' @title doAncova -#' -#' @param dose tbd -#' @param response tbd -#' @param prior tbd -#' -#' @return tbd - -doAncova <- function( - dose, - response, - prior) { - mixnorm <- RBesT::mixnorm - - - anova <- stats::lm(response ~ factor(dose) - 1) - mean.sim <- anova$coefficients - se.mean.sim <- summary(anova)$coefficients[, 2] - - n.dose <- unique(prior[, 1]) - - doseAct <- n.dose[-1] - - k <- length(n.dose) - - prior_dose <- lapply(1:k, function(x) matrix(subset(prior, prior[, 1] == n.dose[x])[, -1], ncol = 3)) - - comb.prior <- sapply(1:k, function(x) nrow(prior_dose[[x]])) - - - post_sep <- list() - - for (i in 1:k) { - args <- lapply(1:comb.prior[i], function(x) prior_dose[[i]][x, ]) - mixed.prior <- do.call("mixnorm", args) - post_sep[[i]] <- RBesT::postmix(mixed.prior, m = mean.sim[i], se = se.mean.sim[i]) - } - - comb.post <- prod(comb.prior) - - args <- lapply(1:k, function(x) 1:comb.prior[x]) - comb.ind <- do.call("expand.grid", args) - - post.weight <- matrix(sapply(1:k, function(x) sapply(1:comb.post, function(y) post_sep[[x]][1, comb.ind[y, x]])), nrow = comb.post) - post.weight <- apply(post.weight, 1, prod) - post.mean <- matrix(sapply(1:k, function(x) sapply(1:comb.post, function(y) post_sep[[x]][2, comb.ind[y, x]])), nrow = comb.post) - post.sd <- matrix(sapply(1:k, function(x) sapply(1:comb.post, function(y) post_sep[[x]][3, comb.ind[y, x]])), nrow = comb.post) - post.sd.overall <- sapply(1:k, function(x) summary(post_sep[[x]])[2]) - post.mean.pbo <- RBesT::qmix(post_sep[[1]], p = 0.5) - post.mean.act1 <- RBesT::qmix(post_sep[[2]], p = 0.5) - post.mean.act2 <- RBesT::qmix(post_sep[[3]], p = 0.5) - post.mean.act3 <- RBesT::qmix(post_sep[[4]], p = 0.5) - diff <- c( - post.mean.act1 - post.mean.pbo, - post.mean.act2 - post.mean.pbo, - post.mean.act3 - post.mean.pbo - ) - post.mean.max.diff <- c(max(diff), which.max(diff)) - - output <- list() - - - output[["obs."]] <- post.weight - output[["est."]] <- post.mean - output[["cov."]] <- post.sd - output[["post.sd."]] <- post.sd.overall - output[["pbomedian."]] <- post.mean.pbo - output[["act1median."]] <- post.mean.act1 - output[["act2median."]] <- post.mean.act2 - output[["act3median."]] <- post.mean.act3 - output[["maxdiffmedian+which."]] <- post.mean.max.diff - - return(output) -} \ No newline at end of file diff --git a/R/modelling.R b/R/modelling.R index e6abc27..7895095 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -1,238 +1,259 @@ -#' @title doFit +#' @title getModelFits +#' +#' @param models tbd #' @param dose_levels tbd -#' -#' @param shapes tbd #' @param posterior tbd -#' +#' @param simple tbd +#' #' @export -doFit <- function( - dose_levels, - shapes, - posterior) { - mD <- max(dose_levels) - bounds <- DoseFinding::defBnds( - mD = mD, - emax = c(0.001, 1.5) * mD, - exponential = c(0.1, 2) * mD, - logistic = matrix(c(0.001, 0.01, 1.5, 1 / 2) * mD, 2) - ) - - summary <- list() - - resp <- c( - posterior$pbomedian, - posterior$act1median, - posterior$act2median, - posterior$act3median - ) - - cov <- diag(posterior$post.sd) - - for (j in seq_along(shapes)) { - sublist <- list() - shape <- shapes[j] - fit <- DoseFinding::fitMod( - dose = dose_levels, - resp = resp, - S = cov, - model = shape, - type = "general", - bnds = bounds[[shape]] - ) - fit_vec <- stats::predict(fit, predType = "ls-means") - gAIC <- DoseFinding::gAIC(fit) - max_effect <- max(fit_vec) - min(fit_vec) +getModelFits <- function ( + + models, + dose_levels, + posterior, + simple = FALSE + +) { + + getModelFit <- ifelse(simple, getModelFitSimple, getModelFitOpt) + model_fits <- lapply(models, getModelFit, dose_levels, posterior) + + model_fits <- addModelWeights(model_fits) + + names(model_fits) <- models + attr(model_fits, "posterior") <- posterior + class(model_fits) <- "modelFits" + + return (model_fits) + +} - sublist[["model"]] <- shape - sublist[["fit"]] <- fit$coefs - sublist[["pred"]] <- fit_vec - sublist[["gAIC"]] <- gAIC - sublist[["max_effect"]] <- max_effect +## ignoring mixture posterior +getModelFitSimple <- function ( + + model, + dose_levels, + posterior + +) { + + fit <- DoseFinding::fitMod( + dose = dose_levels, + resp = summary.postList(posterior)[, 1], + S = diag(summary.postList(posterior)[, 2]^2), + 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, + dose_levels = dose_levels, + pred_values = pred_vals, + max_effect = max_effect, + gAIC = gAIC) + + return (model_fit) + +} - summary[[j]] <- sublist +## taking mixture posterior into account +getModelFitOpt <- function ( + + model, + dose_levels, + posterior + +) { + + getOptParams <- function ( + + model, + dose_levels, + posterior + + ) { + + 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 (GENERAL$ERROR$MODEL_OPTIONS)}) + + 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), + params = list( + x0 = simple_fit$coeffs, + lb = lb, + ub = ub), + c_names = names(simple_fit$coeffs)) + + return (param_list) + } - - return(summary) + + 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, + lb = param_list$params$lb, + ub = param_list$params$ub, + opts = param_list$opts, + dose_levels = dose_levels, + post_combs = post_combs, + 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$max_effect <- max(model_fit$pred_values) - min(model_fit$pred_values) + model_fit$gAIC <- getGenAIC(model_fit, post_combs) + + return (model_fit) + } -#' @title estimateModels -#' @param models tbd -#' -#' @param post tbd -#' @param doses tbd -#' -#' @export -estimateModels <- function( - models, - post, - doses) { - fitModel <- function(model, post, doses) { - getF <- function( - theta, - post, - dose, - expression_i) { - L <- length(post$obs) - comp <- rep(0, L) - - for (i in 1:L) { - comp[i] <- eval(expression_i) - } - - return(sum(post$obs * comp)) - } - - switch(model, - "emax" = { - x0 <- c(0, 1, 0.5) - maxeval <- 5000 - lb <- c(-Inf, -Inf, 0.001 * max(doses)) - ub <- c(Inf, Inf, 1.5 * max(doses)) - expression_i <- quote(sum((post$est[i, ] - (theta[1] + (theta[2] * dose) / (theta[3] + dose)))^2 / (post$cov[i, ]^2))) - }, - "sigEmax" = { - x0 <- c(0, 1, 0.5, 1) - maxeval <- 10000 - lb <- c(-Inf, -Inf, 0.001 * max(doses), 0.5) - ub <- c(Inf, Inf, 1.5 * max(doses), 10) - expression_i <- quote(sum((post$Mean[i, ] - (theta[1] + (theta[2] * dose^theta[4]) / (theta[3]^theta[4] + dose^theta[4])))^2 / (post$cov[i, ]^2))) - }, - "exponential" = { - x0 <- c(-1, 1, 20) - maxeval <- 10000 - lb <- c(-Inf, -Inf, 0.1 * max(doses)) - ub <- c(Inf, Inf, 2 * max(doses)) - expression_i <- quote(sum((post$est[i, ] - (theta[1] + theta[2] * (exp(dose / theta[3]) - 1)))^2 / (post$cov[i, ]^2))) - }, - "quadratic" = { - x0 <- c(0, 1, -1) - maxeval <- 10000 - lb <- NULL - ub <- c(Inf, Inf, 0) - expression_i <- quote(sum((post$est[i, ] - (theta[1] + theta[2] * dose + theta[3] * dose^2))^2 / (post$cov[i, ]^2))) - }, - "linear" = { - x0 <- c(0, 1) - maxeval <- 10000 - lb <- NULL - ub <- NULL - expression_i <- quote(sum((post$est[i, ] - (theta[1] + theta[2] * dose))^2 / (post$cov[i, ]^2))) - }, - "logistic" = { - x0 <- c(0, 1, 0.5, 5) - maxeval <- 10000 - lb <- c(-Inf, -Inf, 0.001 * max(doses), 0.01 * max(doses)) - ub <- c(Inf, Inf, 1.5 * max(doses), 0.5 * max(doses)) - expression_i <- quote(sum((post$est[i, ] - (theta[1] + theta[2] / (1 + exp((theta[3] - dose) / theta[4]))))^2 / (post$cov[i, ]^2))) - }, - { - stop("model must be one of 'emax', 'sigEmax', 'exponential', 'linear', 'logistic', 'quadratic'") - } - ) - res <- nloptr::nloptr( - x0 = x0, - eval_f = function(theta) { - getF(theta, - post = post, - dose = doses, - expression_i = expression_i - ) - }, - opts = list( - "algorithm" = "NLOPT_LN_NELDERMEAD", - maxeval = maxeval - ), - lb = lb, - ub = ub - ) - - return(list( - model = model, - fit = res - )) +predictModelFit <- function ( + + model_fit, + doses = NULL + +) { + + if (is.null(doses)) { + + doses <- model_fit$dose_levels + } - - results <- lapply(models, function(model) { - fitModel(model, post, doses) - }) - - return(results) + + pred_vals <- switch ( + model_fit$model, + "emax" = {DoseFinding::emax( + doses, + model_fit$coeffs["e0"], + model_fit$coeffs["eMax"], + model_fit$coeffs["ed50"])}, + "sigEmax" = {DoseFinding::sigEmax( + doses, + model_fit$coeffs["e0"], + model_fit$coeffs["eMax"], + model_fit$coeffs["ed50"], + model_fit$coeffs["h"])}, + "exponential" = {DoseFinding::exponential( + doses, + model_fit$coeffs["e0"], + model_fit$coeffs["e1"], + model_fit$coeffs["delta"])}, + "quadratic" = {DoseFinding::quadratic( + doses, + model_fit$coeffs["e0"], + model_fit$coeffs["b1"], + model_fit$coeffs["b2"])}, + "linear" = {DoseFinding::linear( + doses, + model_fit$coeffs["e0"], + model_fit$coeffs["delta"])}, + "logistic" = {DoseFinding::logistic( + doses, + model_fit$coeffs["e0"], + model_fit$coeffs["eMax"], + model_fit$coeffs["ed50"], + model_fit$coeffs["delta"])}, + {stop(GENERAL$ERROR$MODEL_OPTIONS)}) + + return (pred_vals) + } -#' @title getGenAICs -#' @param post tbd -#' -#' @param fit_mods tbd -#' @param doses tbd -#' -#' @export -getGenAICs <- function( - post, - fit_mods, - doses) { - getGenAIC <- function( - post, - fit_mod, - doses) { - model <- fit_mod$model - p <- length(fit_mod$fit$solution) - theta <- fit_mod$fit$solution - L <- length(post$obs) - - fitted <- switch(model, - "emax" = { - DoseFinding::emax(doses, theta[1], theta[2], theta[3]) - }, - "sigEmax" = { - DoseFinding::sigEmax(doses, theta[1], theta[2], theta[3], theta[4]) - }, - "exponential" = { - DoseFinding::exponential(doses, theta[1], theta[2], theta[3]) - }, - "quadratic" = { - DoseFinding::quadratic(doses, theta[1], theta[2], theta[3]) - }, - "linear" = { - DoseFinding::linear(doses, theta[1], theta[2]) - }, - "logistic" = { - DoseFinding::logistic(doses, theta[1], theta[2], theta[3], theta[4]) - }, - { - stop("model must be one of 'emax', 'sigEmax', 'exponential', 'linear', 'logistic', 'quadratic'") - } - ) - - gAIC <- rep(0, L) - - for (i in 1:L) { - gAIC[i] <- sum(1 / post$cov[i, ]^2 * (post$est[i, ] - fitted)^2) + 2 * p - } - - avgAIC <- stats::weighted.mean(gAIC, w = post$obs) - - res <- list( - model = model, - fitted = fitted, - gAIC = gAIC, - avgAIC = avgAIC, - weightedAIC = NULL - ) - - return(res) - } +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) + +} - results <- lapply(fit_mods, function(fit_mod) { - getGenAIC(post, fit_mod, doses) +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]) + }) - - avgAIC_values <- sapply(results, function(res) res$avgAIC) - exp_values <- exp(-0.5 * avgAIC_values) - denominator <- sum(exp_values) - - for (i in seq_along(results)) { - results[[i]]$weightedAIC <- exp_values[i] / denominator - } - - return(results) + + return (model_fits_out) + } diff --git a/R/plot.R b/R/plot.R index e69de29..9f8278b 100644 --- a/R/plot.R +++ b/R/plot.R @@ -0,0 +1,184 @@ +#' @title plot.modelFits +#' +#' @param x tbd An object of type modelFits +#' @param gAIC tbd +#' @param avg_fit tbd +#' @param cr_intv tbd +#' @param alpha_CrI tbd +#' @param cr_bands tbd +#' @param alpha_CrB tbd +#' @param n_bs_smpl tbd +#' @param acc_color tbd +#' @param ... tbd +#' +#' @return tbd +#' @export +plot.modelFits <- function ( + + x, + gAIC = TRUE, + avg_fit = TRUE, + cr_intv = TRUE, + alpha_CrI = 0.05, + cr_bands = FALSE, + alpha_CrB = c(0.05, 0.5), + n_bs_smpl = 1e3, + acc_color = "orange", + ... + +) { + + plot_res <- 1e2 + model_fits <- x + + dose_levels <- model_fits[[1]]$dose_levels + 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) + + preds_models <- sapply(model_fits, predictModelFit, doses = doses) + model_names <- names(model_fits) + + if (avg_fit) { + + mod_weights <- sapply(model_fits, function (x) x$model_weight) + avg_mod <- preds_models %*% mod_weights + + preds_models <- cbind(preds_models, avg_mod) + model_names <- c(model_names, "avgFit") + + } + + gg_data <- data.frame( + dose_seqs = rep(doses, 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)) + + 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 = _) + + label_avg <- paste0(paste_names, "=", round(mod_weights, 1), + collapse = ", ") + label_gAUC <- c(label_gAUC, label_avg) + + } + + } + + if (cr_bands) { + + crB_data <- getBootstrapBands( + model_fits = model_fits, + n_samples = n_bs_smpl, + alpha = alpha_CrB, + avg_fit = avg_fit, + dose_seq = doses) + + getInx <- function (alpha_CrB) { + n <- length(alpha_CrB) + inx_list <- lapply(seq_len(n), function (i) c(i, 2 * n - i + 1) + 3) + return (inx_list)} + + } + + plts <- ggplot2::ggplot() + + ## Layout etc. + 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), + label = label_gAUC), + 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 = dose_seqs, + ymin = crB_data[, ", inx[1], "], + ymax = crB_data[, ", inx[2], "]), + 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(dose_seqs, `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], + ymax = post_summary[, 5]), + mapping = ggplot2::aes(x = x, + ymin = ymin, + ymax = ymax), + width = 0, + alpha = 0.5) + + } + } + + ## Posterior Medians + ggplot2::geom_point( + data = data.frame( + dose_levels = dose_levels, + medians = post_summary[, 4]), + mapping = ggplot2::aes(dose_levels, medians), + size = 2) + + ## Fitted Models + ggplot2::geom_line( + data = gg_data, + mapping = ggplot2::aes(dose_seqs, fits)) + + ## Faceting + ggplot2::facet_wrap(~ models) + + return (plts) + +} \ No newline at end of file diff --git a/R/postShape.R b/R/postShape.R deleted file mode 100644 index 8da36a6..0000000 --- a/R/postShape.R +++ /dev/null @@ -1,26 +0,0 @@ -#' @title postShape -#' @param data tbd -#' -#' @param prior tbd -#' @param n_simulations tbd -#' -#' @export -postShape <- function( - data, - prior, - n_simulations) { - ancova <- list() - for (i in seq_len(n_simulations)) { - datai <- data[data$simulation == i, ] - dosei <- datai$dose - responsei <- datai$response - ancovai <- doAncova( - dose = dosei, - response = responsei, - prior = prior - ) - ancova[[i]] <- ancovai - } - - return(ancova) -} diff --git a/R/posterior.R b/R/posterior.R new file mode 100644 index 0000000..df37f5e --- /dev/null +++ b/R/posterior.R @@ -0,0 +1,156 @@ +#' @title getPriorList +#' +#' @param hist_data tbd +#' @param dose_levels tbd +#' @param dose_names prior_list +#' @param robustify_weight tbd +#' +#' @export +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 +#' @param prior_list prior_list +#' @param mu_hat tbd +#' @param sd_hat tbd +#' +#' @export +getPosterior <- function( + data, + prior_list, + mu_hat = NULL, + sd_hat = NULL + +) { + + posterior_list <- lapply(split(data, data$simulation), getPosteriorI, + prior_list = prior_list, + mu_hat = mu_hat, + sd_hat = sd_hat) + + if (length(posterior_list) == 1) { + + posterior_list <- posterior_list[[1]] + + } + + return (posterior_list) + +} + +getPosteriorI <- function( + + data_i, + prior_list, + mu_hat = NULL, + sd_hat = NULL + +) { + + if (is.null(mu_hat) && is.null(sd_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] + + } else if (!is.null(mu_hat) && !is.null(sd_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)) + + } else { + + stop ("Both mu_hat and sd_hat must be provided.") + + } + + post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = sd_hat) + + if (is.null(names(prior_list))) { + + names(prior_list) <- c("Ctr", paste0("DG_", seq_along(post_list[-1]))) + + } + + names(post_list) <- names(prior_list) + class(post_list) <- "postList" + + return (post_list) + +} + +getPostCombsI <- function ( + + posterior_i + +) { + + post_params <- list( + weights = lapply(posterior_i, function (x) x[1, ]), + means = lapply(posterior_i, function (x) x[2, ]), + vars = lapply(posterior_i, function (x) x[3, ]^2)) + + post_combs <- lapply(post_params, expand.grid) + post_combs$weights <- apply(post_combs$weights, 1, prod) + + return (post_combs) + +} + + diff --git a/R/s3methods.R b/R/s3methods.R new file mode 100644 index 0000000..14aa663 --- /dev/null +++ b/R/s3methods.R @@ -0,0 +1,190 @@ +## BayesianMCPMod ----------------------------------------- + +#' @export +print.BayesianMCPMod <- function ( + + x, + ... + +) { + + n_models <- ncol(x$BayesianMCP) - 2L + model_names <- colnames(x$BayesianMCP)[-c(1, 2)] |> + sub(pattern = "post_probs.", replacement = "", x = _) + + model_success <- colMeans(do.call(rbind, lapply(x$Mod, function (y) { + + if (!is.null(y)) { + + sapply(y, function (z) z$significant) + + } else { + + model_signs <- rep(FALSE, n_models) + names(model_signs) <- model_names + + return (model_signs) + + } + + }))) + + print(x$BayesianMCP) + cat("\n") + cat("Model Significance Frequencies\n") + print(model_success, ...) + +} + +## BayesianMCP -------------------------------------------- + +#' @export +print.BayesianMCP <- function ( + + x, + ... + +) { + + power <- mean(x[, 1]) + n_sim <- nrow(x) + + cat("Bayesian Multiple Comparison Procedure\n") + cat(" Estimated Success Rate: ", power, "\n") + cat(" N Simulations: ", n_sim) + +} + +## ModelFits ---------------------------------------------- + +#' @export +predict.ModelFits <- function ( + + object, + doses = NULL, + ... + +) { + + lapply(object, predictModelFit, doses = doses, ...) + +} + +#' @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)) + 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 = _) |> + gsub("linear", "linear ", x = _) |> + 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", + paste(dose_names, round(dose_levels, n_digits), sep = " = "), "\n") + cat("\n") + cat("Predictions, Maximum Effect, gAIC & Significance\n") + print(out_table, ...) + +} + +## plot.ModelFits() +## see file R/plot.R + +## postList ----------------------------------------------- + +#' @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/R/simulation.R b/R/simulation.R new file mode 100644 index 0000000..7c60082 --- /dev/null +++ b/R/simulation.R @@ -0,0 +1,64 @@ +#' @title simulateData +#' +#' @param n_patients tbd +#' @param dose_levels tbd +#' @param sd tbd +#' @param mods tbd +#' @param n_sim tbd +#' @param true_model tbd +#' +#' @export +simulateData <- function( + + n_patients, + dose_levels, + sd, + mods, + n_sim = 1e3, + true_model = NULL + +) { + + if (!is.null(true_model)) { + + n_sim <- 1 + mods_attr <- attributes(mods) + mods_attr$names <- true_model + mods <- mods[true_model] + attributes(mods) <- mods_attr + + } + + sim_info <- data.frame( + simulation = rep(seq_len(n_sim), each = sum(n_patients)), + 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) + + sim_data <- cbind(sim_info, model_responses + random_noise) + + if (!is.null(true_model)) { + + sim_data <- getModelData(sim_data, true_model) + + } + + return (sim_data) + +} + +getModelData <- function ( + + sim_data, + model_name + +) { + + model_data <- sim_data[, c("simulation", "dose", model_name)] + colnames(model_data)[3] <- "response" + + return (model_data) + +} \ No newline at end of file diff --git a/README.md b/README.md index 66ebd6e..b4cc3fb 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,10 @@ + +[![Linux](https://github.com/Boehringer-Ingelheim/BayesianMCPMod/actions/workflows/linux.yml/badge.svg?branch=main)](https://github.com/Boehringer-Ingelheim/BayesianMCPMod/actions/workflows/linux.yml) +[![Windows](https://github.com/Boehringer-Ingelheim/BayesianMCPMod/actions/workflows/windows.yml/badge.svg?branch=main)](https://github.com/Boehringer-Ingelheim/BayesianMCPMod/actions/workflows/windows.yml) +[![Macos](https://github.com/Boehringer-Ingelheim/BayesianMCPMod/actions/workflows/macos.yml/badge.svg?branch=main)](https://github.com/Boehringer-Ingelheim/BayesianMCPMod/actions/workflows/macos.yml) +[![Codecov test coverage](https://codecov.io/gh/Boehringer-Ingelheim/BayesianMCPMod/branch/main/graph/badge.svg)](https://codecov.io/gh/Boehringer-Ingelheim/BayesianMCPMod?branch=main) + + # The `{BayesianMCPmod}` package ... diff --git a/_pkgdown.yml b/_pkgdown.yml new file mode 100644 index 0000000..1b0d774 --- /dev/null +++ b/_pkgdown.yml @@ -0,0 +1,15 @@ +url: https://boehringer-ingelheim.github.io/BayesianMCPMod/ +template: + bootstrap: 5 + +navbar: + structure: + left: [intro, articles, reference, news] + right: [search, github] + +articles: +- title: Get started + navbar: ~ + contents: + - analysis_normal + - Simulation_Example diff --git a/man/BMCPMod.Rd b/man/BMCPMod.Rd deleted file mode 100644 index 1d82663..0000000 --- a/man/BMCPMod.Rd +++ /dev/null @@ -1,20 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/BMCPMod.R -\name{BMCPMod} -\alias{BMCPMod} -\title{BMCPMod} -\usage{ -BMCPMod(ancova1, cont_Mat1, crit_prob, n_simulations) -} -\arguments{ -\item{ancova1}{tbd} - -\item{cont_Mat1}{tbd} - -\item{crit_prob}{tbd} - -\item{n_simulations}{tbd} -} -\description{ -BMCPMod -} diff --git a/man/BayesMCPMod.Rd b/man/BayesMCPMod.Rd deleted file mode 100644 index d6df56f..0000000 --- a/man/BayesMCPMod.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/helper.R -\name{BayesMCPMod} -\alias{BayesMCPMod} -\title{BayesMCPMod} -\usage{ -BayesMCPMod(ancova, contMat1, crit_prob) -} -\arguments{ -\item{ancova}{tbd} - -\item{contMat1}{tbd} - -\item{crit_prob}{tbd} -} -\value{ -tbd -} -\description{ -BayesMCPMod -} diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd new file mode 100644 index 0000000..6883d23 --- /dev/null +++ b/man/assessDesign.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/BMCPMod.R +\name{assessDesign} +\alias{assessDesign} +\title{assessDesign} +\usage{ +assessDesign( + n_patients, + mods, + prior_list, + n_sim = 1000, + alpha_crit_val = 0.05, + simple = TRUE +) +} +\arguments{ +\item{n_patients}{tbd} + +\item{mods}{tbd} + +\item{prior_list}{tbd} + +\item{n_sim}{tbd} + +\item{alpha_crit_val}{tbd} + +\item{simple}{tbd} +} +\description{ +assessDesign +} diff --git a/man/doAncova.Rd b/man/doAncova.Rd deleted file mode 100644 index 7540b9a..0000000 --- a/man/doAncova.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/helper.R -\name{doAncova} -\alias{doAncova} -\title{doAncova} -\usage{ -doAncova(dose, response, prior) -} -\arguments{ -\item{dose}{tbd} - -\item{response}{tbd} - -\item{prior}{tbd} -} -\value{ -tbd -} -\description{ -doAncova -} diff --git a/man/estimateModels.Rd b/man/estimateModels.Rd deleted file mode 100644 index 59e0200..0000000 --- a/man/estimateModels.Rd +++ /dev/null @@ -1,18 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/modelling.R -\name{estimateModels} -\alias{estimateModels} -\title{estimateModels} -\usage{ -estimateModels(models, post, doses) -} -\arguments{ -\item{models}{tbd} - -\item{post}{tbd} - -\item{doses}{tbd} -} -\description{ -estimateModels -} diff --git a/man/getBootstrapBands.Rd b/man/getBootstrapBands.Rd new file mode 100644 index 0000000..de4bd80 --- /dev/null +++ b/man/getBootstrapBands.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bootstrapping.R +\name{getBootstrapBands} +\alias{getBootstrapBands} +\title{getBootstrapBands} +\usage{ +getBootstrapBands( + model_fits, + n_samples = 1000, + alpha = c(0.05, 0.5), + avg_fit = TRUE, + dose_seq = NULL +) +} +\arguments{ +\item{model_fits}{tbd} + +\item{n_samples}{tbd} + +\item{alpha}{tbd} + +\item{avg_fit}{tbd} + +\item{dose_seq}{tbd} +} +\value{ +tbd +} +\description{ +getBootstrapBands +} diff --git a/man/getContrMat.Rd b/man/getContrMat.Rd new file mode 100644 index 0000000..d5393e3 --- /dev/null +++ b/man/getContrMat.Rd @@ -0,0 +1,20 @@ +% 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/getCritProb.Rd b/man/getCritProb.Rd new file mode 100644 index 0000000..eeb1b2c --- /dev/null +++ b/man/getCritProb.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/BMCPMod.R +\name{getCritProb} +\alias{getCritProb} +\title{getCritProb} +\usage{ +getCritProb(mods, dose_levels, dose_weights, alpha_crit_val) +} +\arguments{ +\item{mods}{tbd} + +\item{dose_levels}{tbd} + +\item{dose_weights}{tbd} + +\item{alpha_crit_val}{tbd} +} +\description{ +getCritProb +} diff --git a/man/getGenAICs.Rd b/man/getGenAICs.Rd deleted file mode 100644 index a81fa2a..0000000 --- a/man/getGenAICs.Rd +++ /dev/null @@ -1,18 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/modelling.R -\name{getGenAICs} -\alias{getGenAICs} -\title{getGenAICs} -\usage{ -getGenAICs(post, fit_mods, doses) -} -\arguments{ -\item{post}{tbd} - -\item{fit_mods}{tbd} - -\item{doses}{tbd} -} -\description{ -getGenAICs -} diff --git a/man/doFit.Rd b/man/getModelFits.Rd similarity index 50% rename from man/doFit.Rd rename to man/getModelFits.Rd index 90395c6..f186384 100644 --- a/man/doFit.Rd +++ b/man/getModelFits.Rd @@ -1,18 +1,20 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/modelling.R -\name{doFit} -\alias{doFit} -\title{doFit} +\name{getModelFits} +\alias{getModelFits} +\title{getModelFits} \usage{ -doFit(dose_levels, shapes, posterior) +getModelFits(models, dose_levels, posterior, simple = FALSE) } \arguments{ -\item{dose_levels}{tbd} +\item{models}{tbd} -\item{shapes}{tbd} +\item{dose_levels}{tbd} \item{posterior}{tbd} + +\item{simple}{tbd} } \description{ -doFit +getModelFits } diff --git a/man/getPosterior.Rd b/man/getPosterior.Rd new file mode 100644 index 0000000..501f5de --- /dev/null +++ b/man/getPosterior.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/posterior.R +\name{getPosterior} +\alias{getPosterior} +\title{getPosterior} +\usage{ +getPosterior(data, prior_list, mu_hat = NULL, sd_hat = NULL) +} +\arguments{ +\item{data}{tbd} + +\item{prior_list}{prior_list} + +\item{mu_hat}{tbd} + +\item{sd_hat}{tbd} +} +\description{ +getPosterior +} diff --git a/man/getPriorList.Rd b/man/getPriorList.Rd new file mode 100644 index 0000000..bb7e346 --- /dev/null +++ b/man/getPriorList.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/posterior.R +\name{getPriorList} +\alias{getPriorList} +\title{getPriorList} +\usage{ +getPriorList(hist_data, dose_levels, dose_names = NULL, robustify_weight = 0.5) +} +\arguments{ +\item{hist_data}{tbd} + +\item{dose_levels}{tbd} + +\item{dose_names}{prior_list} + +\item{robustify_weight}{tbd} +} +\description{ +getPriorList +} diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd new file mode 100644 index 0000000..e6b3543 --- /dev/null +++ b/man/performBayesianMCP.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/BMCPMod.R +\name{performBayesianMCP} +\alias{performBayesianMCP} +\title{BayesianMCP} +\usage{ +performBayesianMCP(posteriors_list, contr_mat, crit_prob) +} +\arguments{ +\item{posteriors_list}{tbd} + +\item{contr_mat}{tbd} + +\item{crit_prob}{tbd} +} +\description{ +BayesianMCP +} diff --git a/man/performBayesianMCPMod.Rd b/man/performBayesianMCPMod.Rd new file mode 100644 index 0000000..83a4a06 --- /dev/null +++ b/man/performBayesianMCPMod.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/BMCPMod.R +\name{performBayesianMCPMod} +\alias{performBayesianMCPMod} +\title{performBayesianMCPMod} +\usage{ +performBayesianMCPMod(posteriors_list, contr_mat, crit_prob, simple = FALSE) +} +\arguments{ +\item{posteriors_list}{tbd} + +\item{contr_mat}{tbd} + +\item{crit_prob}{tbd} + +\item{simple}{tbd} +} +\description{ +performBayesianMCPMod +} diff --git a/man/plot.modelFits.Rd b/man/plot.modelFits.Rd new file mode 100644 index 0000000..3b35078 --- /dev/null +++ b/man/plot.modelFits.Rd @@ -0,0 +1,46 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot.R +\name{plot.modelFits} +\alias{plot.modelFits} +\title{plot.modelFits} +\usage{ +\method{plot}{modelFits}( + x, + gAIC = TRUE, + avg_fit = TRUE, + cr_intv = TRUE, + alpha_CrI = 0.05, + cr_bands = FALSE, + alpha_CrB = c(0.05, 0.5), + n_bs_smpl = 1000, + acc_color = "orange", + ... +) +} +\arguments{ +\item{x}{tbd An object of type modelFits} + +\item{gAIC}{tbd} + +\item{avg_fit}{tbd} + +\item{cr_intv}{tbd} + +\item{alpha_CrI}{tbd} + +\item{cr_bands}{tbd} + +\item{alpha_CrB}{tbd} + +\item{n_bs_smpl}{tbd} + +\item{acc_color}{tbd} + +\item{...}{tbd} +} +\value{ +tbd +} +\description{ +plot.modelFits +} diff --git a/man/postShape.Rd b/man/postShape.Rd deleted file mode 100644 index da9a36f..0000000 --- a/man/postShape.Rd +++ /dev/null @@ -1,18 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/postShape.R -\name{postShape} -\alias{postShape} -\title{postShape} -\usage{ -postShape(data, prior, n_simulations) -} -\arguments{ -\item{data}{tbd} - -\item{prior}{tbd} - -\item{n_simulations}{tbd} -} -\description{ -postShape -} diff --git a/man/simulateData.Rd b/man/simulateData.Rd index 606e557..32038df 100644 --- a/man/simulateData.Rd +++ b/man/simulateData.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/dataSimulation.R +% Please edit documentation in R/simulation.R \name{simulateData} \alias{simulateData} \title{simulateData} @@ -7,12 +7,10 @@ simulateData( n_patients, dose_levels, - pcov, - SD, - placebo_effect, - max_effect, - mod, - n_simulations + sd, + mods, + n_sim = 1000, + true_model = NULL ) } \arguments{ @@ -20,17 +18,13 @@ simulateData( \item{dose_levels}{tbd} -\item{pcov}{tbd} +\item{sd}{tbd} -\item{SD}{tbd} +\item{mods}{tbd} -\item{placebo_effect}{tbd} +\item{n_sim}{tbd} -\item{max_effect}{tbd} - -\item{mod}{tbd} - -\item{n_simulations}{tbd} +\item{true_model}{tbd} } \description{ simulateData diff --git a/tests/testthat/test-bootstrapping.R b/tests/testthat/test-bootstrapping.R new file mode 100644 index 0000000..c4caf87 --- /dev/null +++ b/tests/testthat/test-bootstrapping.R @@ -0,0 +1,152 @@ +# 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-global.R b/tests/testthat/test-global.R index 6e40d42..ab74964 100644 --- a/tests/testthat/test-global.R +++ b/tests/testthat/test-global.R @@ -1,14 +1,3 @@ test_that("multiplication works", { expect_equal(2 * 2, 4) -}) - -test_that("dummy tests", { - expect_error(calc_pow()) - expect_error(simulateData()) - expect_error(postShape()) - expect_error(BayesMCPMod()) - expect_error(estimateModel()) - expect_error(BMCPMod()) - expect_error(doFit()) - expect_error(getGenAICs()) }) \ No newline at end of file diff --git a/tests/testthat/test-modelling.R b/tests/testthat/test-modelling.R new file mode 100644 index 0000000..e398d28 --- /dev/null +++ b/tests/testthat/test-modelling.R @@ -0,0 +1,27 @@ +# 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 new file mode 100644 index 0000000..732580d --- /dev/null +++ b/tests/testthat/test-plot.R @@ -0,0 +1,171 @@ +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 new file mode 100644 index 0000000..47eeb69 --- /dev/null +++ b/tests/testthat/test-posterior.R @@ -0,0 +1,43 @@ +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") +}) diff --git a/vignettes/.gitignore b/vignettes/.gitignore new file mode 100644 index 0000000..097b241 --- /dev/null +++ b/vignettes/.gitignore @@ -0,0 +1,2 @@ +*.html +*.R diff --git a/vignettes/Simulation_Example.Rmd b/vignettes/Simulation_Example.Rmd new file mode 100644 index 0000000..bcd8cef --- /dev/null +++ b/vignettes/Simulation_Example.Rmd @@ -0,0 +1,115 @@ +--- +title: "Simulation Example of Bayesian MCPMod for Continuous Data" +output: rmarkdown::html_vignette +number_sections: true +vignette: > + %\VignetteIndexEntry{Simulation Example of Bayesian MCPMod for Continuous Data} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +```{r setup} +library(BayesianMCPMod) +library(clinDR) +library(dplyr) + +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. +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. + + +```{r Calculation of a MAP prior} +data("metaData") +testdata <- as.data.frame(metaData) +dataset <- filter(testdata, bname == "BRINTELLIX") +histcontrol <- filter(dataset, dose == 0, primtime == 8, indication == "MAJOR DEPRESSIVE DISORDER") + +##Create MAP Prior +hist_data <- data.frame( + trial = histcontrol$nctno, + est = histcontrol$rslt, + se = histcontrol$se, + sd = histcontrol$sd, + n = histcontrol$sampsize) + +dose_levels <- c(0, 2.5, 5, 10, 20) + +prior_list <- getPriorList( + hist_data = hist_data, + dose_levels = dose_levels) + +``` + +# 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. +```{r} +#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") +sigemax<- DoseFinding::guesst(d = c(2.5,5), + p = c(0.1,0.6), + model = "sigEmax") +#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, + #betaMod = beta, + doses = dose_levels, + maxEff = -3, + placEff = -12.8) + +n_patients=c(60,80,80,80,80) +``` + +# Calculation of success probabilities + +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 +``` + +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 +``` + +For this specific trial setting the adapted allocation ratio leads to increased success probabilities under all assumed dose response relationships. + + diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd new file mode 100644 index 0000000..e84ebac --- /dev/null +++ b/vignettes/analysis_normal.Rmd @@ -0,0 +1,179 @@ +--- +title: "Analysis Example of Bayesian MCPMod for Continuous Data" +output: rmarkdown::html_vignette +number_sections: true +vignette: > + %\VignetteIndexEntry{Analysis Example of Bayesian MCPMod for Continuous Data} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +```{r setup} +library(BayesianMCPMod) +library(clinDR) +library(dplyr) + +set.seed(7015) +``` +# Background and data + +In this vignette we will show the use of the Bayesian MCPMod package for continuous distributed data. +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... + + +# 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} +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) + +##Create MAP Prior +hist_data <- data.frame( + trial = histcontrol$nctno, + est = histcontrol$rslt, + se = histcontrol$se, + sd = histcontrol$sd, + n = histcontrol$sampsize) + +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. + +```{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") + + +mods <- DoseFinding::Mods( + linear = NULL, + emax = emax, + exponential = exp, + doses = dose_levels, + maxEff = -1, + placEff = -12.8) + +new_trial <- filter(dataset, primtime == 8, indication == "MAJOR DEPRESSIVE DISORDER",protid==6) +n_patients <- c(150, 142, 147, 149) + +``` + +# 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. + +```{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"), + mods = mods, + true_model = "emax") + +posterior <- getPosterior(prior=prior_list,data=data_emax, + mu_hat = new_trial$rslt, + sd_hat = new_trial$se) + +``` + +# 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. + +```{r Execution of Bayesian MCPMod Test step} +crit_pval <- getCritProb( + mods = mods, + dose_levels = dose_levels, + dose_weights = n_patients, + alpha_crit_val = 0.1) + +contr_mat_prior <- getContrMat( + 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( + posteriors_list = 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 + +BMCP_result +``` +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. + +# 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. + +```{r Model fitting} +#Model fit +post_observed <- posterior +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) + +plot(fit,cr_bands=TRUE) +plot(fit_simple, cr_bands = TRUE) + +``` +# Additional notes +TBD, whether certain wrapper functions should be mentioned. + +