From 23460355dd9b8b758f2eca435a622fdc2a41ec36 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Tue, 5 Sep 2023 16:27:15 +0200 Subject: [PATCH 01/43] - initial test commit --- R/plot.R | 53 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/R/plot.R b/R/plot.R index e69de29..d15ca48 100644 --- a/R/plot.R +++ b/R/plot.R @@ -0,0 +1,53 @@ +plot.estMod <- function ( + + est_mod + # dose_levels + # posteriors = posterior_linear[[2]] + +) { + + model <- est_mod$model + theta <- est_mod$fit$solution + + dose <- seq(min(dose_levels), max(dose_levels), length.out = 1e3) + + switch(model, + "emax" = { + resp_expr <- quote(theta[1] + (theta[2] * dose) / (theta[3] + dose))}, + "sigEmax" = { + resp_expr <- quote(theta[1] + (theta[2] * dose^theta[4]) / (theta[3]^theta[4] + dose^theta[4]))}, + "exponential" = { + resp_expr <- quote(theta[1] + theta[2] * (exp(dose / theta[3]) - 1))}, + "quadratic" = { + resp_expr <- quote(theta[1] + theta[2] * dose + theta[3] * dose^2)}, + "linear" = { + resp_expr <- quote(theta[1] + theta[2] * dose)}, + "logistic" = { + resp_expr <- quote(theta[1] + theta[2] / (1 + exp((theta[3] - dose) / theta[4])))}, + { + stop(GENERAL$ERROR$MODEL_OPTIONS)} + ) + + df <- data.frame(dose = dose, response = eval(resp_expr)) + + data.frame(dose = dose_levels, obs = ) + + plt <- ggplot2::ggplot(data = df) + + ggplot2::geom_line(ggplot2::aes(dose, response)) + + ggplot2::geom_point() + + return(plt) + +} + +plot.estMods <- function ( + + est_mods + +) { + + plts <- lapply(est_mods, plot) + + + +} From d2cf056738f7b213216cd6b5425dd1967fa74845 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Fri, 15 Sep 2023 17:42:12 +0200 Subject: [PATCH 02/43] - feature restructuring done --- R/BMCPMod.R | 73 +++++-- R/dataSimulation.R | 43 ----- R/helper.R | 144 +++----------- R/modelling.R | 470 ++++++++++++++++++++++++--------------------- R/postShape.R | 26 --- R/posterior.R | 106 ++++++++++ R/simulation.R | 32 +++ 7 files changed, 469 insertions(+), 425 deletions(-) delete mode 100644 R/dataSimulation.R delete mode 100644 R/postShape.R create mode 100644 R/posterior.R create mode 100644 R/simulation.R diff --git a/R/BMCPMod.R b/R/BMCPMod.R index c9215cf..50a3bf9 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -1,26 +1,61 @@ #' @title BMCPMod -#' @param ancova1 tbd -#' -#' @param cont_Mat1 tbd +#' +#' @param posteriors_list tbd +#' @param contr_mat tbd #' @param crit_prob tbd -#' @param n_simulations tbd -#' +#' #' @export BMCPMod <- function( - ancova1, - cont_Mat1, - crit_prob, - n_simulations) { - adj1_p <- list() + + posteriors_list, + contr_mat, + crit_prob + +) { + + lapply(posteriors_list, BayesMCPMod, contr_mat, crit_prob) + +} - for (i in 1:n_simulations) { - ancova <- ancova1[[i]] - adj1_p[[i]] <- BayesMCPMod( - ancova, - cont_Mat1, - crit_prob - ) +BayesMCPMod <- 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) + } - - return(adj1_p) + + post_combs_i <- getPostCombsI(posterior_i) + post_probs <- apply(contr_mat$contMat, 2, getPostProb, post_combs_i) + + res_list <- list( + sign = ifelse(max_prob > crit_prob, 1, 0), + p_val = max(post_probs), + post_probs = post_probs) + + return (res_list) + } 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 index d5aeab1..97658a7 100644 --- a/R/helper.R +++ b/R/helper.R @@ -1,116 +1,34 @@ -#' @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) +## from DoseFinding +powCalc <- function(alternative, critV, df, corMat, deltaMat, control) { + mvtnorm.control <- DoseFinding::mvtnorm.control + + pmvt <- mvtnorm::pmvt + + nC <- nrow(corMat) + if (alternative[1] == "two.sided") { + lower <- rep(-critV, nC) + } else { + lower <- rep(-Inf, nC) } - - max_prob <- max(probs) - - if (max_prob > crit_prob) { - sign <- 1 + upper <- rep(critV, nC) + if (!missing(control)) { + if (!is.list(control)) { + stop("when specified, 'control' must be a list") + } + ctrl <- do.call("mvtnorm.control", control) + } else { + ctrl <- control } - - 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]) + ctrl$interval <- NULL + nScen <- ncol(deltaMat) + res <- numeric(nScen) + for (i in 1:nScen) { + pmvtCall <- c(list(lower, upper, + df = df, corr = corMat, + delta = deltaMat[, i], algorithm = ctrl + )) + res[i] <- as.vector(1 - do.call("pmvt", pmvtCall)) } - - 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 + names(res) <- colnames(deltaMat) + res +} diff --git a/R/modelling.R b/R/modelling.R index e6abc27..22a33e6 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -1,238 +1,260 @@ -#' @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) - - sublist[["model"]] <- shape - sublist[["fit"]] <- fit$coefs - sublist[["pred"]] <- fit_vec - sublist[["gAIC"]] <- gAIC - sublist[["max_effect"]] <- max_effect - - summary[[j]] <- sublist +getModelFits <- function ( + + models, + dose_levels, + posterior, + simple = FALSE + +) { + + if (simple) { + + model_fits <- lapply(models, getModelFitSimple, dose_levels, posterior) + + } else { + + model_fits <- lapply(models, getModelFitOpt, dose_levels, posterior) + } - - return(summary) + + model_fits <- addModelWeights(model_fits) + names(model_fits) <- models + + attr(model_fits, "posterior") <- posterior + class(model_fits) <- "modelFits" + + return (model_fits) + } -#' @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 - ) +## 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 = doFit_ID$VALUES$FIT_TYPE, + 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) + +} - return(list( - model = model, - fit = res - )) +## 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) + } - - results <- lapply(models, function(model) { - fitModel(model, post, doses) - }) - - return(results) + + 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 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) +predictModelFit <- function ( + + model_fit, + doses = NULL + +) { + + if (is.null(doses)) { + + doses <- model_fit$dose_levels + } + + 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) + +} - results <- lapply(fit_mods, function(fit_mod) { - getGenAIC(post, fit_mod, doses) - }) - - 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 - } +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) + +} - return(results) +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) + } 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..0ed309e --- /dev/null +++ b/R/posterior.R @@ -0,0 +1,106 @@ +#' @title getPosterior +#' +#' @param data tbd +#' @param prior prior_list +#' +#' @export +getPosterior <- function( + + data, + prior_list + +) { + + lapply(split(data, data$simulation), getPosteriorI, prior_list = prior_list) + +} + +getPosteriorI <- function( + + data_i, + prior_list + +) { + + anova_res <- lm(data_i$response ~ factor(data_i$dose) - 1) # take mean & var out in separate function + anova_mean <- summary(anova_res)$coefficients[, 1] # + anova_se <- summary(anova_res)$coefficients[, 2] # + + post_list <- mapply(RBesT::postmix, prior_list, m = anova_mean, se = anova_se) + + names(post_list) <- c("Ctr", paste0("DG_", seq_along(post_list[-1]))) + 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) + +} + +summary.postList <- function ( + + post_list + +) { + + summary_list <- lapply(post_list, summary) + names(summary_list) <- names(post_list) + summary_tab <- do.call(rbind, summary_list) + + return (summary_tab) + +} + +print.postList <- function ( + + post_list + +) { + + 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(post_list) + + names(post_list) <- rownames(summary_tab) + class(post_list) <- NULL + + list_out <- list(summary_tab, getMaxDiff(summary_tab[, 4]), post_list) + 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..c8ee63e --- /dev/null +++ b/R/simulation.R @@ -0,0 +1,32 @@ +#' @title simulateData +#' +#' @param n_patients tbd +#' @param dose_levels tbd +#' @param sd tbd +#' @param mods tbd +#' @param n_sim tbd +#' +#' @export +simulateData <- function( + + n_patients, + dose_levels, + sd, + mods, + n_sim + +) { + + 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) + + return (sim_data) + +} From ad2a638c3a7ba51150e4ac82ea09aa54c4b1326f Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Fri, 15 Sep 2023 17:53:45 +0200 Subject: [PATCH 03/43] - removed helper.R, became superflouos --- R/helper.R | 34 ---------------------------------- 1 file changed, 34 deletions(-) delete mode 100644 R/helper.R diff --git a/R/helper.R b/R/helper.R deleted file mode 100644 index 97658a7..0000000 --- a/R/helper.R +++ /dev/null @@ -1,34 +0,0 @@ -## from DoseFinding -powCalc <- function(alternative, critV, df, corMat, deltaMat, control) { - mvtnorm.control <- DoseFinding::mvtnorm.control - - pmvt <- mvtnorm::pmvt - - nC <- nrow(corMat) - if (alternative[1] == "two.sided") { - lower <- rep(-critV, nC) - } else { - lower <- rep(-Inf, nC) - } - upper <- rep(critV, nC) - if (!missing(control)) { - if (!is.list(control)) { - stop("when specified, 'control' must be a list") - } - ctrl <- do.call("mvtnorm.control", control) - } else { - ctrl <- control - } - ctrl$interval <- NULL - nScen <- ncol(deltaMat) - res <- numeric(nScen) - for (i in 1:nScen) { - pmvtCall <- c(list(lower, upper, - df = df, corr = corMat, - delta = deltaMat[, i], algorithm = ctrl - )) - res[i] <- as.vector(1 - do.call("pmvt", pmvtCall)) - } - names(res) <- colnames(deltaMat) - res -} From 6e12c0c745ad1003c55a75b653012a5db5fafbce Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Fri, 15 Sep 2023 17:56:21 +0200 Subject: [PATCH 04/43] - updated documentation --- NAMESPACE | 6 ++---- man/BMCPMod.Rd | 8 +++----- man/BayesMCPMod.Rd | 21 --------------------- man/doAncova.Rd | 21 --------------------- man/estimateModels.Rd | 18 ------------------ man/getGenAICs.Rd | 18 ------------------ man/{doFit.Rd => getModelFits.Rd} | 16 +++++++++------- man/getPosterior.Rd | 16 ++++++++++++++++ man/postShape.Rd | 18 ------------------ man/simulateData.Rd | 25 +++++-------------------- 10 files changed, 35 insertions(+), 132 deletions(-) delete mode 100644 man/BayesMCPMod.Rd delete mode 100644 man/doAncova.Rd delete mode 100644 man/estimateModels.Rd delete mode 100644 man/getGenAICs.Rd rename man/{doFit.Rd => getModelFits.Rd} (50%) create mode 100644 man/getPosterior.Rd delete mode 100644 man/postShape.Rd diff --git a/NAMESPACE b/NAMESPACE index 2e76e4c..306d25f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,8 +1,6 @@ # Generated by roxygen2: do not edit by hand export(BMCPMod) -export(doFit) -export(estimateModels) -export(getGenAICs) -export(postShape) +export(getModelFits) +export(getPosterior) export(simulateData) diff --git a/man/BMCPMod.Rd b/man/BMCPMod.Rd index 1d82663..874fc43 100644 --- a/man/BMCPMod.Rd +++ b/man/BMCPMod.Rd @@ -4,16 +4,14 @@ \alias{BMCPMod} \title{BMCPMod} \usage{ -BMCPMod(ancova1, cont_Mat1, crit_prob, n_simulations) +BMCPMod(posteriors_list, contr_mat, crit_prob) } \arguments{ -\item{ancova1}{tbd} +\item{posteriors_list}{tbd} -\item{cont_Mat1}{tbd} +\item{contr_mat}{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/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/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..200d3b5 --- /dev/null +++ b/man/getPosterior.Rd @@ -0,0 +1,16 @@ +% 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) +} +\arguments{ +\item{data}{tbd} + +\item{prior}{prior_list} +} +\description{ +getPosterior +} 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..b3572d5 100644 --- a/man/simulateData.Rd +++ b/man/simulateData.Rd @@ -1,36 +1,21 @@ % 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} \usage{ -simulateData( - n_patients, - dose_levels, - pcov, - SD, - placebo_effect, - max_effect, - mod, - n_simulations -) +simulateData(n_patients, dose_levels, sd, mods, n_sim) } \arguments{ \item{n_patients}{tbd} \item{dose_levels}{tbd} -\item{pcov}{tbd} +\item{sd}{tbd} -\item{SD}{tbd} +\item{mods}{tbd} -\item{placebo_effect}{tbd} - -\item{max_effect}{tbd} - -\item{mod}{tbd} - -\item{n_simulations}{tbd} +\item{n_sim}{tbd} } \description{ simulateData From 602dd76e9f2a04960bfdf3dffe1226a45e689c70 Mon Sep 17 00:00:00 2001 From: Bossert Date: Fri, 22 Sep 2023 11:47:24 +0200 Subject: [PATCH 05/43] Create dummy vignette --- vignettes/.gitignore | 2 ++ vignettes/analysis_normal.Rmd | 25 +++++++++++++++++++++++++ 2 files changed, 27 insertions(+) create mode 100644 vignettes/.gitignore create mode 100644 vignettes/analysis_normal.Rmd 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/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd new file mode 100644 index 0000000..f343beb --- /dev/null +++ b/vignettes/analysis_normal.Rmd @@ -0,0 +1,25 @@ +--- +title: "Analysis Example of Bayesian MCPMod for Continuous Data" +output: rmarkdown::html_vignette +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) +``` +This vignette shows how to apply the BayesianMCPMod package +```{r} +##Create MAP Prior +``` + + From 80a80ad469b92b99b4863e6c4cf8b297b2e27ef6 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Fri, 22 Sep 2023 13:21:04 +0200 Subject: [PATCH 06/43] - small bug fix for significance test --- R/BMCPMod.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 50a3bf9..abdc170 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -52,7 +52,7 @@ BayesMCPMod <- function ( post_probs <- apply(contr_mat$contMat, 2, getPostProb, post_combs_i) res_list <- list( - sign = ifelse(max_prob > crit_prob, 1, 0), + sign = ifelse(max(post_probs) > crit_prob, 1, 0), p_val = max(post_probs), post_probs = post_probs) From 6b2beb09496effdfe9aaedfa9fa2e1293aca71dd Mon Sep 17 00:00:00 2001 From: Andersen Date: Mon, 25 Sep 2023 13:14:19 +0200 Subject: [PATCH 07/43] fixed warnings in cmd check --- BayesianMCPMod.Rproj | 4 ++++ DESCRIPTION | 1 + R/posterior.R | 2 +- man/getPosterior.Rd | 2 +- 4 files changed, 7 insertions(+), 2 deletions(-) diff --git a/BayesianMCPMod.Rproj b/BayesianMCPMod.Rproj index 8e3c2eb..21a4da0 100644 --- a/BayesianMCPMod.Rproj +++ b/BayesianMCPMod.Rproj @@ -11,3 +11,7 @@ Encoding: UTF-8 RnwWeave: Sweave LaTeX: pdfLaTeX + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source diff --git a/DESCRIPTION b/DESCRIPTION index 034bfc9..af2a0c9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -22,6 +22,7 @@ Roxygen: list(markdown = TRUE) RoxygenNote: 7.2.3 Imports: DoseFinding, + ggplot2, stats, RBesT, nloptr diff --git a/R/posterior.R b/R/posterior.R index 0ed309e..a47e1ff 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -1,7 +1,7 @@ #' @title getPosterior #' #' @param data tbd -#' @param prior prior_list +#' @param prior_list prior_list #' #' @export getPosterior <- function( diff --git a/man/getPosterior.Rd b/man/getPosterior.Rd index 200d3b5..95c54d4 100644 --- a/man/getPosterior.Rd +++ b/man/getPosterior.Rd @@ -9,7 +9,7 @@ getPosterior(data, prior_list) \arguments{ \item{data}{tbd} -\item{prior}{prior_list} +\item{prior_list}{prior_list} } \description{ getPosterior From 8467e278c59b02925bbca8b24a2ac193428ade52 Mon Sep 17 00:00:00 2001 From: Bossert Date: Fri, 29 Sep 2023 14:48:13 +0200 Subject: [PATCH 08/43] First draft vignette - including only test step --- vignettes/analysis_normal.Rmd | 145 ++++++++++++++++++++++++++++++++++ 1 file changed, 145 insertions(+) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index f343beb..5af387a 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -16,10 +16,155 @@ knitr::opts_chunk$set( ```{r setup} library(BayesianMCPMod) +library(clinDR) +library(dplyr) ``` This vignette shows how to apply the BayesianMCPMod package ```{r} +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]))) +``` + +```{r} +#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) ``` +```{r} +#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) + +``` + + +```{r} +#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 <- BMCPMod( + posteriors_list = posterior_emax, + contr_mat = contr_mat_prior, + crit_prob = crit_pval) + +BMCP_result +``` + + +```{r} +#Model fit +#This part is currently not working +#post_observed <- posterior_emax[[1]] +#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) + +#fit1 <- DoseFinding::fitMod( +# dose = dose_levels, +# resp = fit_simple$linear$pred_values, + # S = diag(summary.postList(post_observed)[, 2]^2), +# model = "linear", + # type = "general") +#plot(fit1, + # plotData = "meansCI", + # CI = TRUE) +``` From 585a2153c585c4ff1bcc32db37d432e9fca205dd Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Wed, 4 Oct 2023 17:48:08 +0200 Subject: [PATCH 09/43] - included plot function --- R/modelling.R | 58 +++++++++++---------- R/plot.R | 138 ++++++++++++++++++++++++++++++++++---------------- R/posterior.R | 5 +- 3 files changed, 130 insertions(+), 71 deletions(-) diff --git a/R/modelling.R b/R/modelling.R index 22a33e6..2e66e4e 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -49,7 +49,7 @@ getModelFitSimple <- function ( resp = summary.postList(posterior)[, 1], S = diag(summary.postList(posterior)[, 2]^2), model = model, - type = doFit_ID$VALUES$FIT_TYPE, + type = "general", bnds = DoseFinding::defBnds(mD = max(dose_levels))[[model]]) pred_vals <- stats::predict(fit, predType = "ls-means") @@ -190,31 +190,37 @@ predictModelFit <- function ( 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"])}, + "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) diff --git a/R/plot.R b/R/plot.R index d15ca48..4d1939c 100644 --- a/R/plot.R +++ b/R/plot.R @@ -1,53 +1,105 @@ -plot.estMod <- function ( +plot.modelFits <- function ( - est_mod - # dose_levels - # posteriors = posterior_linear[[2]] + model_fits, + CrI = FALSE, + gAIC = TRUE, + avg_fit = TRUE ) { - model <- est_mod$model - theta <- est_mod$fit$solution - - dose <- seq(min(dose_levels), max(dose_levels), length.out = 1e3) - - switch(model, - "emax" = { - resp_expr <- quote(theta[1] + (theta[2] * dose) / (theta[3] + dose))}, - "sigEmax" = { - resp_expr <- quote(theta[1] + (theta[2] * dose^theta[4]) / (theta[3]^theta[4] + dose^theta[4]))}, - "exponential" = { - resp_expr <- quote(theta[1] + theta[2] * (exp(dose / theta[3]) - 1))}, - "quadratic" = { - resp_expr <- quote(theta[1] + theta[2] * dose + theta[3] * dose^2)}, - "linear" = { - resp_expr <- quote(theta[1] + theta[2] * dose)}, - "logistic" = { - resp_expr <- quote(theta[1] + theta[2] / (1 + exp((theta[3] - dose) / theta[4])))}, - { - stop(GENERAL$ERROR$MODEL_OPTIONS)} - ) - - df <- data.frame(dose = dose, response = eval(resp_expr)) - - data.frame(dose = dose_levels, obs = ) - - plt <- ggplot2::ggplot(data = df) + - ggplot2::geom_line(ggplot2::aes(dose, response)) + - ggplot2::geom_point() - - return(plt) - -} - -plot.estMods <- function ( + plot_resolution <- 1e3 + + dose_levels <- model_fits[[1]]$dose_levels + post_summary <- summary.postList(attr(model_fits, "posterior")) + doses <- seq(from = min(dose_levels), + to = max(dose_levels), length.out = plot_resolution) + + preds_models <- sapply(model_fits, predictModelFit, doses = doses) + model_names <- names(model_fits) + + if (avg_fit) { + + mod_weigts <- sapply(model_fits, function (x) x$model_weight) + avg_mod <- preds_models %*% mod_weigts + + preds_models <- cbind(preds_models, avg_mod) + model_names <- c(model_names, "averageModel") - est_mods + } -) { + gg_data <- data.frame( + dose_levels = rep(doses, length(model_names)), + fits = as.vector(preds_models), + models = rep(factor(model_names, + levels = c("linear", "emax", "exponential", + "sigEmax", "logistic", "quadratic", + "averageModel")), + each = plot_resolution)) - plts <- lapply(est_mods, plot) + if (gAIC) { + + g_AICs <- sapply(model_fits, function (x) x$gAIC) + label_gAUC <- paste("AIC:", round(g_AICs), digits = 1) + + if (avg_fit) { + + ## replace with pipe once native pipe => becomes available + paste_names <- names(model_fits) |> + 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_weigts, 1), + collapse = ", ") + label_gAUC <- c(label_gAUC, label_avg) + + } + + } + 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)} + } + + ## Posterior Credible Intervals + {if (CrI) { + 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, + fits = post_summary[, 4]), + mapping = ggplot2::aes(dose_levels, fits), + size = 2) + + ## Fitted Models + ggplot2::geom_line( + data = gg_data, + mapping = ggplot2::aes(dose_levels, fits)) + + ## Faceting + ggplot2::facet_wrap(~ models) + return (plts) -} +} \ No newline at end of file diff --git a/R/posterior.R b/R/posterior.R index 0ed309e..1d66001 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -55,11 +55,12 @@ getPostCombsI <- function ( summary.postList <- function ( - post_list + post_list, + ... ) { - summary_list <- lapply(post_list, summary) + summary_list <- lapply(post_list, summary, ...) names(summary_list) <- names(post_list) summary_tab <- do.call(rbind, summary_list) From 644ec5ce9d23a98a8bd24a45b4af788598c81b53 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Wed, 4 Oct 2023 17:54:44 +0200 Subject: [PATCH 10/43] -fixed typo in plot function --- R/plot.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/plot.R b/R/plot.R index 4d1939c..5821af3 100644 --- a/R/plot.R +++ b/R/plot.R @@ -7,7 +7,7 @@ plot.modelFits <- function ( ) { - plot_resolution <- 1e3 + plot_resolution <- 1e3 dose_levels <- model_fits[[1]]$dose_levels post_summary <- summary.postList(attr(model_fits, "posterior")) @@ -39,7 +39,7 @@ plot.modelFits <- function ( if (gAIC) { g_AICs <- sapply(model_fits, function (x) x$gAIC) - label_gAUC <- paste("AIC:", round(g_AICs), digits = 1) + label_gAUC <- paste("AIC:", round(g_AICs, digits = 1)) if (avg_fit) { From cee4363ec937fe55ce523a7eda013e05f988a7fd Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Wed, 4 Oct 2023 19:00:15 +0200 Subject: [PATCH 11/43] - getPosterior can now take custom vectors for mu and sd - other small changes --- R/modelling.R | 4 ++-- R/plot.R | 8 ++++---- R/posterior.R | 36 +++++++++++++++++++++++++++++------- 3 files changed, 35 insertions(+), 13 deletions(-) diff --git a/R/modelling.R b/R/modelling.R index 2e66e4e..635436a 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -17,11 +17,11 @@ getModelFits <- function ( if (simple) { - model_fits <- lapply(models, getModelFitSimple, dose_levels, posterior) + model_fits <- lapply(models, getModelFitSimple, dose_levels, posterior[[1]]) } else { - model_fits <- lapply(models, getModelFitOpt, dose_levels, posterior) + model_fits <- lapply(models, getModelFitOpt, dose_levels, posterior[[1]]) } diff --git a/R/plot.R b/R/plot.R index 5821af3..3179773 100644 --- a/R/plot.R +++ b/R/plot.R @@ -43,16 +43,16 @@ plot.modelFits <- function ( if (avg_fit) { - ## replace with pipe once native pipe => becomes available - paste_names <- names(model_fits) |> + mod_weigts <- sort(mod_weigts, decreasing = TRUE) + paste_names <- names(mod_weigts) |> 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_weigts, 1), - collapse = ", ") + label_avg <- paste0(paste_names, "=", round(mod_weigts, 1), + collapse = ", ") label_gAUC <- c(label_gAUC, label_avg) } diff --git a/R/posterior.R b/R/posterior.R index 1d66001..1ef955d 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -7,26 +7,48 @@ getPosterior <- function( data, - prior_list + prior_list, + mu_hat = NULL, + sd_hat = NULL ) { - lapply(split(data, data$simulation), getPosteriorI, prior_list = prior_list) + lapply(split(data, data$simulation), getPosteriorI, + prior_list = prior_list, + mu_hat = mu_hat, + sd_hat = sd_hat) } getPosteriorI <- function( data_i, - prior_list + prior_list, + mu_hat = NULL, + sd_hat = NULL ) { - anova_res <- lm(data_i$response ~ factor(data_i$dose) - 1) # take mean & var out in separate function - anova_mean <- summary(anova_res)$coefficients[, 1] # - anova_se <- summary(anova_res)$coefficients[, 2] # + if (is.null(mu_hat) && is.null(sd_hat)) { + + anova_res <- 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) == nrow(sd_hat)) + + } else { + + stop ("Both mu_hat and S_hat must be provided.") + + } - post_list <- mapply(RBesT::postmix, prior_list, m = anova_mean, se = anova_se) + post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = sd_hat) names(post_list) <- c("Ctr", paste0("DG_", seq_along(post_list[-1]))) class(post_list) <- "postList" From 63dfec0fe4eb2fd63a4759bc92a5b8901dc21c57 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Wed, 4 Oct 2023 19:18:33 +0200 Subject: [PATCH 12/43] - for single set of simulated data (e.g. observed data), getPosterior returns the posterior directly, not the list --- R/modelling.R | 4 ++-- R/plot.R | 10 ++++++++-- R/posterior.R | 16 ++++++++++++---- 3 files changed, 22 insertions(+), 8 deletions(-) diff --git a/R/modelling.R b/R/modelling.R index 635436a..2e66e4e 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -17,11 +17,11 @@ getModelFits <- function ( if (simple) { - model_fits <- lapply(models, getModelFitSimple, dose_levels, posterior[[1]]) + model_fits <- lapply(models, getModelFitSimple, dose_levels, posterior) } else { - model_fits <- lapply(models, getModelFitOpt, dose_levels, posterior[[1]]) + model_fits <- lapply(models, getModelFitOpt, dose_levels, posterior) } diff --git a/R/plot.R b/R/plot.R index 3179773..64d675c 100644 --- a/R/plot.R +++ b/R/plot.R @@ -68,16 +68,20 @@ plot.modelFits <- function ( 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)} + size = 3) + + } } + ## Posterior Credible Intervals {if (CrI) { + ggplot2::geom_errorbar( data = data.frame(x = dose_levels, ymin = post_summary[, 3], @@ -85,7 +89,9 @@ plot.modelFits <- function ( mapping = ggplot2::aes(x = x, ymin = ymin, ymax = ymax), - width = 0, alpha = 0.5)} + width = 0, alpha = 0.5) + + } } + ## Posterior Medians ggplot2::geom_point( diff --git a/R/posterior.R b/R/posterior.R index 1ef955d..760f297 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -13,10 +13,18 @@ getPosterior <- function( ) { - lapply(split(data, data$simulation), getPosteriorI, - prior_list = prior_list, - mu_hat = mu_hat, - sd_hat = sd_hat) + 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) } From 4222181b728ab167557c484dccd9e367775a14f9 Mon Sep 17 00:00:00 2001 From: Bossert Date: Thu, 5 Oct 2023 09:56:25 +0200 Subject: [PATCH 13/43] Inclusion of model fits. --- vignettes/analysis_normal.Rmd | 37 ++++++++++++++++------------------- 1 file changed, 17 insertions(+), 20 deletions(-) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 5af387a..e6ab219 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -146,25 +146,22 @@ BMCP_result ```{r} #Model fit #This part is currently not working -#post_observed <- posterior_emax[[1]] -#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) - -#fit1 <- DoseFinding::fitMod( -# dose = dose_levels, -# resp = fit_simple$linear$pred_values, - # S = diag(summary.postList(post_observed)[, 2]^2), -# model = "linear", - # type = "general") -#plot(fit1, - # plotData = "meansCI", - # CI = TRUE) +post_observed <- posterior_emax[[1]] +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) ``` From 5e4327744f1587d9da89ac5485b28cfd8e680561 Mon Sep 17 00:00:00 2001 From: Bossert Date: Thu, 5 Oct 2023 12:59:09 +0200 Subject: [PATCH 14/43] Addition of (minimum) explanation text to vignette --- vignettes/analysis_normal.Rmd | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index e6ab219..2ed8389 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -19,7 +19,10 @@ library(BayesianMCPMod) library(clinDR) library(dplyr) ``` -This vignette shows how to apply the BayesianMCPMod package +This vignette shows how to apply the BayesianMCPMod package. + +# i) Calculation of a MAP prior +In a first step a meta analytic prior will be calculated using historical data. ```{r} data("metaData") testdata <- as.data.frame(metaData) @@ -64,7 +67,7 @@ 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]))) ``` - +# ii) Pre-Specification of candidate models ```{r} #Pre-Specification (B)MCPMod @@ -87,7 +90,7 @@ mods <- DoseFinding::Mods( maxEff = 10, placEff = 1.4) ``` - +# iii) Trial results ```{r} #Simulation of new trial ##Note: This part will be simplified and direct results from one trial will be used @@ -114,7 +117,7 @@ posterior_emax <- getPosterior( ``` - +# iv) Execution of Bayesian MCPMod Test step ```{r} #Evaluation of Bayesian MCPMod @@ -142,7 +145,8 @@ BMCP_result <- BMCPMod( BMCP_result ``` - +# v) Model fitting and Plotting +In the model fitting step the posterior distribution is used as basis. ```{r} #Model fit #This part is currently not working @@ -165,3 +169,8 @@ fit<- getModelFits( simple = FALSE) ``` +```{r} +plot(fit,CrI= TRUE) +plot(fit_simple, CrI= TRUE) +``` +For the plotting the credible intervals are shown as well. From 13ebcb20fe22b0195bf23293487c53de2d263923 Mon Sep 17 00:00:00 2001 From: Bossert Date: Fri, 6 Oct 2023 14:53:53 +0200 Subject: [PATCH 15/43] Correction of two minor bugs in the vignette file. --- vignettes/analysis_normal.Rmd | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 2ed8389..57a74f6 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -138,7 +138,7 @@ contr_mat_prior <- DoseFinding::optContr( w = n_patients + ess_prior) BMCP_result <- BMCPMod( - posteriors_list = posterior_emax, + posteriors_list = list(posterior_emax), contr_mat = contr_mat_prior, crit_prob = crit_pval) @@ -150,7 +150,7 @@ In the model fitting step the posterior distribution is used as basis. ```{r} #Model fit #This part is currently not working -post_observed <- posterior_emax[[1]] +post_observed <- posterior_emax model_shapes <- c("linear", "emax", "exponential") @@ -170,7 +170,7 @@ fit<- getModelFits( ``` ```{r} -plot(fit,CrI= TRUE) -plot(fit_simple, CrI= TRUE) +plot.modelFits(fit,CrI= TRUE) +plot.modelFits(fit_simple, CrI= TRUE) ``` For the plotting the credible intervals are shown as well. From 7dc310b06ce5ea56d0a5a3100b56c4f5948aae10 Mon Sep 17 00:00:00 2001 From: Andersen Date: Mon, 9 Oct 2023 14:46:35 +0200 Subject: [PATCH 16/43] fixed github action errors --- DESCRIPTION | 9 +++++++-- NAMESPACE | 1 + R/plot.R | 11 ++++++++++- R/posterior.R | 3 ++- man/getPosterior.Rd | 6 +++++- man/plot_modelFits.Rd | 23 +++++++++++++++++++++++ vignettes/analysis_normal.Rmd | 4 ++-- 7 files changed, 50 insertions(+), 7 deletions(-) create mode 100644 man/plot_modelFits.Rd diff --git a/DESCRIPTION b/DESCRIPTION index af2a0c9..e699523 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -25,9 +25,14 @@ Imports: ggplot2, stats, RBesT, - nloptr + nloptr, + clinDR, + knitr 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 306d25f..75af723 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,4 +3,5 @@ export(BMCPMod) export(getModelFits) export(getPosterior) +export(plot_modelFits) export(simulateData) diff --git a/R/plot.R b/R/plot.R index 64d675c..4ccdfc1 100644 --- a/R/plot.R +++ b/R/plot.R @@ -1,4 +1,13 @@ -plot.modelFits <- function ( +#' @title plot_modelFits +#' +#' @param model_fits tbd +#' @param CrI tbd +#' @param gAIC tbd +#' @param avg_fit tbd +#' +#' @return tbd +#' @export +plot_modelFits <- function ( model_fits, CrI = FALSE, diff --git a/R/posterior.R b/R/posterior.R index 42d1026..afa4f80 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -2,10 +2,11 @@ #' #' @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, diff --git a/man/getPosterior.Rd b/man/getPosterior.Rd index 95c54d4..501f5de 100644 --- a/man/getPosterior.Rd +++ b/man/getPosterior.Rd @@ -4,12 +4,16 @@ \alias{getPosterior} \title{getPosterior} \usage{ -getPosterior(data, prior_list) +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/plot_modelFits.Rd b/man/plot_modelFits.Rd new file mode 100644 index 0000000..4b0b7a2 --- /dev/null +++ b/man/plot_modelFits.Rd @@ -0,0 +1,23 @@ +% 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{ +plot_modelFits(model_fits, CrI = FALSE, gAIC = TRUE, avg_fit = TRUE) +} +\arguments{ +\item{model_fits}{tbd} + +\item{CrI}{tbd} + +\item{gAIC}{tbd} + +\item{avg_fit}{tbd} +} +\value{ +tbd +} +\description{ +plot_modelFits +} diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 57a74f6..60069f9 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -170,7 +170,7 @@ fit<- getModelFits( ``` ```{r} -plot.modelFits(fit,CrI= TRUE) -plot.modelFits(fit_simple, CrI= TRUE) +plot_modelFits(fit,CrI= TRUE) +plot_modelFits(fit_simple, CrI= TRUE) ``` For the plotting the credible intervals are shown as well. From 497a709dff483e4fed79cb89eb4c69ddcaada387 Mon Sep 17 00:00:00 2001 From: Andersen Date: Mon, 9 Oct 2023 15:01:02 +0200 Subject: [PATCH 17/43] added missing dependency --- DESCRIPTION | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index e699523..d96cda2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -27,7 +27,8 @@ Imports: RBesT, nloptr, clinDR, - knitr + knitr, + rmarkdown Suggests: testthat (>= 3.0.0) Depends: From a5cd79382e4737dfb9bc2a71f77bb47c2c85126a Mon Sep 17 00:00:00 2001 From: Andersen Date: Mon, 9 Oct 2023 16:50:06 +0200 Subject: [PATCH 18/43] styled vignette! --- vignettes/analysis_normal.Rmd | 108 +++++++++++++++++++--------------- 1 file changed, 62 insertions(+), 46 deletions(-) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 60069f9..ccf3f91 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -26,16 +26,17 @@ In a first step a meta analytic prior will be calculated using historical data. ```{r} data("metaData") testdata <- as.data.frame(metaData) -dataset<-filter(testdata, bname == "VIAGRA") -histcontrol<-filter(dataset, dose==0, primtime==12,indication=="ERECTILE DYSFUNCTION") +dataset <- filter(testdata, bname == "VIAGRA") +histcontrol <- filter(dataset, dose == 0, primtime == 12, indication == "ERECTILE DYSFUNCTION") -##Create MAP Prior +## 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) + 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)) @@ -47,39 +48,46 @@ gmap <- RBesT::gMAP( family = gaussian, beta.prior = cbind(0, 100 * sd_tot), tau.dist = "HalfNormal", - tau.prior = cbind(0, sd_tot / 4)) + tau.prior = cbind(0, sd_tot / 4) +) -prior_ctr <- RBesT::robustify( +prior_ctr <- RBesT::robustify( priormix = RBesT::automixfit(gmap), weight = 0.5, mean = 1.4, - sigma = sd_tot) + sigma = sd_tot +) -#RBesT::ess(prior_ctr) +# 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") + 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]))) ``` # ii) Pre-Specification of candidate models ```{r} -#Pre-Specification (B)MCPMod +# 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") +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( @@ -88,46 +96,50 @@ mods <- DoseFinding::Mods( exponential = exp, doses = dose_levels, maxEff = 10, - placEff = 1.4) + placEff = 1.4 +) ``` # iii) Trial results ```{r} -#Simulation of new trial -##Note: This part will be simplified and direct results from one trial will be used +# 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) + placEff = 1.4 +) -n_patients<-c(50,50,50,50) +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) + n_sim = 1 +) -data_emax <- data[, c("simulation", "dose", "emax")] +data_emax <- data[, c("simulation", "dose", "emax")] names(data_emax)[3] <- "response" posterior_emax <- getPosterior( data = data_emax, - prior_list = prior_list) - + prior_list = prior_list +) ``` # iv) Execution of Bayesian MCPMod Test step ```{r} -#Evaluation of Bayesian MCPMod +# 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) + 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))) @@ -135,12 +147,14 @@ ess_prior <- round(unlist(lapply(prior_list, RBesT::ess))) contr_mat_prior <- DoseFinding::optContr( models = mods, doses = dose_levels, - w = n_patients + ess_prior) + w = n_patients + ess_prior +) BMCP_result <- BMCPMod( posteriors_list = list(posterior_emax), - contr_mat = contr_mat_prior, - crit_prob = crit_pval) + contr_mat = contr_mat_prior, + crit_prob = crit_pval +) BMCP_result ``` @@ -148,10 +162,10 @@ BMCP_result # v) Model fitting and Plotting In the model fitting step the posterior distribution is used as basis. ```{r} -#Model fit -#This part is currently not working +# Model fit +# This part is currently not working post_observed <- posterior_emax -model_shapes <- c("linear", "emax", "exponential") +model_shapes <- c("linear", "emax", "exponential") # Option a) Simplified approach by using frequentist idea @@ -159,18 +173,20 @@ fit_simple <- getModelFits( models = model_shapes, dose_levels = dose_levels, posterior = post_observed, - simple = TRUE) + simple = TRUE +) # Option b) Making use of the complete posterior distribution -fit<- getModelFits( +fit <- getModelFits( models = model_shapes, dose_levels = dose_levels, posterior = post_observed, - simple = FALSE) + simple = FALSE +) ``` ```{r} -plot_modelFits(fit,CrI= TRUE) -plot_modelFits(fit_simple, CrI= TRUE) +plot_modelFits(fit, CrI = TRUE) +plot_modelFits(fit_simple, CrI = TRUE) ``` For the plotting the credible intervals are shown as well. From 973c313782d4ccc82b3792554c89ccd435dd7b39 Mon Sep 17 00:00:00 2001 From: Andersen Date: Mon, 9 Oct 2023 18:32:08 +0200 Subject: [PATCH 19/43] added tests --- tests/testthat/test-modelling.R | 27 +++++ tests/testthat/test-plot.R | 175 ++++++++++++++++++++++++++++++++ tests/testthat/test-posterior.R | 41 ++++++++ 3 files changed, 243 insertions(+) create mode 100644 tests/testthat/test-modelling.R create mode 100644 tests/testthat/test-plot.R create mode 100644 tests/testthat/test-posterior.R diff --git a/tests/testthat/test-modelling.R b/tests/testthat/test-modelling.R new file mode 100644 index 0000000..2de1ba4 --- /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_is(pred_vals, "numeric") +}) diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R new file mode 100644 index 0000000..bd33665 --- /dev/null +++ b/tests/testthat/test-plot.R @@ -0,0 +1,175 @@ +# Additional test cases for plot_modelFits function + +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 <- 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 <- BMCPMod( + 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 + plot6 <- plot_modelFits(fit) + expect_is(plot6, "ggplot") + + # Test with CrI = TRUE and more models + plot7 <- plot_modelFits(fit, CrI = TRUE) + expect_is(plot7, "ggplot") + + # Test with gAIC = FALSE and more models + plot8 <- plot_modelFits(fit, gAIC = FALSE) + expect_is(plot8, "ggplot") + + # Test with avg_fit = FALSE and more models + plot9 <- plot_modelFits(fit, avg_fit = FALSE) + expect_is(plot9, "ggplot") + + # Test with all non-default parameters and more models + plot10 <- plot_modelFits(fit, CrI = TRUE, gAIC = FALSE, avg_fit = FALSE) + expect_is(plot10, "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..e089f9a --- /dev/null +++ b/tests/testthat/test-posterior.R @@ -0,0 +1,41 @@ +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_is(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_is(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_is(summary_tab, "matrix") +}) From 5c852d72a38d5d51867e70d56d6aee03fd4baa14 Mon Sep 17 00:00:00 2001 From: Andersen Date: Mon, 9 Oct 2023 20:28:30 +0200 Subject: [PATCH 20/43] added missing dependency --- DESCRIPTION | 3 ++- tests/testthat/test-global.R | 11 ----------- tests/testthat/test-plot.R | 6 +----- 3 files changed, 3 insertions(+), 17 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index d96cda2..06e6dd6 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -28,7 +28,8 @@ Imports: nloptr, clinDR, knitr, - rmarkdown + rmarkdown, + dplyr Suggests: testthat (>= 3.0.0) Depends: 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-plot.R b/tests/testthat/test-plot.R index bd33665..a6f37c0 100644 --- a/tests/testthat/test-plot.R +++ b/tests/testthat/test-plot.R @@ -1,11 +1,7 @@ -# Additional test cases for plot_modelFits function - 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") @@ -114,7 +110,7 @@ test_that("Test plot_modelFits with different model_fits input", { ) ## 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) + crit_pval <- stats::pnorm(crit_val_equal) ess_prior <- round(unlist(lapply(prior_list, RBesT::ess))) From ff4bb34e3fb443b91a9a56917761b94dbae3c94d Mon Sep 17 00:00:00 2001 From: Andersen Date: Tue, 10 Oct 2023 08:18:16 +0200 Subject: [PATCH 21/43] added covr workflow and badges --- .github/workflows/test-coverage.yaml | 32 ++++++++++++++++++++++++++++ README.md | 8 +++++++ 2 files changed, 40 insertions(+) create mode 100644 .github/workflows/test-coverage.yaml diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml new file mode 100644 index 0000000..006fee5 --- /dev/null +++ b/.github/workflows/test-coverage.yaml @@ -0,0 +1,32 @@ +# 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 + + - 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/README.md b/README.md index 66ebd6e..29de37a 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,11 @@ + +[![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) +[![R-CMD-check](https://github.com/Boehringer-Ingelheim/BayesianMCPMod/actions/workflows/check-package.yaml/badge.svg)](https://github.com/Boehringer-Ingelheim/BayesianMCPMod/actions/workflows/check-package.yaml) +[![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 ... From 209aed4e8de064dccf250e8cf54d01caef674722 Mon Sep 17 00:00:00 2001 From: Andersen Date: Tue, 10 Oct 2023 08:31:08 +0200 Subject: [PATCH 22/43] adjusted covr yaml --- .github/workflows/test-coverage.yaml | 6 ++++++ README.md | 1 - 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 006fee5..fab4029 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -22,6 +22,12 @@ jobs: - 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 diff --git a/README.md b/README.md index 29de37a..b4cc3fb 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,6 @@ [![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) -[![R-CMD-check](https://github.com/Boehringer-Ingelheim/BayesianMCPMod/actions/workflows/check-package.yaml/badge.svg)](https://github.com/Boehringer-Ingelheim/BayesianMCPMod/actions/workflows/check-package.yaml) [![Codecov test coverage](https://codecov.io/gh/Boehringer-Ingelheim/BayesianMCPMod/branch/main/graph/badge.svg)](https://codecov.io/gh/Boehringer-Ingelheim/BayesianMCPMod?branch=main) From 4893c8bd5c19e7d7810aba6066ea46035622ca39 Mon Sep 17 00:00:00 2001 From: Andersen Date: Tue, 10 Oct 2023 10:57:18 +0200 Subject: [PATCH 23/43] added pkgdown yaml --- _pkgdown.yml | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 _pkgdown.yml diff --git a/_pkgdown.yml b/_pkgdown.yml new file mode 100644 index 0000000..f51063b --- /dev/null +++ b/_pkgdown.yml @@ -0,0 +1,26 @@ +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 + +reference: +- title: "Package documentation" + contents: + - BayesianMCPMod +- title: "All functions" + contents: + - BMCPMod + - getModelFits + - getPosterior + - plot_modelsFits + - simulateData From 107b9a0af126670c9814cd6b9d992dda14724628 Mon Sep 17 00:00:00 2001 From: Andersen Date: Tue, 10 Oct 2023 15:13:22 +0200 Subject: [PATCH 24/43] removed error --- _pkgdown.yml | 3 --- 1 file changed, 3 deletions(-) diff --git a/_pkgdown.yml b/_pkgdown.yml index f51063b..2c20849 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -14,9 +14,6 @@ articles: - analysis_normal reference: -- title: "Package documentation" - contents: - - BayesianMCPMod - title: "All functions" contents: - BMCPMod From a1a34ac0f736ba37f4997b214162110ad1127f07 Mon Sep 17 00:00:00 2001 From: Andersen Date: Tue, 10 Oct 2023 15:27:19 +0200 Subject: [PATCH 25/43] alter ref --- _pkgdown.yml | 9 --------- 1 file changed, 9 deletions(-) diff --git a/_pkgdown.yml b/_pkgdown.yml index 2c20849..4cfd6f8 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -12,12 +12,3 @@ articles: navbar: ~ contents: - analysis_normal - -reference: -- title: "All functions" - contents: - - BMCPMod - - getModelFits - - getPosterior - - plot_modelsFits - - simulateData From 773e1385b04c53757bf21b15544036e53a1b00b1 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Tue, 10 Oct 2023 18:08:31 +0200 Subject: [PATCH 26/43] - interim commit to enable pulling latest updated --- R/bootstrapping.R | 17 ++++++++++++++++ R/modelling.R | 11 ++-------- R/plot.R | 23 ++++++++++++--------- R/posterior.R | 2 +- vignettes/analysis_normal.Rmd | 38 +++++++++++++++++------------------ 5 files changed, 52 insertions(+), 39 deletions(-) diff --git a/R/bootstrapping.R b/R/bootstrapping.R index e69de29..ee355a8 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -0,0 +1,17 @@ +getMCMCSample <- function ( + + model_fits, + n_samples = 1e4 + +) { + + post_list <- attr(model_fits, "posterior") + + mu_hat_samples <- sapply(post_list, RBesT::rmix, n = n_samples) + sd_hat <- summary.postList(post_observed)[, 2]^2 + + getModelFitOpt() + + return () + +} \ No newline at end of file diff --git a/R/modelling.R b/R/modelling.R index 2e66e4e..b6f7f72 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -15,15 +15,8 @@ getModelFits <- function ( ) { - if (simple) { - - model_fits <- lapply(models, getModelFitSimple, dose_levels, posterior) - - } else { - - model_fits <- lapply(models, getModelFitOpt, dose_levels, posterior) - - } + getModelFit <- ifelse(simple, getModelFitSimple, getModelFitOpt) + model_fits <- lapply(models, getModelFit, dose_levels, posterior) model_fits <- addModelWeights(model_fits) names(model_fits) <- models diff --git a/R/plot.R b/R/plot.R index 64d675c..88641f2 100644 --- a/R/plot.R +++ b/R/plot.R @@ -1,16 +1,19 @@ plot.modelFits <- function ( model_fits, - CrI = FALSE, - gAIC = TRUE, - avg_fit = TRUE + gAIC = TRUE, + avg_fit = TRUE, + CrI = FALSE, + alpha_CrI = 0.05 ) { plot_resolution <- 1e3 dose_levels <- model_fits[[1]]$dose_levels - post_summary <- summary.postList(attr(model_fits, "posterior")) + post_summary <- summary.postList( + post_list = 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_resolution) @@ -19,15 +22,15 @@ plot.modelFits <- function ( if (avg_fit) { - mod_weigts <- sapply(model_fits, function (x) x$model_weight) - avg_mod <- preds_models %*% mod_weigts + 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, "averageModel") } - gg_data <- data.frame( + gg_data <- data.frame( dose_levels = rep(doses, length(model_names)), fits = as.vector(preds_models), models = rep(factor(model_names, @@ -43,15 +46,15 @@ plot.modelFits <- function ( if (avg_fit) { - mod_weigts <- sort(mod_weigts, decreasing = TRUE) - paste_names <- names(mod_weigts) |> + 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_weigts, 1), + label_avg <- paste0(paste_names, "=", round(mod_weights, 1), collapse = ", ") label_gAUC <- c(label_gAUC, label_avg) diff --git a/R/posterior.R b/R/posterior.R index 760f297..eebc2a7 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -52,7 +52,7 @@ getPosteriorI <- function( } else { - stop ("Both mu_hat and S_hat must be provided.") + stop ("Both mu_hat and sd_hat must be provided.") } diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 2ed8389..326337c 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -25,13 +25,13 @@ This vignette shows how to apply the BayesianMCPMod package. In a first step a meta analytic prior will be calculated using historical data. ```{r} data("metaData") -testdata <- as.data.frame(metaData) -dataset<-filter(testdata, bname == "VIAGRA") -histcontrol<-filter(dataset, dose==0, primtime==12,indication=="ERECTILE DYSFUNCTION") +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), + trial = c(1, 2, 3, 4), est = histcontrol$rslt, se = histcontrol$se, sd = histcontrol$sd, @@ -49,7 +49,7 @@ gmap <- RBesT::gMAP( tau.dist = "HalfNormal", tau.prior = cbind(0, sd_tot / 4)) -prior_ctr <- RBesT::robustify( +prior_ctr <- RBesT::robustify( priormix = RBesT::automixfit(gmap), weight = 0.5, mean = 1.4, @@ -73,13 +73,13 @@ prior_list <- c(list(prior_ctr), rep(list(prior_trt), times = length(dose_levels ## 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") +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( @@ -100,7 +100,7 @@ mods_sim <- DoseFinding::Mods( maxEff = 12, placEff = 1.4) -n_patients<-c(50,50,50,50) +n_patients <- c(50, 50, 50, 50) data <- simulateData( n_patients = n_patients, dose_levels = dose_levels, @@ -126,8 +126,8 @@ contr_mat <- DoseFinding::optContr( 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) +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))) @@ -138,7 +138,7 @@ contr_mat_prior <- DoseFinding::optContr( w = n_patients + ess_prior) BMCP_result <- BMCPMod( - posteriors_list = posterior_emax, + posteriors_list = list(posterior_emax), contr_mat = contr_mat_prior, crit_prob = crit_pval) @@ -150,7 +150,7 @@ In the model fitting step the posterior distribution is used as basis. ```{r} #Model fit #This part is currently not working -post_observed <- posterior_emax[[1]] +post_observed <- posterior_emax model_shapes <- c("linear", "emax", "exponential") @@ -170,7 +170,7 @@ fit<- getModelFits( ``` ```{r} -plot(fit,CrI= TRUE) -plot(fit_simple, CrI= TRUE) +plot.modelFits(fit, CrI = TRUE) +plot.modelFits(fit_simple, CrI = TRUE) ``` For the plotting the credible intervals are shown as well. From 3088a10536dcd25e5448c985fd57c366375904ec Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Thu, 12 Oct 2023 17:16:43 +0200 Subject: [PATCH 27/43] - adapted vignette to new plot function --- vignettes/analysis_normal.Rmd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 326337c..c190d1f 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -170,7 +170,7 @@ fit<- getModelFits( ``` ```{r} -plot.modelFits(fit, CrI = TRUE) -plot.modelFits(fit_simple, CrI = TRUE) +plot_modelFits(fit, cr_intv = TRUE) +plot_modelFits(fit_simple, cr_intv = TRUE, cr_bands = TRUE) ``` For the plotting the credible intervals are shown as well. From c375ffeff06c386d84ad7cf2c951d935eca7d6d1 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Thu, 12 Oct 2023 17:22:34 +0200 Subject: [PATCH 28/43] - fixed changes in plot test file due to re-naming --- tests/testthat/test-plot.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R index a6f37c0..d2ab7a5 100644 --- a/tests/testthat/test-plot.R +++ b/tests/testthat/test-plot.R @@ -153,8 +153,8 @@ test_that("Test plot_modelFits with different model_fits input", { plot6 <- plot_modelFits(fit) expect_is(plot6, "ggplot") - # Test with CrI = TRUE and more models - plot7 <- plot_modelFits(fit, CrI = TRUE) + # Test with cr_intv = TRUE and more models + plot7 <- plot_modelFits(fit, cr_intv = TRUE) expect_is(plot7, "ggplot") # Test with gAIC = FALSE and more models @@ -166,6 +166,6 @@ test_that("Test plot_modelFits with different model_fits input", { expect_is(plot9, "ggplot") # Test with all non-default parameters and more models - plot10 <- plot_modelFits(fit, CrI = TRUE, gAIC = FALSE, avg_fit = FALSE) + plot10 <- plot_modelFits(fit, cr_intv = TRUE, gAIC = FALSE, avg_fit = FALSE) expect_is(plot10, "ggplot") }) \ No newline at end of file From ccd270b5bc15f6d20bbebd5b23ddf8f2851eab6e Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Thu, 12 Oct 2023 17:31:28 +0200 Subject: [PATCH 29/43] - devtools::check() pleasing --- R/bootstrapping.R | 4 ++-- R/posterior.R | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/bootstrapping.R b/R/bootstrapping.R index 99eb507..3f7faf1 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -45,7 +45,7 @@ getBootsrapBands <- function ( bnds = DoseFinding::defBnds( mD = max(model_fits[[1]]$dose_levels))[[model]]) - preds <- predict(fit, doseSeq = dose_seq, predType = "ls-means") + preds <- stats::predict(fit, doseSeq = dose_seq, predType = "ls-means") attr(preds, "gAIC") <- DoseFinding::gAIC(fit) return (preds) @@ -74,7 +74,7 @@ getBootsrapBands <- function ( sort_indx <- order(rep(seq_along(model_names), length(dose_seq))) quant_mat <- t(apply(X = preds[sort_indx, ], MARGIN = 1, - FUN = quantile, + FUN = stats::quantile, probs = quantile_probs)) cr_bounds_data <- cbind( diff --git a/R/posterior.R b/R/posterior.R index 0c85d86..7b673f9 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -40,7 +40,7 @@ getPosteriorI <- function( if (is.null(mu_hat) && is.null(sd_hat)) { - anova_res <- lm(data_i$response ~ factor(data_i$dose) - 1) + 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] From 89fd571d38903fd484a49dcea3b91990d011c7cb Mon Sep 17 00:00:00 2001 From: Andersen Date: Fri, 13 Oct 2023 10:36:54 +0200 Subject: [PATCH 30/43] adjusted bootstrapping.R line, error in Vignette as well as created tests for bootstrapping file --- DESCRIPTION | 2 +- R/bootstrapping.R | 4 +- tests/testthat/test-bootstrapping.R | 152 ++++++++++++++++++++++++++++ tests/testthat/test-modelling.R | 2 +- tests/testthat/test-plot.R | 20 ++-- tests/testthat/test-posterior.R | 8 +- vignettes/analysis_normal.Rmd | 2 +- 7 files changed, 171 insertions(+), 19 deletions(-) create mode 100644 tests/testthat/test-bootstrapping.R diff --git a/DESCRIPTION b/DESCRIPTION index 06e6dd6..c1cc9bd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BayesianMCPMod Title: Bayesian MCPMod -Version: 0.1.3 +Version: 0.1.3-1 Authors@R: c( person("Sebastian", "Bossert", role = "aut", diff --git a/R/bootstrapping.R b/R/bootstrapping.R index 3f7faf1..b73bd4c 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -9,7 +9,6 @@ #' @return tbd #' @export getBootsrapBands <- function ( - model_fits, n_samples = 1e3, alpha = c(0.05, 0.5), @@ -17,10 +16,9 @@ getBootsrapBands <- function ( dose_seq = NULL ) { - mu_hat_samples <- sapply(attr(model_fits, "posterior"), RBesT::rmix, n = n_samples) - sd_hat <- summary.postList(post_observed)[, 2] + sd_hat <- summary.postList(model_fits)[, 2] dose_levels <- model_fits[[1]]$dose_levels model_names <- names(model_fits) diff --git a/tests/testthat/test-bootstrapping.R b/tests/testthat/test-bootstrapping.R new file mode 100644 index 0000000..00b5bdc --- /dev/null +++ b/tests/testthat/test-bootstrapping.R @@ -0,0 +1,152 @@ +# # Test cases +# test_that("test getBootsrapBands", { +# 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 <- BMCPMod( +# 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 <- getBootsrapBands(fit_simple) +# result <- getBootsrapBands(fit) +# expect_type(result_simple, "data.frame") +# expect_type(result, "data.frame") +# +# result_2_simple <- getBootsrapBands(fit_simple, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3)) +# result_2 <- getBootsrapBands(fit, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3)) +# expect_type(result_2_simple, "data.frame") +# expect_type(result_2, "data.frame") +# +# result_3_simple <- getBootsrapBands(fit_simple, dose_seq = NULL) +# result_3 <- getBootsrapBands(fit, dose_seq = NULL) +# expect_type(result_3_simple, "data.frame") +# expect_type(result_3, "data.frame") +# }) +# diff --git a/tests/testthat/test-modelling.R b/tests/testthat/test-modelling.R index 2de1ba4..e398d28 100644 --- a/tests/testthat/test-modelling.R +++ b/tests/testthat/test-modelling.R @@ -23,5 +23,5 @@ test_that("predictModelFit works correctly", { ) pred_vals <- predictModelFit(model_fit) - expect_is(pred_vals, "numeric") + expect_type(pred_vals, "double") }) diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R index d2ab7a5..b1c1381 100644 --- a/tests/testthat/test-plot.R +++ b/tests/testthat/test-plot.R @@ -150,22 +150,22 @@ test_that("Test plot_modelFits with different model_fits input", { ) # Test with default parameters and more models - plot6 <- plot_modelFits(fit) - expect_is(plot6, "ggplot") + plot1 <- plot_modelFits(fit) + expect_s3_class(plot1, "ggplot") # Test with cr_intv = TRUE and more models - plot7 <- plot_modelFits(fit, cr_intv = TRUE) - expect_is(plot7, "ggplot") + plot2 <- plot_modelFits(fit, cr_intv = TRUE) + expect_s3_class(plot2, "ggplot") # Test with gAIC = FALSE and more models - plot8 <- plot_modelFits(fit, gAIC = FALSE) - expect_is(plot8, "ggplot") + plot3 <- plot_modelFits(fit, gAIC = FALSE) + expect_s3_class(plot3, "ggplot") # Test with avg_fit = FALSE and more models - plot9 <- plot_modelFits(fit, avg_fit = FALSE) - expect_is(plot9, "ggplot") + plot4 <- plot_modelFits(fit, avg_fit = FALSE) + expect_s3_class(plot4, "ggplot") # Test with all non-default parameters and more models - plot10 <- plot_modelFits(fit, cr_intv = TRUE, gAIC = FALSE, avg_fit = FALSE) - expect_is(plot10, "ggplot") + 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 index e089f9a..47eeb69 100644 --- a/tests/testthat/test-posterior.R +++ b/tests/testthat/test-posterior.R @@ -9,7 +9,8 @@ test_that("getPosterior works correctly", { # Test getPosterior function posterior_list <- getPosterior(data, prior_list, mu_hat, sd_hat) - expect_is(posterior_list, "postList") + expect_type(posterior_list, "character") + expect_s3_class(posterior_list, "postList") }) test_that("getPosteriorI works correctly", { @@ -22,7 +23,8 @@ test_that("getPosteriorI works correctly", { # Test getPosteriorI function post_list <- getPosteriorI(data_i, prior_list, mu_hat, sd_hat) - expect_is(post_list, "postList") + expect_type(post_list, "character") + expect_s3_class(post_list, "postList") }) test_that("summary.postList works correctly", { @@ -37,5 +39,5 @@ test_that("summary.postList works correctly", { # Test summary.postList function summary_tab <- summary.postList(post_list) - expect_is(summary_tab, "matrix") + expect_type(summary_tab, "character") }) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index c190d1f..424448c 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -171,6 +171,6 @@ fit<- getModelFits( ```{r} plot_modelFits(fit, cr_intv = TRUE) -plot_modelFits(fit_simple, cr_intv = TRUE, cr_bands = TRUE) +# plot_modelFits(fit_simple, cr_intv = TRUE, cr_bands = TRUE) ``` For the plotting the credible intervals are shown as well. From f5f49132a0a4ebb55c86437cb0a7aa7eddebc65c Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Fri, 13 Oct 2023 10:45:16 +0200 Subject: [PATCH 31/43] - typo correction in bootstrapping function --- R/bootstrapping.R | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/R/bootstrapping.R b/R/bootstrapping.R index b73bd4c..db97781 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -9,6 +9,7 @@ #' @return tbd #' @export getBootsrapBands <- function ( + model_fits, n_samples = 1e3, alpha = c(0.05, 0.5), @@ -16,9 +17,10 @@ getBootsrapBands <- function ( dose_seq = NULL ) { + mu_hat_samples <- sapply(attr(model_fits, "posterior"), RBesT::rmix, n = n_samples) - sd_hat <- summary.postList(model_fits)[, 2] + sd_hat <- summary.postList(attr(model_fits, "posterior"))[, 2] dose_levels <- model_fits[[1]]$dose_levels model_names <- names(model_fits) From cedb075bff8f40e47fe8d8b343f58e4082ba683d Mon Sep 17 00:00:00 2001 From: Andersen Date: Fri, 13 Oct 2023 11:08:43 +0200 Subject: [PATCH 32/43] adjusted tests --- tests/testthat/test-bootstrapping.R | 304 ++++++++++++++-------------- vignettes/analysis_normal.Rmd | 2 +- 2 files changed, 153 insertions(+), 153 deletions(-) diff --git a/tests/testthat/test-bootstrapping.R b/tests/testthat/test-bootstrapping.R index 00b5bdc..0addcb5 100644 --- a/tests/testthat/test-bootstrapping.R +++ b/tests/testthat/test-bootstrapping.R @@ -1,152 +1,152 @@ -# # Test cases -# test_that("test getBootsrapBands", { -# 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 <- BMCPMod( -# 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 <- getBootsrapBands(fit_simple) -# result <- getBootsrapBands(fit) -# expect_type(result_simple, "data.frame") -# expect_type(result, "data.frame") -# -# result_2_simple <- getBootsrapBands(fit_simple, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3)) -# result_2 <- getBootsrapBands(fit, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3)) -# expect_type(result_2_simple, "data.frame") -# expect_type(result_2, "data.frame") -# -# result_3_simple <- getBootsrapBands(fit_simple, dose_seq = NULL) -# result_3 <- getBootsrapBands(fit, dose_seq = NULL) -# expect_type(result_3_simple, "data.frame") -# expect_type(result_3, "data.frame") -# }) -# +# Test cases +test_that("test getBootsrapBands", { + 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 <- BMCPMod( + 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 <- getBootsrapBands(fit_simple) + result <- getBootsrapBands(fit) + expect_type(result_simple, "list") + expect_type(result, "list") + + result_2_simple <- getBootsrapBands(fit_simple, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3)) + result_2 <- getBootsrapBands(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 <- getBootsrapBands(fit_simple, dose_seq = NULL) + result_3 <- getBootsrapBands(fit, dose_seq = NULL) + expect_type(result_3_simple, "list") + expect_type(result_3, "list") +}) + diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 424448c..c190d1f 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -171,6 +171,6 @@ fit<- getModelFits( ```{r} plot_modelFits(fit, cr_intv = TRUE) -# plot_modelFits(fit_simple, cr_intv = TRUE, cr_bands = TRUE) +plot_modelFits(fit_simple, cr_intv = TRUE, cr_bands = TRUE) ``` For the plotting the credible intervals are shown as well. From acd3f82b1c88f80a09914c2f0f1e678866acdfbf Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Fri, 13 Oct 2023 14:27:56 +0200 Subject: [PATCH 33/43] - added color for bootstrap credible bands (and median) in plot_modelFits() - changed cr_intv = TRUE as default in plot_modelFits() - added bootstrap median - added predict.ModelFits() - removed superfluous arguments from getModelFitSimple() - corrected typos --- NAMESPACE | 2 +- R/bootstrapping.R | 6 +- R/modelling.R | 31 ++++---- R/plot.R | 79 +++++++++++-------- ...tBootsrapBands.Rd => getBootstrapBands.Rd} | 10 +-- man/plot_modelFits.Rd | 7 +- vignettes/analysis_normal.Rmd | 5 +- 7 files changed, 76 insertions(+), 64 deletions(-) rename man/{getBootsrapBands.Rd => getBootstrapBands.Rd} (75%) diff --git a/NAMESPACE b/NAMESPACE index 20e67aa..97c3aee 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,7 +1,7 @@ # Generated by roxygen2: do not edit by hand export(BMCPMod) -export(getBootsrapBands) +export(getBootstrapBands) export(getModelFits) export(getPosterior) export(plot_modelFits) diff --git a/R/bootstrapping.R b/R/bootstrapping.R index db97781..8585b83 100644 --- a/R/bootstrapping.R +++ b/R/bootstrapping.R @@ -1,4 +1,4 @@ -#' @title getBootsrapBands +#' @title getBootstrapBands #' #' @param model_fits tbd #' @param n_samples tbd @@ -8,7 +8,7 @@ #' #' @return tbd #' @export -getBootsrapBands <- function ( +getBootstrapBands <- function ( model_fits, n_samples = 1e3, @@ -24,7 +24,7 @@ getBootsrapBands <- function ( dose_levels <- model_fits[[1]]$dose_levels model_names <- names(model_fits) - quantile_probs <- sort(unique(c(alpha / 2, 1 - alpha / 2))) + quantile_probs <- c(0.5, sort(unique(c(alpha / 2, 1 - alpha / 2)))) if (is.null(dose_seq)) { diff --git a/R/modelling.R b/R/modelling.R index fcad488..e070986 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -33,28 +33,14 @@ getModelFitSimple <- function ( model, dose_levels, - posterior, - mu_hat = NULL, - sd_hat = NULL + posterior ) { - if (is.null(mu_hat) && is.null(sd_hat)) { - - mu_hat <- summary.postList(posterior)[, 1] - sd_hat <- summary.postList(posterior)[, 2] - - } else { - - stopifnot(length(mu_hat) == length(sd_hat), - length(mu_hat) == length(dose_levels)) - - } - fit <- DoseFinding::fitMod( dose = dose_levels, - resp = mu_hat, - S = diag(sd_hat^2), + resp = summary.postList(posterior)[, 1], + S = diag(summary.postList(posterior)[, 2]^2), model = model, type = "general", bnds = DoseFinding::defBnds(mD = max(dose_levels))[[model]]) @@ -182,6 +168,17 @@ getModelFitOpt <- function ( } +predict.ModelFits <- function ( + + model_fits, + doses = NULL + +) { + + lapply(model_fits, predictModelFit, doses = doses) + +} + predictModelFit <- function ( model_fit, diff --git a/R/plot.R b/R/plot.R index 6fcaeab..0e07abd 100644 --- a/R/plot.R +++ b/R/plot.R @@ -8,7 +8,8 @@ #' @param cr_bands tbd #' @param alpha_CrB tbd #' @param n_bs_smpl tbd -#' +#' @param acc_color tbd +#' #' @return tbd #' @export plot_modelFits <- function ( @@ -16,15 +17,16 @@ plot_modelFits <- function ( model_fits, gAIC = TRUE, avg_fit = TRUE, - cr_intv = FALSE, + cr_intv = TRUE, alpha_CrI = 0.05, cr_bands = FALSE, alpha_CrB = c(0.05, 0.5), - n_bs_smpl = 1e3 + n_bs_smpl = 1e3, + acc_color = "orange" ) { - plot_resolution <- 1e2 + plot_res <- 1e2 dose_levels <- model_fits[[1]]$dose_levels post_summary <- summary.postList( @@ -32,7 +34,7 @@ plot_modelFits <- function ( probs = c(alpha_CrI / 2, 0.5, 1 - alpha_CrI / 2)) doses <- seq(from = min(dose_levels), to = max(dose_levels), - length.out = plot_resolution) + length.out = plot_res) preds_models <- sapply(model_fits, predictModelFit, doses = doses) model_names <- names(model_fits) @@ -54,7 +56,7 @@ plot_modelFits <- function ( levels = c("linear", "emax", "exponential", "sigEmax", "logistic", "quadratic", "avgFit")), - each = plot_resolution)) + each = plot_res)) if (gAIC) { @@ -81,7 +83,7 @@ plot_modelFits <- function ( if (cr_bands) { - crB_data <- getBootsrapBands( + crB_data <- getBootstrapBands( model_fits = model_fits, n_samples = n_bs_smpl, alpha = alpha_CrB, @@ -90,9 +92,8 @@ plot_modelFits <- function ( getInx <- function (alpha_CrB) { n <- length(alpha_CrB) - inx_list <- lapply(seq_len(n), function (i) c(i, 2 * n - i + 1) + 2) - return (inx_list) - } + inx_list <- lapply(seq_len(n), function (i) c(i, 2 * n - i + 1) + 3) + return (inx_list)} } @@ -115,7 +116,37 @@ plot_modelFits <- function ( 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) { @@ -133,8 +164,9 @@ plot_modelFits <- function ( } + ## Posterior Medians ggplot2::geom_point( - data = data.frame(dose_levels = dose_levels, - medians = post_summary[, 4]), + data = data.frame( + dose_levels = dose_levels, + medians = post_summary[, 4]), mapping = ggplot2::aes(dose_levels, medians), size = 2) + ## Fitted Models @@ -144,27 +176,6 @@ plot_modelFits <- function ( ## Faceting ggplot2::facet_wrap(~ models) - ## Bootstrapped Credible Bands - if (cr_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], "]), - alpha = 0.2)") - - plts <- plts + eval(parse(text = loop_txt)) - - } - rm(inx) - - } - return (plts) } \ No newline at end of file diff --git a/man/getBootsrapBands.Rd b/man/getBootstrapBands.Rd similarity index 75% rename from man/getBootsrapBands.Rd rename to man/getBootstrapBands.Rd index d036274..de4bd80 100644 --- a/man/getBootsrapBands.Rd +++ b/man/getBootstrapBands.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/bootstrapping.R -\name{getBootsrapBands} -\alias{getBootsrapBands} -\title{getBootsrapBands} +\name{getBootstrapBands} +\alias{getBootstrapBands} +\title{getBootstrapBands} \usage{ -getBootsrapBands( +getBootstrapBands( model_fits, n_samples = 1000, alpha = c(0.05, 0.5), @@ -27,5 +27,5 @@ getBootsrapBands( tbd } \description{ -getBootsrapBands +getBootstrapBands } diff --git a/man/plot_modelFits.Rd b/man/plot_modelFits.Rd index 3f80bc5..9a5f6b8 100644 --- a/man/plot_modelFits.Rd +++ b/man/plot_modelFits.Rd @@ -8,11 +8,12 @@ plot_modelFits( model_fits, gAIC = TRUE, avg_fit = TRUE, - cr_intv = FALSE, + cr_intv = TRUE, alpha_CrI = 0.05, cr_bands = FALSE, alpha_CrB = c(0.05, 0.5), - n_bs_smpl = 1000 + n_bs_smpl = 1000, + acc_color = "orange" ) } \arguments{ @@ -31,6 +32,8 @@ plot_modelFits( \item{alpha_CrB}{tbd} \item{n_bs_smpl}{tbd} + +\item{acc_color}{tbd} } \value{ tbd diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 424448c..0d1329c 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -18,6 +18,7 @@ knitr::opts_chunk$set( library(BayesianMCPMod) library(clinDR) library(dplyr) +set.seed(7015) ``` This vignette shows how to apply the BayesianMCPMod package. @@ -120,11 +121,11 @@ posterior_emax <- getPosterior( # iv) Execution of Bayesian MCPMod Test step ```{r} #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) @@ -171,6 +172,6 @@ fit<- getModelFits( ```{r} plot_modelFits(fit, cr_intv = TRUE) -# plot_modelFits(fit_simple, cr_intv = TRUE, cr_bands = TRUE) +plot_modelFits(fit_simple, cr_intv = TRUE, cr_bands = TRUE) ``` For the plotting the credible intervals are shown as well. From 9eba0fc3a97d0e80db9c1b46146eb62386ad0332 Mon Sep 17 00:00:00 2001 From: Andersen Date: Fri, 13 Oct 2023 14:46:20 +0200 Subject: [PATCH 34/43] adjusted test calls --- DESCRIPTION | 2 +- tests/testthat/test-bootstrapping.R | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index c1cc9bd..8744241 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BayesianMCPMod Title: Bayesian MCPMod -Version: 0.1.3-1 +Version: 0.1.3-2 Authors@R: c( person("Sebastian", "Bossert", role = "aut", diff --git a/tests/testthat/test-bootstrapping.R b/tests/testthat/test-bootstrapping.R index 0addcb5..c9e57a7 100644 --- a/tests/testthat/test-bootstrapping.R +++ b/tests/testthat/test-bootstrapping.R @@ -134,18 +134,18 @@ posterior = post_observed, posterior = post_observed, simple = FALSE) - result_simple <- getBootsrapBands(fit_simple) - result <- getBootsrapBands(fit) + result_simple <- getBootstrapBands(fit_simple) + result <- getBootstrapBands(fit) expect_type(result_simple, "list") expect_type(result, "list") - result_2_simple <- getBootsrapBands(fit_simple, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3)) - result_2 <- getBootsrapBands(fit, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3)) + 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 <- getBootsrapBands(fit_simple, dose_seq = NULL) - result_3 <- getBootsrapBands(fit, dose_seq = NULL) + 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") }) From aa1f665594dd3a29722e8ba0f0f2330e6becf8f8 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Fri, 13 Oct 2023 15:07:23 +0200 Subject: [PATCH 35/43] - corrected typo on bootstrap tests - updated plot function in vignette --- tests/testthat/test-bootstrapping.R | 14 +++++++------- vignettes/analysis_normal.Rmd | 4 ++-- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/testthat/test-bootstrapping.R b/tests/testthat/test-bootstrapping.R index 0addcb5..0772093 100644 --- a/tests/testthat/test-bootstrapping.R +++ b/tests/testthat/test-bootstrapping.R @@ -1,5 +1,5 @@ # Test cases -test_that("test getBootsrapBands", { +test_that("test getBootstrapBands", { library(clinDR) library(dplyr) @@ -134,18 +134,18 @@ posterior = post_observed, posterior = post_observed, simple = FALSE) - result_simple <- getBootsrapBands(fit_simple) - result <- getBootsrapBands(fit) + result_simple <- getBootstrapBands(fit_simple) + result <- getBootstrapBands(fit) expect_type(result_simple, "list") expect_type(result, "list") - result_2_simple <- getBootsrapBands(fit_simple, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3)) - result_2 <- getBootsrapBands(fit, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3)) + 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 <- getBootsrapBands(fit_simple, dose_seq = NULL) - result_3 <- getBootsrapBands(fit, dose_seq = NULL) + 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/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 0d1329c..2609346 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -171,7 +171,7 @@ fit<- getModelFits( ``` ```{r} -plot_modelFits(fit, cr_intv = TRUE) -plot_modelFits(fit_simple, cr_intv = TRUE, cr_bands = TRUE) +plot_modelFits(fit) +plot_modelFits(fit_simple, cr_bands = TRUE) ``` For the plotting the credible intervals are shown as well. From eddb11152d1dc89b00481f05baa593e757fffcff Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Fri, 13 Oct 2023 16:03:35 +0200 Subject: [PATCH 36/43] - introduced S3 methods for modelFits and postList: plot and predict for the former, print and summary for the latter --- NAMESPACE | 5 ++++- R/modelling.R | 8 +++++--- R/plot.R | 19 +++++++++++-------- R/posterior.R | 15 ++++++++++----- man/{plot_modelFits.Rd => plot.modelFits.Rd} | 19 +++++++++++-------- tests/testthat/test-plot.R | 12 ++++++------ vignettes/analysis_normal.Rmd | 4 ++-- 7 files changed, 49 insertions(+), 33 deletions(-) rename man/{plot_modelFits.Rd => plot.modelFits.Rd} (69%) diff --git a/NAMESPACE b/NAMESPACE index 97c3aee..3cd5f34 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,8 +1,11 @@ # Generated by roxygen2: do not edit by hand +S3method(plot,modelFits) +S3method(predict,ModelFits) +S3method(print,postList) +S3method(summary,postList) export(BMCPMod) export(getBootstrapBands) export(getModelFits) export(getPosterior) -export(plot_modelFits) export(simulateData) diff --git a/R/modelling.R b/R/modelling.R index e070986..d4b1c04 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -168,14 +168,16 @@ getModelFitOpt <- function ( } +#' @export predict.ModelFits <- function ( - model_fits, - doses = NULL + object, + doses = NULL, + ... ) { - lapply(model_fits, predictModelFit, doses = doses) + lapply(x, predictModelFit, doses = doses, ...) } diff --git a/R/plot.R b/R/plot.R index 0e07abd..9f8278b 100644 --- a/R/plot.R +++ b/R/plot.R @@ -1,6 +1,6 @@ -#' @title plot_modelFits +#' @title plot.modelFits #' -#' @param model_fits tbd +#' @param x tbd An object of type modelFits #' @param gAIC tbd #' @param avg_fit tbd #' @param cr_intv tbd @@ -9,12 +9,13 @@ #' @param alpha_CrB tbd #' @param n_bs_smpl tbd #' @param acc_color tbd +#' @param ... tbd #' #' @return tbd #' @export -plot_modelFits <- function ( +plot.modelFits <- function ( - model_fits, + x, gAIC = TRUE, avg_fit = TRUE, cr_intv = TRUE, @@ -22,16 +23,18 @@ plot_modelFits <- function ( cr_bands = FALSE, alpha_CrB = c(0.05, 0.5), n_bs_smpl = 1e3, - acc_color = "orange" + acc_color = "orange", + ... ) { - plot_res <- 1e2 + plot_res <- 1e2 + model_fits <- x dose_levels <- model_fits[[1]]$dose_levels post_summary <- summary.postList( - post_list = attr(model_fits, "posterior"), - probs = c(alpha_CrI / 2, 0.5, 1 - alpha_CrI / 2)) + 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) diff --git a/R/posterior.R b/R/posterior.R index 7b673f9..fa36e08 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -84,27 +84,32 @@ getPostCombsI <- function ( } +#' @export summary.postList <- function ( - post_list, + object, ... ) { - summary_list <- lapply(post_list, summary, ...) - names(summary_list) <- names(post_list) + summary_list <- lapply(object, summary, ...) + names(summary_list) <- names(object) summary_tab <- do.call(rbind, summary_list) return (summary_tab) } +#' @export print.postList <- function ( - post_list + x, + ... ) { + post_list <- x + getMaxDiff <- function ( medians @@ -133,6 +138,6 @@ print.postList <- function ( "Maximum Difference to Control and Dose Group", "Posterior Distributions") - print(list_out) + print(list_out, ...) } diff --git a/man/plot_modelFits.Rd b/man/plot.modelFits.Rd similarity index 69% rename from man/plot_modelFits.Rd rename to man/plot.modelFits.Rd index 9a5f6b8..3b35078 100644 --- a/man/plot_modelFits.Rd +++ b/man/plot.modelFits.Rd @@ -1,11 +1,11 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/plot.R -\name{plot_modelFits} -\alias{plot_modelFits} -\title{plot_modelFits} +\name{plot.modelFits} +\alias{plot.modelFits} +\title{plot.modelFits} \usage{ -plot_modelFits( - model_fits, +\method{plot}{modelFits}( + x, gAIC = TRUE, avg_fit = TRUE, cr_intv = TRUE, @@ -13,11 +13,12 @@ plot_modelFits( cr_bands = FALSE, alpha_CrB = c(0.05, 0.5), n_bs_smpl = 1000, - acc_color = "orange" + acc_color = "orange", + ... ) } \arguments{ -\item{model_fits}{tbd} +\item{x}{tbd An object of type modelFits} \item{gAIC}{tbd} @@ -34,10 +35,12 @@ plot_modelFits( \item{n_bs_smpl}{tbd} \item{acc_color}{tbd} + +\item{...}{tbd} } \value{ tbd } \description{ -plot_modelFits +plot.modelFits } diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R index b1c1381..a41749c 100644 --- a/tests/testthat/test-plot.R +++ b/tests/testthat/test-plot.R @@ -1,4 +1,4 @@ -test_that("Test plot_modelFits with different model_fits input", { +test_that("Test plot.modelFits with different model_fits input", { library(BayesianMCPMod) library(clinDR) library(dplyr) @@ -150,22 +150,22 @@ test_that("Test plot_modelFits with different model_fits input", { ) # Test with default parameters and more models - plot1 <- plot_modelFits(fit) + plot1 <- plot.modelFits(fit) expect_s3_class(plot1, "ggplot") # Test with cr_intv = TRUE and more models - plot2 <- plot_modelFits(fit, cr_intv = TRUE) + 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) + 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) + 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) + 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/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 2609346..6f1c0e0 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -171,7 +171,7 @@ fit<- getModelFits( ``` ```{r} -plot_modelFits(fit) -plot_modelFits(fit_simple, cr_bands = TRUE) +plot(fit) +plot(fit_simple, cr_bands = TRUE) ``` For the plotting the credible intervals are shown as well. From 7f516af9512d4ab46b9002faae09dcd6272d5a7e Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Fri, 13 Oct 2023 16:54:17 +0200 Subject: [PATCH 37/43] - corrected typo in predict.ModelFit() --- R/modelling.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/modelling.R b/R/modelling.R index d4b1c04..d28d179 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -177,7 +177,7 @@ predict.ModelFits <- function ( ) { - lapply(x, predictModelFit, doses = doses, ...) + lapply(object, predictModelFit, doses = doses, ...) } From db33d70c0658798e7cfcd2094e12a9c0b3733c66 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Sat, 14 Oct 2023 13:37:18 +0200 Subject: [PATCH 38/43] - added S3 printing methods for BayesianMCP, BayesianMCPMod, modelFits - added separate file for S3 methods - added function performMayesianMCPMod - renamed function BMCPMod to perform BayesianMCP - renamed function BayesMCPMod to BayesMCPi - unified handling of posterior_list vs posterior in performBayesianMCP and getModelFits --- NAMESPACE | 6 +- R/BMCPMod.R | 109 +++++++++++-- R/modelling.R | 13 -- R/posterior.R | 56 ------- R/s3methods.R | 185 ++++++++++++++++++++++ man/{BMCPMod.Rd => performBayesianMCP.Rd} | 10 +- man/performBayesianMCPMod.Rd | 20 +++ tests/testthat/test-bootstrapping.R | 2 +- tests/testthat/test-plot.R | 2 +- vignettes/analysis_normal.Rmd | 4 +- 10 files changed, 319 insertions(+), 88 deletions(-) create mode 100644 R/s3methods.R rename man/{BMCPMod.Rd => performBayesianMCP.Rd} (58%) create mode 100644 man/performBayesianMCPMod.Rd diff --git a/NAMESPACE b/NAMESPACE index 3cd5f34..3f6d834 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,10 +2,14 @@ S3method(plot,modelFits) S3method(predict,ModelFits) +S3method(print,BayesianMCP) +S3method(print,BayesianMCPMod) +S3method(print,modelFits) S3method(print,postList) S3method(summary,postList) -export(BMCPMod) export(getBootstrapBands) export(getModelFits) export(getPosterior) +export(performBayesianMCP) +export(performBayesianMCPMod) export(simulateData) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index abdc170..dc20bbe 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -1,11 +1,92 @@ -#' @title BMCPMod +#' @title performBayesianMCPMod #' #' @param posteriors_list tbd #' @param contr_mat tbd #' @param crit_prob tbd +#' @param simple tbd #' #' @export -BMCPMod <- function( +performBayesianMCPMod <- function ( + + posteriors_list, + contr_mat, + crit_prob, + simple = FALSE + +) { + + if (class(posteriors_list) == "postList") { + + posteriors_list <- list(posteriors_list) + + } + + b_mcp <- BayesianMCP( + 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) + +} + +#' @title BayesianMCP +#' +#' @param posteriors_list tbd +#' @param contr_mat tbd +#' @param crit_prob tbd +#' +#' @export +performBayesianMCP <- function( posteriors_list, contr_mat, @@ -13,11 +94,22 @@ BMCPMod <- function( ) { - lapply(posteriors_list, BayesMCPMod, 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) } -BayesMCPMod <- function ( +BayesMCPi <- function ( posterior_i, contr_mat, @@ -51,11 +143,10 @@ BayesMCPMod <- function ( post_combs_i <- getPostCombsI(posterior_i) post_probs <- apply(contr_mat$contMat, 2, getPostProb, post_combs_i) - res_list <- list( - sign = ifelse(max(post_probs) > crit_prob, 1, 0), - p_val = max(post_probs), - post_probs = post_probs) + res <- c(sign = ifelse(max(post_probs) > crit_prob, 1, 0), + p_val = max(post_probs), + post_probs = post_probs) - return (res_list) + return (res) } diff --git a/R/modelling.R b/R/modelling.R index d28d179..7895095 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -168,19 +168,6 @@ getModelFitOpt <- function ( } -#' @export -predict.ModelFits <- function ( - - object, - doses = NULL, - ... - -) { - - lapply(object, predictModelFit, doses = doses, ...) - -} - predictModelFit <- function ( model_fit, diff --git a/R/posterior.R b/R/posterior.R index fa36e08..8a22a58 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -84,60 +84,4 @@ getPostCombsI <- 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, - ... - -) { - - post_list <- 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(post_list) - - names(post_list) <- rownames(summary_tab) - class(post_list) <- NULL - - list_out <- list(summary_tab, getMaxDiff(summary_tab[, 4]), post_list) - names(list_out) <- c("Summary of Posterior Distributions", - "Maximum Difference to Control and Dose Group", - "Posterior Distributions") - - print(list_out, ...) - -} diff --git a/R/s3methods.R b/R/s3methods.R new file mode 100644 index 0000000..2ec698f --- /dev/null +++ b/R/s3methods.R @@ -0,0 +1,185 @@ +## BayesianMCPMod ----------------------------------------- + +#' @export +print.BayesianMCPMod <- function ( + + x, + ... + +) { + + n_models <- ncol(x$BayesianMCP) - 2L + + model_success <- colMeans(do.call(rbind, lapply(x$Mod, function (y) { + + if (!is.null(y)) { + + sapply(y, function (z) z$significant) + + } else { + + rep(FALSE, n_models) + + } + + }))) + + 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/man/BMCPMod.Rd b/man/performBayesianMCP.Rd similarity index 58% rename from man/BMCPMod.Rd rename to man/performBayesianMCP.Rd index 874fc43..e6b3543 100644 --- a/man/BMCPMod.Rd +++ b/man/performBayesianMCP.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/BMCPMod.R -\name{BMCPMod} -\alias{BMCPMod} -\title{BMCPMod} +\name{performBayesianMCP} +\alias{performBayesianMCP} +\title{BayesianMCP} \usage{ -BMCPMod(posteriors_list, contr_mat, crit_prob) +performBayesianMCP(posteriors_list, contr_mat, crit_prob) } \arguments{ \item{posteriors_list}{tbd} @@ -14,5 +14,5 @@ BMCPMod(posteriors_list, contr_mat, crit_prob) \item{crit_prob}{tbd} } \description{ -BMCPMod +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/tests/testthat/test-bootstrapping.R b/tests/testthat/test-bootstrapping.R index 0772093..c4caf87 100644 --- a/tests/testthat/test-bootstrapping.R +++ b/tests/testthat/test-bootstrapping.R @@ -108,7 +108,7 @@ test_that("test getBootstrapBands", { doses = dose_levels, w = n_patients + ess_prior) - BMCP_result <- BMCPMod( + BMCP_result <- performBayesianMCP( posteriors_list = list(posterior_emax), contr_mat = contr_mat_prior, crit_prob = crit_pval) diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R index a41749c..732580d 100644 --- a/tests/testthat/test-plot.R +++ b/tests/testthat/test-plot.R @@ -121,7 +121,7 @@ test_that("Test plot.modelFits with different model_fits input", { w = n_patients + ess_prior ) - BMCP_result <- BMCPMod( + BMCP_result <- performBayesianMCP( posteriors_list = list(posterior_emax), contr_mat = contr_mat_prior, crit_prob = crit_pval diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 6f1c0e0..3a2ece6 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -138,8 +138,8 @@ contr_mat_prior <- DoseFinding::optContr( doses = dose_levels, w = n_patients + ess_prior) -BMCP_result <- BMCPMod( - posteriors_list = list(posterior_emax), +BMCP_result <- performBayesianMCP( + posteriors_list = posterior_emax, contr_mat = contr_mat_prior, crit_prob = crit_pval) From c332b7a4990685c4f3324878ed7c31a32028078d Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Sat, 14 Oct 2023 15:08:56 +0200 Subject: [PATCH 39/43] - added function assessDesign() - fixed typo in performBayesianMCPMod - added functionality to print.BayesianMCPMod - added getModelData - to vignette: added automatic numbering, and assessment of trial design --- NAMESPACE | 1 + R/BMCPMod.R | 74 ++++++++++++++++++++++++++++++++++- R/s3methods.R | 7 +++- R/simulation.R | 24 ++++++++++++ man/assessDesign.Rd | 37 ++++++++++++++++++ vignettes/analysis_normal.Rmd | 26 ++++++++---- 6 files changed, 160 insertions(+), 9 deletions(-) create mode 100644 man/assessDesign.Rd diff --git a/NAMESPACE b/NAMESPACE index 3f6d834..f78c639 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,6 +7,7 @@ S3method(print,BayesianMCPMod) S3method(print,modelFits) S3method(print,postList) S3method(summary,postList) +export(assessDesign) export(getBootstrapBands) export(getModelFits) export(getPosterior) diff --git a/R/BMCPMod.R b/R/BMCPMod.R index dc20bbe..375b1ba 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -1,3 +1,75 @@ +#' @title assessDesign +#' +#' @param n_patients tbd +#' @param dose_levels tbd +#' @param sd 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, + dose_levels, + sd, + mods, + prior_list, + + n_sim = 1e3, + alpha_crit_val = 0.05, + simple = TRUE + +) { + + data <- simulateData( + n_patients = n_patients, + dose_levels = dose_levels, + sd = sd, + 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) + + contr_mat <- DoseFinding::optContr( + models = mods, + doses = dose_levels, + w = n_patients) + + crit_pval <- pnorm(DoseFinding:::critVal( + corMat = contr_mat$corMat, + alpha = alpha_crit_val, + df = 0, + alternative = "one.sided")) + + ess_prior <- suppressMessages(round(unlist(lapply(prior_list, RBesT::ess)))) + contr_mat_prior <- DoseFinding::optContr( + models = mods, + doses = dose_levels, + w = n_patients + ess_prior) + + 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 performBayesianMCPMod #' #' @param posteriors_list tbd @@ -21,7 +93,7 @@ performBayesianMCPMod <- function ( } - b_mcp <- BayesianMCP( + b_mcp <- performBayesianMCP( posteriors_list = posteriors_list, contr_mat = contr_mat, crit_prob = crit_prob) diff --git a/R/s3methods.R b/R/s3methods.R index 2ec698f..14aa663 100644 --- a/R/s3methods.R +++ b/R/s3methods.R @@ -9,6 +9,8 @@ print.BayesianMCPMod <- function ( ) { 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) { @@ -18,7 +20,10 @@ print.BayesianMCPMod <- function ( } else { - rep(FALSE, n_models) + model_signs <- rep(FALSE, n_models) + names(model_signs) <- model_names + + return (model_signs) } diff --git a/R/simulation.R b/R/simulation.R index c8ee63e..ea38651 100644 --- a/R/simulation.R +++ b/R/simulation.R @@ -30,3 +30,27 @@ simulateData <- function( 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) + +} + +getTrialOutcome <- function ( + + + +) { + + + +} \ No newline at end of file diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd new file mode 100644 index 0000000..930d34f --- /dev/null +++ b/man/assessDesign.Rd @@ -0,0 +1,37 @@ +% 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, + dose_levels, + sd, + mods, + prior_list, + n_sim = 1000, + alpha_crit_val = 0.05, + simple = TRUE +) +} +\arguments{ +\item{n_patients}{tbd} + +\item{dose_levels}{tbd} + +\item{sd}{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/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index 3a2ece6..d808aab 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -1,6 +1,7 @@ --- 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} @@ -18,11 +19,12 @@ knitr::opts_chunk$set( library(BayesianMCPMod) library(clinDR) library(dplyr) + set.seed(7015) ``` This vignette shows how to apply the BayesianMCPMod package. -# i) Calculation of a MAP prior +# Calculation of a MAP prior In a first step a meta analytic prior will be calculated using historical data. ```{r} data("metaData") @@ -68,7 +70,7 @@ 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]))) ``` -# ii) Pre-Specification of candidate models +# Pre-Specification of candidate models ```{r} #Pre-Specification (B)MCPMod @@ -82,7 +84,6 @@ emax <- DoseFinding::guesst(d = 100, p = c(0.9), model = "emax") - mods <- DoseFinding::Mods( linear = NULL, emax = emax, @@ -91,7 +92,19 @@ mods <- DoseFinding::Mods( maxEff = 10, placEff = 1.4) ``` -# iii) Trial results +# Assessing the Trial Design +```{r} +n_patients <- c(50, 50, 50, 50) + +assessDesign( + n_patients = n_patients, + dose_levels = dose_levels, + sd = sd_tot, + mods = mods, + prior_list = prior_list, + n_sim = 1e1) +``` +# Trial results ```{r} #Simulation of new trial ##Note: This part will be simplified and direct results from one trial will be used @@ -101,7 +114,6 @@ mods_sim <- DoseFinding::Mods( maxEff = 12, placEff = 1.4) -n_patients <- c(50, 50, 50, 50) data <- simulateData( n_patients = n_patients, dose_levels = dose_levels, @@ -118,7 +130,7 @@ posterior_emax <- getPosterior( ``` -# iv) Execution of Bayesian MCPMod Test step +# Execution of Bayesian MCPMod Test step ```{r} #Evaluation of Bayesian MCPMod contr_mat <- DoseFinding::optContr( @@ -146,7 +158,7 @@ BMCP_result <- performBayesianMCP( BMCP_result ``` -# v) Model fitting and Plotting +# Model fitting and Plotting In the model fitting step the posterior distribution is used as basis. ```{r} #Model fit From a1ed96df37c800c677eb7ce6b9969d2167d79116 Mon Sep 17 00:00:00 2001 From: Wojciekowski Date: Sat, 14 Oct 2023 17:22:22 +0200 Subject: [PATCH 40/43] - added new functions getContrMat, getCritProb, getPriorList - changed simulation to output a single data set if !is.null(true_model) - updated vignette --- NAMESPACE | 3 + R/BMCPMod.R | 93 +++++++++++++++++++++------ R/posterior.R | 71 ++++++++++++++++++++- R/simulation.R | 30 +++++---- man/assessDesign.Rd | 6 -- man/getContrMat.Rd | 20 ++++++ man/getCritProb.Rd | 20 ++++++ man/getPriorList.Rd | 20 ++++++ man/simulateData.Rd | 11 +++- vignettes/analysis_normal.Rmd | 117 +++++++++++----------------------- 10 files changed, 273 insertions(+), 118 deletions(-) create mode 100644 man/getContrMat.Rd create mode 100644 man/getCritProb.Rd create mode 100644 man/getPriorList.Rd diff --git a/NAMESPACE b/NAMESPACE index f78c639..f89ff4d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,8 +9,11 @@ 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 375b1ba..0a3f56c 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -1,8 +1,6 @@ #' @title assessDesign #' #' @param n_patients tbd -#' @param dose_levels tbd -#' @param sd tbd #' @param mods tbd #' @param prior_list tbd #' @param n_sim tbd @@ -13,8 +11,6 @@ assessDesign <- function ( n_patients, - dose_levels, - sd, mods, prior_list, @@ -24,10 +20,12 @@ assessDesign <- function ( ) { + dose_levels <- attr(prior_list, "dose_levels") + data <- simulateData( n_patients = n_patients, dose_levels = dose_levels, - sd = sd, + sd = attr(prior_list, "sd_tot"), mods = mods, n_sim = n_sim) @@ -39,22 +37,17 @@ assessDesign <- function ( data = getModelData(data, model_name), prior_list = prior_list) - contr_mat <- DoseFinding::optContr( - models = mods, - doses = dose_levels, - w = n_patients) + crit_pval <- getCritProb( + mods = mods, + dose_levels = dose_levels, + dose_weights = n_patients, + alpha_crit_val = alpha_crit_val) - crit_pval <- pnorm(DoseFinding:::critVal( - corMat = contr_mat$corMat, - alpha = alpha_crit_val, - df = 0, - alternative = "one.sided")) - - ess_prior <- suppressMessages(round(unlist(lapply(prior_list, RBesT::ess)))) - contr_mat_prior <- DoseFinding::optContr( - models = mods, - doses = dose_levels, - w = n_patients + ess_prior) + 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, @@ -70,6 +63,66 @@ assessDesign <- function ( } +#' @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 diff --git a/R/posterior.R b/R/posterior.R index 8a22a58..df37f5e 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -1,3 +1,66 @@ +#' @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 @@ -59,7 +122,13 @@ getPosteriorI <- function( post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = sd_hat) - names(post_list) <- c("Ctr", paste0("DG_", seq_along(post_list[-1]))) + 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) diff --git a/R/simulation.R b/R/simulation.R index ea38651..7c60082 100644 --- a/R/simulation.R +++ b/R/simulation.R @@ -5,6 +5,7 @@ #' @param sd tbd #' @param mods tbd #' @param n_sim tbd +#' @param true_model tbd #' #' @export simulateData <- function( @@ -13,10 +14,21 @@ simulateData <- function( dose_levels, sd, mods, - n_sim + 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), @@ -27,6 +39,12 @@ simulateData <- function( sim_data <- cbind(sim_info, model_responses + random_noise) + if (!is.null(true_model)) { + + sim_data <- getModelData(sim_data, true_model) + + } + return (sim_data) } @@ -43,14 +61,4 @@ getModelData <- function ( return (model_data) -} - -getTrialOutcome <- function ( - - - -) { - - - } \ No newline at end of file diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd index 930d34f..6883d23 100644 --- a/man/assessDesign.Rd +++ b/man/assessDesign.Rd @@ -6,8 +6,6 @@ \usage{ assessDesign( n_patients, - dose_levels, - sd, mods, prior_list, n_sim = 1000, @@ -18,10 +16,6 @@ assessDesign( \arguments{ \item{n_patients}{tbd} -\item{dose_levels}{tbd} - -\item{sd}{tbd} - \item{mods}{tbd} \item{prior_list}{tbd} 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/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/simulateData.Rd b/man/simulateData.Rd index b3572d5..32038df 100644 --- a/man/simulateData.Rd +++ b/man/simulateData.Rd @@ -4,7 +4,14 @@ \alias{simulateData} \title{simulateData} \usage{ -simulateData(n_patients, dose_levels, sd, mods, n_sim) +simulateData( + n_patients, + dose_levels, + sd, + mods, + n_sim = 1000, + true_model = NULL +) } \arguments{ \item{n_patients}{tbd} @@ -16,6 +23,8 @@ simulateData(n_patients, dose_levels, sd, mods, n_sim) \item{mods}{tbd} \item{n_sim}{tbd} + +\item{true_model}{tbd} } \description{ simulateData diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index d808aab..eabf560 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -26,7 +26,7 @@ This vignette shows how to apply the BayesianMCPMod package. # Calculation of a MAP prior In a first step a meta analytic prior will be calculated using historical data. -```{r} +```{r Calculation of a MAP prior} data("metaData") testdata <- as.data.frame(metaData) dataset <- filter(testdata, bname == "VIAGRA") @@ -34,55 +34,34 @@ histcontrol <- filter(dataset, dose == 0, primtime == 12, indication == "ERECTIL ##Create MAP Prior hist_data <- data.frame( - trial = c(1, 2, 3, 4), + trial = seq_len(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) +dose_levels <- c(0, 50, 100, 200) -#RBesT::ess(prior_ctr) +prior_list <- getPriorList( + hist_data = hist_data, + dose_levels = dose_levels) -## 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 of candidate models -```{r} +```{r Pre-Specification of candidate models} #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") +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, @@ -93,62 +72,44 @@ mods <- DoseFinding::Mods( placEff = 1.4) ``` # Assessing the Trial Design -```{r} +```{r Assessing the Trial Design} n_patients <- c(50, 50, 50, 50) assessDesign( n_patients = n_patients, - dose_levels = dose_levels, - sd = sd_tot, mods = mods, prior_list = prior_list, - n_sim = 1e1) + n_sim = 1e2) ``` # Trial results -```{r} +```{r Trial results} #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) - -data <- simulateData( +data_emax <- 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" + sd = attr(prior_list, "sd_tot"), + mods = mods, + true_model = "emax") posterior_emax <- getPosterior( data = data_emax, prior_list = prior_list) - ``` # Execution of Bayesian MCPMod Test step -```{r} -#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) +```{r Execution of Bayesian MCPMod Test step} +crit_pval <- getCritProb( + mods = mods, + dose_levels = dose_levels, + dose_weights = n_patients, + alpha_crit_val = 0.05) + +contr_mat_prior <- getContrMat( + mods = mods, + dose_levels = dose_levels, + dose_weights = n_patients, + prior_list = prior_list) BMCP_result <- performBayesianMCP( posteriors_list = posterior_emax, @@ -160,13 +121,11 @@ BMCP_result # Model fitting and Plotting In the model fitting step the posterior distribution is used as basis. -```{r} +```{r Model fitting} #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, @@ -175,14 +134,14 @@ fit_simple <- getModelFits( simple = TRUE) # Option b) Making use of the complete posterior distribution -fit<- getModelFits( +fit <- getModelFits( models = model_shapes, dose_levels = dose_levels, posterior = post_observed, simple = FALSE) ``` -```{r} +```{r Plotting} plot(fit) plot(fit_simple, cr_bands = TRUE) ``` From b45e0b29534c1bcc83022c5a08934c7d45fcd464 Mon Sep 17 00:00:00 2001 From: Bossert Date: Mon, 16 Oct 2023 14:23:55 +0200 Subject: [PATCH 41/43] Vignette update (other trial example, making use of newly available functions,...) --- vignettes/analysis_normal.Rmd | 115 +++++++++++++++++++++------------- 1 file changed, 73 insertions(+), 42 deletions(-) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index eabf560..1e8098c 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -22,69 +22,78 @@ library(dplyr) set.seed(7015) ``` -This vignette shows how to apply the BayesianMCPMod package. +# 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. +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. +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 == "VIAGRA") -histcontrol <- filter(dataset, dose == 0, primtime == 12, indication == "ERECTILE DYSFUNCTION") +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 = seq_len(4), + trial = histcontrol$nctno, est = histcontrol$rslt, se = histcontrol$se, sd = histcontrol$sd, n = histcontrol$sampsize) -dose_levels <- c(0, 50, 100, 200) +dose_levels <- c(0, 2.5, 5, 10) prior_list <- getPriorList( hist_data = hist_data, dose_levels = dose_levels) ``` -# Pre-Specification of candidate models + +# 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 = 50, - p = c(0.2), - model = "exponential", - Maxd = max(dose_levels)) -emax <- DoseFinding::guesst( - d = 100, - p = c(0.9), - model = "emax") +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 = 10, - placEff = 1.4) -``` -# Assessing the Trial Design -```{r Assessing the Trial Design} -n_patients <- c(50, 50, 50, 50) + maxEff = -1, + placEff = -12.8) + +new_trial <- filter(dataset, primtime == 8, indication == "MAJOR DEPRESSIVE DISORDER",protid==6) +n_patients <- c(150, 142, 147, 149) -assessDesign( - n_patients = n_patients, - mods = mods, - prior_list = prior_list, - n_sim = 1e2) ``` -# Trial results + +# 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 part will be simplified and direct results from one trial will be used +##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, @@ -92,18 +101,22 @@ data_emax <- simulateData( mods = mods, true_model = "emax") -posterior_emax <- getPosterior( - data = data_emax, - prior_list = prior_list) +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.05) + alpha_crit_val = 0.1) contr_mat_prior <- getContrMat( mods = mods, @@ -111,19 +124,35 @@ contr_mat_prior <- getContrMat( 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_emax, + posteriors_list = posterior, contr_mat = contr_mat_prior, - crit_prob = crit_pval) + 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. -# Model fitting and Plotting -In the model fitting step the posterior distribution is used as basis. ```{r Model fitting} #Model fit -post_observed <- posterior_emax +post_observed <- posterior model_shapes <- c("linear", "emax", "exponential") # Option a) Simplified approach by using frequentist idea @@ -139,10 +168,12 @@ fit <- getModelFits( dose_levels = dose_levels, posterior = post_observed, simple = FALSE) -``` -```{r Plotting} -plot(fit) +plot(fit,cr_bands=TRUE) plot(fit_simple, cr_bands = TRUE) + ``` -For the plotting the credible intervals are shown as well. +# Additional notes +TBD, whether certain wrapper functions should be mentioned. + + From 31ebd7fa6b6f27fb8127b7937d78ea07ede98783 Mon Sep 17 00:00:00 2001 From: Bossert Date: Wed, 18 Oct 2023 15:38:46 +0200 Subject: [PATCH 42/43] New draft vignette around simulation added --- _pkgdown.yml | 1 + vignettes/Simulation_Example.Rmd | 58 ++++++++++++++++++++++++++++++++ 2 files changed, 59 insertions(+) create mode 100644 vignettes/Simulation_Example.Rmd diff --git a/_pkgdown.yml b/_pkgdown.yml index 4cfd6f8..1b0d774 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -12,3 +12,4 @@ articles: navbar: ~ contents: - analysis_normal + - Simulation_Example diff --git a/vignettes/Simulation_Example.Rmd b/vignettes/Simulation_Example.Rmd new file mode 100644 index 0000000..ade8ccc --- /dev/null +++ b/vignettes/Simulation_Example.Rmd @@ -0,0 +1,58 @@ +--- +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 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") + +##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) + +``` + From 09b214b5ccdffd045993cbe131e85327c5cd3d8a Mon Sep 17 00:00:00 2001 From: Bossert Date: Thu, 19 Oct 2023 10:37:44 +0200 Subject: [PATCH 43/43] Further updates in vignettes. --- vignettes/Simulation_Example.Rmd | 59 +++++++++++++++++++++++++++++++- vignettes/analysis_normal.Rmd | 2 +- 2 files changed, 59 insertions(+), 2 deletions(-) diff --git a/vignettes/Simulation_Example.Rmd b/vignettes/Simulation_Example.Rmd index ade8ccc..bcd8cef 100644 --- a/vignettes/Simulation_Example.Rmd +++ b/vignettes/Simulation_Example.Rmd @@ -30,7 +30,7 @@ As in [link other vignette] we focus on the indication MDD and make use of histo 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 4 trials (with main endpoint CfB in MADRS score after 8 weeks). +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. @@ -56,3 +56,60 @@ prior_list <- getPriorList( ``` +# 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 index 1e8098c..e84ebac 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -25,7 +25,7 @@ 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. +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...