Nothing
R/BMCPMod.R
In BayesianMCPMod: Simulate, Evaluate, and Analyze Dose Finding Trials with Bayesian MCPMod
Defines functions
performBayesianMCP
Documented in
performBayesianMCP
addSignificance <- function (
model_fits,
sign_models
) {
sign_models_fits <- sapply(names(model_fits), function (model_name) {
if (model_name != "avgFit") {
indx <- grepl(model_name, names(sign_models))
sum(sign_models[indx]) > 0
} else {
NA
}
})
model_fits_out <- lapply(seq_along(model_fits), function (i) {
c(model_fits[[i]], significant = sign_models_fits[i])
})
attributes(model_fits_out) <- attributes(model_fits)
return (model_fits_out)
}
#' @title assessDesign
#'
#' @description This function performs simulation based trial design evaluations for a set of specified dose-response models
#'
#' @param n_patients Vector specifying the planned number of patients per dose group. A minimum of 2 patients are required in each group.
#' @param mods An object of class "Mods" as specified in the DoseFinding package.
#' @param prior_list A prior_list object specifying the utilized prior for the different dose groups
#' @param sd A positive value, specification of assumed sd
#' @param n_sim Number of simulations to be performed
#' @param alpha_crit_val (Un-adjusted) Critical value to be used for the MCP testing step. Passed to the getCritProb() function for the calculation of adjusted critical values (on the probability scale). Default is 0.05.
#' @param modeling Boolean variable defining whether the Mod part of Bayesian MCP-Mod will be performed in the assessment. More heavy on resources. Default FALSE.
#' @param simple Boolean variable defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.
#' @param avg_fit Boolean variable, defining whether an average fit (based on generalized AIC weights) should be performed in addition to the individual models. Default TRUE.
#' @param reestimate Boolean variable defining whether critical value should be calculated with re-estimated contrasts (see getCritProb function for more details). Default FALSE
#' @param contr An object of class 'optContr' as created by the getContr() function. Allows specification of a fixed contrasts matrix. Default NULL
#' @param dr_means A vector, allows specification of individual (not model based) assumed effects per dose group. Default NULL
#' @param delta A numeric value for the threshold Delta for the MED assessment. If NULL, no MED assessment is performed. Default NULL.
#' @param evidence_level A numeric value between 0 and 1 for the evidence level gamma for the MED assessment. Only required for Bayesian MED assessment, see ?getMED for details. Default NULL.
#' @param med_selection A string, either "avgFit" or "bestFit", for the method of MED selection. Default "avgFit".
#'
#' @return Returns success probabilities for the different assumed dose-response shapes, attributes also includes information around average success rate (across all assumed models) and prior Effective sample size
#'
#' @examples
#' mods <- DoseFinding::Mods(linear = NULL,
#' emax = c(0.5, 1.2),
#' exponential = 2,
#' doses = c(0, 0.5, 2,4, 8),
#' maxEff = 6)
#' sd <- 12
#' prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2),
#' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,
#' DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
#' DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2))
#' n_patients <- c(40, 60, 60, 60, 60)
#'
#' success_probabilities <- assessDesign(
#' n_patients = n_patients,
#' mods = mods,
#' prior_list = prior_list,
#' sd = sd,
#' n_sim = 1e2) # speed up example run time
#'
#' success_probabilities
#'
#' if (interactive()) { # takes typically > 5 seconds
#'
#' # with MED estimation without bootstrapping
#' # see ?getMED for details
#'
#' success_probabilities <- assessDesign(
#' n_patients = n_patients,
#' mods = mods,
#' prior_list = prior_list,
#' sd = sd,
#' modeling = TRUE,
#' n_sim = 10, # speed up example run time
#' delta = 7)
#'
#' success_probabilities
#'
#' # with MED estimation with bootstrapping
#'
#' success_probabilities <- assessDesign(
#' n_patients = n_patients,
#' mods = mods,
#' prior_list = prior_list,
#' sd = sd,
#' modeling = TRUE,
#' n_sim = 10, # speed up example run time
#' delta = 7,
#' evidence_level = 0.8)
#'
#' success_probabilities
#'
#' }
#'
#' @export
assessDesign <- function (
n_patients,
mods,
prior_list,
sd,
n_sim = 1e3,
alpha_crit_val = 0.05,
modeling = FALSE,
simple = TRUE,
avg_fit = TRUE,
reestimate = FALSE,
contr = NULL,
dr_means = NULL,
delta = NULL,
evidence_level = NULL,
med_selection = c("avgFit", "bestFit")
) {
checkmate::assert_vector(n_patients, len = length(attr(mods, "doses")), any.missing = FALSE)
checkmate::assert_double(n_patients, lower = 2, upper = Inf)
checkmate::assert_class(mods, classes = "Mods")
checkmate::assert_list(prior_list, names = "named", len = length(attr(mods, "doses")), any.missing = FALSE)
# sensitive to how DoseFinding labels their attributes for "Mods" class
checkmate::assert_integerish(n_sim, lower = 1, upper = Inf)
checkmate::assert_double(alpha_crit_val, lower = 0, upper = 1)
checkmate::assert_logical(modeling)
# TODO: check that prior_list has 'sd_tot' attribute, and that it's numeric # this is not applicable at the moment
modeling <- ifelse(!is.null(delta), TRUE, modeling)
dose_levels <- attr(mods, "doses")
data <- simulateData(
n_patients = n_patients,
dose_levels = dose_levels,
sd = sd,
mods = mods,
n_sim = n_sim,
dr_means = dr_means)
model_names <- colnames(data)[-c(1:3)]
crit_prob_adj <- getCritProb(
mods = mods,
dose_levels = dose_levels,
dose_weights = n_patients,
alpha_crit_val = alpha_crit_val)
if (!reestimate & is.null(contr)) {
contr <- getContr(
mods = mods,
dose_levels = dose_levels,
dose_weights = n_patients,
prior_list = prior_list)
}
eval_design <- lapply(model_names, function (model_name) {
posterior_list <- getPosterior(
data = getModelData(data, model_name),
prior_list = prior_list)
if (reestimate & is.null(contr)) {
post_sds <- sapply(posterior_list, function (post) summary(post)[, 2])
contr <- apply(post_sds, 2, function (post_sd) getContr(
mods = mods,
dose_levels = dose_levels,
cov_posterior = diag(post_sd^2)))
}
if (modeling) {
b_mcp_mod <- performBayesianMCPMod(
posterior_list = posterior_list,
contr = contr,
crit_prob_adj = crit_prob_adj,
simple = simple,
avg_fit = avg_fit,
delta = delta,
evidence_level = evidence_level,
med_selection = med_selection,
n_samples = n_sim)
} else {
b_mcp <- performBayesianMCP(
posterior_list = posterior_list,
contr = contr,
crit_prob_adj = crit_prob_adj)
}
})
names(eval_design) <- model_names
avg_success_rate <- mean(sapply(eval_design, function (x) {
ifelse(identical(modeling, FALSE),
attr(x, "successRate"),
attr(x$BayesianMCP, "successRate"))
}))
attr(eval_design, "avgSuccessRate") <- avg_success_rate
if (modeling & !is.null(delta)) {
avg_med_id_rate <- mean(sapply(eval_design, function (x) {
mean(attr(x, "MED")[, "med_reached"])
}))
attr(eval_design, "avgMEDIdentificationRate") <- avg_med_id_rate
}
attr(eval_design, "placEff") <- ifelse(test = is.null(dr_means),
yes = attr(mods, "placEff"),
no = dr_means[1])
attr(eval_design, "maxEff") <- ifelse(test = is.null(dr_means),
yes = attr(mods, "maxEff"),
no = diff(range(dr_means)))
attr(eval_design, "sampleSize") <- n_patients
attr(eval_design, "priorESS") <- round(getESS(prior_list), 1)
return (eval_design)
}
BayesMCPi <- function (
posterior_i,
contr,
crit_prob_adj
) {
post_comps_i <- attr(posterior_i, "posteriorInfo")
## Pinheiro, José, et al.
## "Model‐based dose finding under model uncertainty using general parametric models."
## Statistics in medicine 33.10 (2014): 1646-1661. doi:10.1002/sim.6052
## Section 2.2
mu_mat <- do.call(cbind, post_comps_i$means) # [N candidate models, N posterior components]
S_list <- post_comps_i$covMats # list[[N posterior components]] of cov matrices
# with [N candidate models, N candidate models]
# Optimal contrast matrix
contr_mat <- contr$contMat # [N dose levels, N candidate models]
## Contrast estimates for each component
contr_est_mu <- t(contr_mat) %*% mu_mat # [N candidate models, N posterior components]
contr_est_cov <-
lapply(S_list, function(S) # list[[N posterior components]] of cov matrices
t(contr_mat) %*% S %*% contr_mat) # with [N candidate models, N candidate models]
## Test statistic for each component # [N candidate models, N posterior components]
z_m <- contr_est_mu / do.call(cbind, lapply(contr_est_cov, diag))^0.5
## Fleischer, Frank, et al.
## "Bayesian MCPMod."
## Pharmaceutical Statistics 21.3 (2022): 654-670. doi:10.1002/pst.2193
## Section 2.3
## Posterior probabilities # [N candidate models]
post_probs <- (stats::pnorm(z_m) %*% unlist(post_comps_i$weights))[, 1]
res <- c(sign = ifelse(max(post_probs) > crit_prob_adj, 1, 0),
crit_prob_adj = crit_prob_adj,
max_post_prob = max(post_probs),
post_probs = post_probs)
return (res)
# getPostProb <- function (
#
# contr_j, # j: dose level
# post_combs_i # i: simulation outcome
#
# ) {
#
# ## Fleischer, Frank, et al. "Bayesian MCPMod."
# ## Pharmaceutical Statistics 21.3 (2022): 654-670. doi:10.1002/pst.2193
# ## Sections 2.2 and 2.3
#
# ## Test statistic = sum over all components of
# ## posterior weight * normal probability distribution of
# ## critical values for doses * estimated mean / sqrt(product of critical values for doses)
#
# ## Calculation for each component of the posterior
# contr_theta <- apply(post_combs_i$means, 1, `%*%`, contr_j)
# contr_var <- apply(post_combs_i$vars, 1, `%*%`, contr_j^2)
# contr_weights <- post_combs_i$weights
#
# ## P(c_m * theta > 0 | Y = y) for a shape m (and dose j)
# post_probs <- sum(contr_weights * stats::pnorm(contr_theta / sqrt(contr_var)))
#
# return (post_probs)
#
# }
#
# post_combs_i <- getPostCompsI(posterior_i)
# post_probs <- apply(contr$contMat, 2, getPostProb, post_combs_i)
#
# res <- c(sign = ifelse(max(post_probs) > crit_prob_adj, 1, 0),
# crit_prob_adj = crit_prob_adj,
# max_post_prob = max(post_probs),
# post_probs = post_probs)
#
# return (res)
}
#' @title getContr
#'
#' @description This function calculates contrast vectors that are optimal for detecting certain alternatives via applying the function optContr() of the DoseFinding package.
#' Hereby, 4 different options can be distinguished that are automatically executed based on the input that is provided
#' 1) Bayesian approach: If dose_weights and a prior_list are provided an optimized contrasts for the posterior sample size is calculated.
#' In detail, in a first step the dose_weights (typically the number of patients per dose group) and the prior information is combined by calculating for
#' each dose group a posterior effective sample. Based on this posterior effective sample sizes the allocation ratio is derived, which allows for a calculation on
#' pseudo-optimal contrasts via regular MCPMod are calculated from the
#' regular MCPMod for these specific weights
#' 2) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the
#' regular MCPMod for these specific weights
#' 3) Bayesian approach + re-estimation: If only a cov_posterior (i.e. variability of the posterior distribution) is provided, pseudo-optimal contrasts based on these posterior weights will be calculated
#' 4) Frequentist approach+re-estimation: If only a cov_new_trial (i.e. the estimated variability of a new trial) is provided, optimal contrast vectors are calculated from the
#' regular MCPMod for this specific covariance matrix.
#'
#' @param mods An object of class 'Mods' as created by the function 'DoseFinding::Mods()'
#' @param dose_levels Vector containing the different dosage levels.
#' @param dose_weights Vector specifying weights for the different doses. Please note that in case this information is provided together with a prior (i.e. Option 1) is planned these two inputs should be provided on the same scale (e.g. patient numbers). Default NULL
#' @param prior_list A list of objects of class 'normMix' as created with 'RBesT::mixnorm()'. Only required as input for Option 1. Default NULL
#' @param cov_posterior A covariance matrix with information about the variability of the posterior distribution, only required for Option 3. Default NULL
#' @param cov_new_trial A covariance matrix with information about the observed variability, only required for Option 4. Default NULL
#'
#' @examples
#' dose_levels <- c(0, 0.5, 2, 4, 8)
#' mods <- DoseFinding::Mods(
#' linear = NULL,
#' emax = c(0.5, 1.2),
#' exponential = 2,
#' doses = dose_levels,
#' maxEff = 6)
#' cov_posterior <- diag(c(2.8, 3, 2.5, 3.5, 4)^2)
#'
#' contr_mat <- getContr(
#' mods = mods,
#' dose_levels = dose_levels,
#' cov_posterior = cov_posterior)
#'
#' @return An object of class 'optContr' as provided by the function 'DoseFinding::optContr()'.
#'
#' @export
getContr <- function (
mods,
dose_levels,
dose_weights = NULL,
prior_list = NULL,
cov_posterior = NULL,
cov_new_trial = NULL
) {
checkmate::assert_class(mods, "Mods")
checkmate::assert_numeric(dose_levels, lower = 0, any.missing = FALSE, len = length(attr(mods, "doses")))
if (!is.null(dose_weights)) {
checkmate::assert_numeric(dose_weights, any.missing = FALSE, len = length(attr(mods, "doses")))
}
if (!is.null(prior_list)) {
checkmate::assert_list(prior_list, names = "named", len = length(attr(mods, "doses")), any.missing = FALSE)
}
# Determine scenario
if (!is.null(cov_new_trial) &&
is.null(dose_weights) && is.null(prior_list) && is.null(cov_posterior)) {
w <- NULL
S <- cov_new_trial
} else if (!is.null(dose_weights) &&
is.null(cov_new_trial) && is.null(prior_list) && is.null(cov_posterior)) {
w <- dose_weights
S <- NULL
} else if (!is.null(cov_posterior) &&
is.null(cov_new_trial) && is.null(prior_list) && is.null(dose_weights)) {
w <- NULL
S <- cov_posterior
} else if (!is.null(dose_weights) && !is.null(prior_list) &&
is.null(cov_new_trial) && is.null(cov_posterior)) {
w <- dose_weights +
suppressMessages(round(unlist(lapply(prior_list, RBesT::ess))))
S <- NULL
} else {
stop("Invalid combination of inputs. Valid combinations are:
1. cov_new_trial only (frequentist with re-estimation)
2. dose_weights only (frequentist without re-estimation)
3. cov_posterior only (Bayesian with re-estimation)
4. dose_weights + prior_list (Bayesian without re-estimation)")
}
# Compute contrasts
if (is.null(w)) {
contr <- DoseFinding::optContr(
models = mods,
doses = dose_levels,
S = S
)
} else {
contr <- DoseFinding::optContr(
models = mods,
doses = dose_levels,
w = w
)
}
return(contr)
}
#' @title getCritProb
#'
#' @description This function calculates multiplicity adjusted critical values. The critical values are calculated in such a way that
#' when using non-informative priors the actual error level for falsely declaring a significant trial in the Bayesian MCPMod is controlled (by the specified alpha level).
#' Hereby optimal contrasts of the frequentist MCPMod are applied and two options can be distinguished
#' 1) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the
#' regular MCPMod for these specific weights and the corresponding critical value for this set of contrasts is calculated via the critVal() function of the DoseFinding package.
#' 2) Frequentist approach + re-estimation: If only a cov_new_trial (i.e. the covariance matrix of a new trial) is provided, optimal contrast vectors are calculated from the
#' regular MCPMod for this specific matrix. Here as well the critical value for this set of contrasts is calculated via the critVal() function of the DoseFinding package.
#'
#' @param mods An object of class "Mods" as specified in the DoseFinding package.
#' @param dose_levels Vector containing the different dosage levels.
#' @param dose_weights Vector specifying weights for the different doses, only required for Option i). Default NULL
#' @param cov_new_trial A covariance matrix, only required for Option ii). Default NULL
#' @param alpha_crit_val Significance level. Default set to 0.025.
#'
#' @examples
#' mods <- DoseFinding::Mods(linear = NULL,
#' emax = c(0.5, 1.2),
#' exponential = 2,
#' doses = c(0, 0.5, 2,4, 8))
#' dose_levels <- c(0, 0.5, 2, 4, 8)
#' critVal <- getCritProb(
#' mods = mods,
#' dose_weights = c(50,50,50,50,50), #reflecting the planned sample size
#' dose_levels = dose_levels,
#' alpha_crit_val = 0.05)
#' @return Multiplicity adjusted critical value on the probability scale.
#'
#' @export
getCritProb <- function (
mods,
dose_levels,
dose_weights = NULL,
cov_new_trial = NULL,
alpha_crit_val = 0.025
) {
checkmate::assert_class(mods, classes = "Mods")
checkmate::assert_double(dose_levels, lower = 0, any.missing = FALSE)
checkmate::assert_double(dose_weights, any.missing = FALSE, len = length(dose_levels), null.ok = TRUE)
checkmate::assert_double(alpha_crit_val, lower = 0, upper = 1)
# Get contrast using updated getContr
contr <- getContr(
mods = mods,
dose_levels = dose_levels,
dose_weights = dose_weights,
cov_new_trial = cov_new_trial
)
# Calculate critical probability
crit_prob <- stats::pnorm(DoseFinding::critVal(
corMat = contr$corMat,
alpha = alpha_crit_val,
df = 0,
alternative = "one.sided"
))
return(crit_prob)
}
getModelSuccesses <- function (b_mcp) {
stopifnot(inherits(b_mcp, "BayesianMCP"))
model_indices <- grepl("post_probs.", colnames(b_mcp))
model_names <- colnames(b_mcp)[model_indices] |>
sub(pattern = "post_probs.", replacement = "", x = _)
model_successes <- colMeans(b_mcp[, model_indices] > b_mcp[, "crit_prob_adj"])
names(model_successes) <- model_names
return (model_successes)
}
#' @title performBayesianMCPMod
#'
#' @description Performs Bayesian MCP Test step and modeling in a combined fashion. See performBayesianMCP() function for MCP Test step and getModelFits() for the modeling step
#'
#' @param posterior_list An object of class 'postList' or a list of 'postList' objects as created by getPosterior() containing information about the (mixture) posterior distribution per dose group
#' @param contr An object of class 'optContr' as created by the getContr() function. It contains the contrast matrix to be used for the testing step.
#' @param crit_prob_adj A getCritProb object, specifying the critical value to be used for the testing (on the probability scale).
#' @param simple Boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits() function. Default FALSE.
#' @param avg_fit Boolean variable, defining whether an average fit (based on generalized AIC weights) should be performed in addition to the individual models. Default TRUE.
#' @param delta A numeric value for the threshold Delta for the MED assessment. If NULL, no MED assessment is performed. Default NULL.
#' @param evidence_level A numeric value between 0 and 1 for the evidence level gamma for the MED assessment. Only required for Bayesian MED assessment, see ?getMED for details. Default NULL.
#' @param med_selection A string, either "avgFit" or "bestFit" based on the lowest gAIC, for the method of MED selection. Default "avgFit".
#' @param n_samples A numerical for the number of bootstrapped samples in case the Bayesian MED assessment is performed. Default 1e3.
#' @examples
#' mods <- DoseFinding::Mods(linear = NULL,
#' emax = c(0.5, 1.2),
#' exponential = 2,
#' doses = c(0, 0.5, 2,4, 8))
#' dose_levels <- c(0, 0.5, 2, 4, 8)
#' sd_posterior <- c(2.8, 3, 2.5, 3.5, 4)
#' contr_mat <- getContr(
#' mods = mods,
#' dose_levels = dose_levels,
#' cov_posterior = diag(sd_posterior)^2)
#' critVal <- getCritProb(
#' mods = mods,
#' dose_weights = c(50, 50, 50, 50, 50), #reflecting the planned sample size
#' dose_levels = dose_levels,
#' alpha_crit_val = 0.6) # unreasonable alpha for this example, rather choose 0.05
#' prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2),
#' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,
#' DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
#' DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2))
#' mu <- c(0, 1, 1.5, 2, 2.5)
#' S_hat <- diag(c(5, 4, 6, 7, 8)^2)
#' posterior_list <- getPosterior(
#' prior_list = prior_list,
#' mu_hat = mu,
#' S_hat = S_hat,
#' calc_ess = TRUE)
#'
#' performBayesianMCPMod(posterior_list = posterior_list,
#' contr = contr_mat,
#' crit_prob_adj = critVal,
#' simple = FALSE,
#' delta = 1)
#'
#' @return Bayesian MCP test result as well as modeling result.
#'
#' @export
performBayesianMCPMod <- function (
posterior_list,
contr,
crit_prob_adj,
simple = FALSE,
avg_fit = TRUE,
delta = NULL,
evidence_level = NULL,
med_selection = c("avgFit", "bestFit"),
n_samples = 1e3
) {
if (inherits(posterior_list, "postList")) {
posterior_list <- list(posterior_list)
}
checkmate::assert_list(posterior_list, types = "postList")
checkmate::assert_class(contr, "optContr")
checkmate::assert_class(crit_prob_adj, "numeric")
checkmate::assert_logical(simple)
checkmate::assert_logical(avg_fit)
checkmate::assert_double(delta, null.ok = TRUE)
checkmate::assert_double(evidence_level, lower = 0, upper = 1, null.ok = TRUE)
checkmate::assert_character(med_selection)
checkmate::assert_double(n_samples, lower = 1)
med_selection <- match.arg(med_selection, choices = c("avgFit", "bestFit"))
get_med <- ifelse(is.null(delta), FALSE, TRUE)
if (!is.null(evidence_level)) {
stopifnot("delta must not be NULL if evidence_level is not NULL" = !is.null(delta))
}
if (inherits(contr, "optContr")) {
model_names <- colnames(contr$contMat)
dose_levels <- as.numeric(rownames(contr$contMat))
} else if (length(contr) == length(posterior_list)) {
model_names <- colnames(contr[[1]]$contMat)
dose_levels <- as.numeric(rownames(contr[[1]]$contMat))
} else {
stop ("Argument 'contr' must be of type 'optContr'")
}
b_mcp <- performBayesianMCP(
posterior_list = posterior_list,
contr = contr,
crit_prob_adj = crit_prob_adj)
b_mod <- performBayesianMod(
b_mcp = b_mcp,
posterior_list = posterior_list,
model_names = model_names,
dose_levels = dose_levels,
simple = simple,
avg_fit = avg_fit)
b_mcp_mod <- list(BayesianMCP = b_mcp, Mod = b_mod)
if (get_med) {
med_info <- t(sapply(b_mcp_mod$Mod, function (model_fits) {
if (!is.null(model_fits)) {
if (!is.null(evidence_level)) {
bs_quantiles <- getBootstrapQuantiles(
model_fits = model_fits,
quantiles = 1 - evidence_level,
n_samples = n_samples)
med_info <- getMED(
delta = delta,
evidence_level = evidence_level,
bs_quantiles = bs_quantiles)
} else {
med_info <- getMED(
delta = delta,
model_fits = model_fits)
}
if (med_selection == "avgFit") {
return (med_info[, "avgFit"])
} else if (med_selection == "bestFit") {
model_min <- which.min(sapply(model_fits, function (x) x$gAIC))
return (med_info[, model_min])
}
} else {
return (c(med_reached = 0, med = NA))
}
}))
attr(b_mcp_mod, "MED") <- med_info
attr(b_mcp_mod, "MEDSelection") <- med_selection
}
class(b_mcp_mod) <- "BayesianMCPMod"
return (b_mcp_mod)
}
#' @title performBayesianMCP
#'
#' @description Performs Bayesian MCP Test step, as described in Fleischer et al. (2022).
#' Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate the posterior probability is calculated that the contrast is bigger than 0 (based on the posterior distribution of the dose groups).
#' In order to obtain significant test decision we consider the maximum of the posterior probabilities across the different models. This maximum is compared with a (multiplicity adjusted) critical value (on the probability scale).
#' @references Fleischer F, Bossert S, Deng Q, Loley C, Gierse J. 2022. Bayesian MCPMod. Pharmaceutical Statistics. 21(3): 654-670. doi:10.1002/pst.2193
#' @param posterior_list An object derived with getPosterior with information about the (mixture) posterior distribution per dose group
#' @param contr An object of class 'optContr' as created by the getContr() function. It contains the contrast matrix to be used for the testing step.
#' @param crit_prob_adj A getCritProb object, specifying the critical value to be used for the testing (on the probability scale)
#'
#' @examples
#' mods <- DoseFinding::Mods(linear = NULL,
#' emax = c(0.5, 1.2),
#' exponential = 2,
#' doses = c(0, 0.5, 2,4, 8))
#' dose_levels <- c(0, 0.5, 2, 4, 8)
#' sd_posterior <- c(2.8,3,2.5,3.5,4)
#' contr_mat <- getContr(
#' mods = mods,
#' dose_levels = dose_levels,
#' cov_posterior = diag(sd_posterior)^2)
#' critVal <- getCritProb(
#' mods = mods,
#' dose_weights = c(50, 50, 50, 50, 50), #reflecting the planned sample size
#' dose_levels = dose_levels,
#' alpha_crit_val = 0.05)
#' prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2),
#' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,
#' DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
#' DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2))
#' mu <- c(0, 1, 1.5, 2, 2.5)
#' S_hat <- diag(c(5, 4, 6, 7, 8)^2)
#' posterior_list <- getPosterior(
#' prior_list = prior_list,
#' mu_hat = mu,
#' S_hat = S_hat,
#' calc_ess = TRUE)
#'
#' performBayesianMCP(posterior_list = posterior_list,
#' contr = contr_mat,
#' crit_prob_adj = critVal)
#'
#' @return Bayesian MCP test result, with information about p-values for the individual dose-response shapes and overall significance
#'
#' @export
performBayesianMCP <- function(
posterior_list,
contr,
crit_prob_adj
) {
if (inherits(posterior_list, "postList")) {
posterior_list <- list(posterior_list)
}
checkmate::assert_list(posterior_list, types = "postList")
checkmate::assert(
checkmate::check_class(contr, "optContr"),
checkmate::check_list(contr, types = "optContr"),
combine = "or")
checkmate::assert_class(crit_prob_adj, "numeric")
checkmate::assert_numeric(crit_prob_adj, lower = 0, upper = Inf)
## TODO: What case is covered by the IF clause?
## Optionally providing custom contrast?
## This was written before the checkmate assertions were included
## -> if contrasts are re-evaluated, contr becomes a list of contrasts
if (inherits(contr, "optContr")) {
b_mcp <- t(sapply(posterior_list, BayesMCPi, contr, crit_prob_adj))
} else {
b_mcp <- t(mapply(BayesMCPi, posterior_list, contr, crit_prob_adj))
}
class(b_mcp) <- "BayesianMCP"
attr(b_mcp, "critProbAdj") <- crit_prob_adj
attr(b_mcp, "successRate") <- mean(b_mcp[, 1])
attr(b_mcp, "essAvg") <- ifelse(
test = is.na(attr(posterior_list[[1]], "ess")),
yes = numeric(0),
no = rowMeans(sapply(posterior_list,
function (posteriors) attr(posteriors, "ess"))))
return (b_mcp)
}
performBayesianMod <- function (
b_mcp,
posterior_list,
model_names,
dose_levels,
simple,
avg_fit
) {
## Make parallel processing optional
if (requireNamespace("future.apply", quietly = TRUE)) {
optPar_lapply <- future.apply::future_lapply
} else {
optPar_lapply <- lapply
}
fits_list <- optPar_lapply(seq_along(posterior_list), function (i) {
if (b_mcp[i, 1]) {
model_fits <- getModelFits(
models = model_names,
dose_levels = dose_levels,
posterior = posterior_list[[i]],
avg_fit = avg_fit,
simple = simple)
sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "critProbAdj")
model_fits <- addSignificance(model_fits, sign_models)
} else {
NULL
}
})
}
Try the BayesianMCPMod package in your browser
Any scripts or data that you put into this service are public.
BayesianMCPMod documentation built on Aug. 29, 2025, 5:13 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions performBayesianMCP
Documented in performBayesianMCP
addSignificance <- function (
model_fits,
sign_models
) {
sign_models_fits <- sapply(names(model_fits), function (model_name) {
if (model_name != "avgFit") {
indx <- grepl(model_name, names(sign_models))
sum(sign_models[indx]) > 0
} else {
NA
}
})
model_fits_out <- lapply(seq_along(model_fits), function (i) {
c(model_fits[[i]], significant = sign_models_fits[i])
})
attributes(model_fits_out) <- attributes(model_fits)
return (model_fits_out)
}
#' @title assessDesign
#'
#' @description This function performs simulation based trial design evaluations for a set of specified dose-response models
#'
#' @param n_patients Vector specifying the planned number of patients per dose group. A minimum of 2 patients are required in each group.
#' @param mods An object of class "Mods" as specified in the DoseFinding package.
#' @param prior_list A prior_list object specifying the utilized prior for the different dose groups
#' @param sd A positive value, specification of assumed sd
#' @param n_sim Number of simulations to be performed
#' @param alpha_crit_val (Un-adjusted) Critical value to be used for the MCP testing step. Passed to the getCritProb() function for the calculation of adjusted critical values (on the probability scale). Default is 0.05.
#' @param modeling Boolean variable defining whether the Mod part of Bayesian MCP-Mod will be performed in the assessment. More heavy on resources. Default FALSE.
#' @param simple Boolean variable defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.
#' @param avg_fit Boolean variable, defining whether an average fit (based on generalized AIC weights) should be performed in addition to the individual models. Default TRUE.
#' @param reestimate Boolean variable defining whether critical value should be calculated with re-estimated contrasts (see getCritProb function for more details). Default FALSE
#' @param contr An object of class 'optContr' as created by the getContr() function. Allows specification of a fixed contrasts matrix. Default NULL
#' @param dr_means A vector, allows specification of individual (not model based) assumed effects per dose group. Default NULL
#' @param delta A numeric value for the threshold Delta for the MED assessment. If NULL, no MED assessment is performed. Default NULL.
#' @param evidence_level A numeric value between 0 and 1 for the evidence level gamma for the MED assessment. Only required for Bayesian MED assessment, see ?getMED for details. Default NULL.
#' @param med_selection A string, either "avgFit" or "bestFit", for the method of MED selection. Default "avgFit".
#'
#' @return Returns success probabilities for the different assumed dose-response shapes, attributes also includes information around average success rate (across all assumed models) and prior Effective sample size
#'
#' @examples
#' mods <- DoseFinding::Mods(linear = NULL,
#' emax = c(0.5, 1.2),
#' exponential = 2,
#' doses = c(0, 0.5, 2,4, 8),
#' maxEff = 6)
#' sd <- 12
#' prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2),
#' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,
#' DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
#' DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2))
#' n_patients <- c(40, 60, 60, 60, 60)
#'
#' success_probabilities <- assessDesign(
#' n_patients = n_patients,
#' mods = mods,
#' prior_list = prior_list,
#' sd = sd,
#' n_sim = 1e2) # speed up example run time
#'
#' success_probabilities
#'
#' if (interactive()) { # takes typically > 5 seconds
#'
#' # with MED estimation without bootstrapping
#' # see ?getMED for details
#'
#' success_probabilities <- assessDesign(
#' n_patients = n_patients,
#' mods = mods,
#' prior_list = prior_list,
#' sd = sd,
#' modeling = TRUE,
#' n_sim = 10, # speed up example run time
#' delta = 7)
#'
#' success_probabilities
#'
#' # with MED estimation with bootstrapping
#'
#' success_probabilities <- assessDesign(
#' n_patients = n_patients,
#' mods = mods,
#' prior_list = prior_list,
#' sd = sd,
#' modeling = TRUE,
#' n_sim = 10, # speed up example run time
#' delta = 7,
#' evidence_level = 0.8)
#'
#' success_probabilities
#'
#' }
#'
#' @export
assessDesign <- function (
n_patients,
mods,
prior_list,
sd,
n_sim = 1e3,
alpha_crit_val = 0.05,
modeling = FALSE,
simple = TRUE,
avg_fit = TRUE,
reestimate = FALSE,
contr = NULL,
dr_means = NULL,
delta = NULL,
evidence_level = NULL,
med_selection = c("avgFit", "bestFit")
) {
checkmate::assert_vector(n_patients, len = length(attr(mods, "doses")), any.missing = FALSE)
checkmate::assert_double(n_patients, lower = 2, upper = Inf)
checkmate::assert_class(mods, classes = "Mods")
checkmate::assert_list(prior_list, names = "named", len = length(attr(mods, "doses")), any.missing = FALSE)
# sensitive to how DoseFinding labels their attributes for "Mods" class
checkmate::assert_integerish(n_sim, lower = 1, upper = Inf)
checkmate::assert_double(alpha_crit_val, lower = 0, upper = 1)
checkmate::assert_logical(modeling)
# TODO: check that prior_list has 'sd_tot' attribute, and that it's numeric # this is not applicable at the moment
modeling <- ifelse(!is.null(delta), TRUE, modeling)
dose_levels <- attr(mods, "doses")
data <- simulateData(
n_patients = n_patients,
dose_levels = dose_levels,
sd = sd,
mods = mods,
n_sim = n_sim,
dr_means = dr_means)
model_names <- colnames(data)[-c(1:3)]
crit_prob_adj <- getCritProb(
mods = mods,
dose_levels = dose_levels,
dose_weights = n_patients,
alpha_crit_val = alpha_crit_val)
if (!reestimate & is.null(contr)) {
contr <- getContr(
mods = mods,
dose_levels = dose_levels,
dose_weights = n_patients,
prior_list = prior_list)
}
eval_design <- lapply(model_names, function (model_name) {
posterior_list <- getPosterior(
data = getModelData(data, model_name),
prior_list = prior_list)
if (reestimate & is.null(contr)) {
post_sds <- sapply(posterior_list, function (post) summary(post)[, 2])
contr <- apply(post_sds, 2, function (post_sd) getContr(
mods = mods,
dose_levels = dose_levels,
cov_posterior = diag(post_sd^2)))
}
if (modeling) {
b_mcp_mod <- performBayesianMCPMod(
posterior_list = posterior_list,
contr = contr,
crit_prob_adj = crit_prob_adj,
simple = simple,
avg_fit = avg_fit,
delta = delta,
evidence_level = evidence_level,
med_selection = med_selection,
n_samples = n_sim)
} else {
b_mcp <- performBayesianMCP(
posterior_list = posterior_list,
contr = contr,
crit_prob_adj = crit_prob_adj)
}
})
names(eval_design) <- model_names
avg_success_rate <- mean(sapply(eval_design, function (x) {
ifelse(identical(modeling, FALSE),
attr(x, "successRate"),
attr(x$BayesianMCP, "successRate"))
}))
attr(eval_design, "avgSuccessRate") <- avg_success_rate
if (modeling & !is.null(delta)) {
avg_med_id_rate <- mean(sapply(eval_design, function (x) {
mean(attr(x, "MED")[, "med_reached"])
}))
attr(eval_design, "avgMEDIdentificationRate") <- avg_med_id_rate
}
attr(eval_design, "placEff") <- ifelse(test = is.null(dr_means),
yes = attr(mods, "placEff"),
no = dr_means[1])
attr(eval_design, "maxEff") <- ifelse(test = is.null(dr_means),
yes = attr(mods, "maxEff"),
no = diff(range(dr_means)))
attr(eval_design, "sampleSize") <- n_patients
attr(eval_design, "priorESS") <- round(getESS(prior_list), 1)
return (eval_design)
}
BayesMCPi <- function (
posterior_i,
contr,
crit_prob_adj
) {
post_comps_i <- attr(posterior_i, "posteriorInfo")
## Pinheiro, José, et al.
## "Model‐based dose finding under model uncertainty using general parametric models."
## Statistics in medicine 33.10 (2014): 1646-1661. doi:10.1002/sim.6052
## Section 2.2
mu_mat <- do.call(cbind, post_comps_i$means) # [N candidate models, N posterior components]
S_list <- post_comps_i$covMats # list[[N posterior components]] of cov matrices
# with [N candidate models, N candidate models]
# Optimal contrast matrix
contr_mat <- contr$contMat # [N dose levels, N candidate models]
## Contrast estimates for each component
contr_est_mu <- t(contr_mat) %*% mu_mat # [N candidate models, N posterior components]
contr_est_cov <-
lapply(S_list, function(S) # list[[N posterior components]] of cov matrices
t(contr_mat) %*% S %*% contr_mat) # with [N candidate models, N candidate models]
## Test statistic for each component # [N candidate models, N posterior components]
z_m <- contr_est_mu / do.call(cbind, lapply(contr_est_cov, diag))^0.5
## Fleischer, Frank, et al.
## "Bayesian MCPMod."
## Pharmaceutical Statistics 21.3 (2022): 654-670. doi:10.1002/pst.2193
## Section 2.3
## Posterior probabilities # [N candidate models]
post_probs <- (stats::pnorm(z_m) %*% unlist(post_comps_i$weights))[, 1]
res <- c(sign = ifelse(max(post_probs) > crit_prob_adj, 1, 0),
crit_prob_adj = crit_prob_adj,
max_post_prob = max(post_probs),
post_probs = post_probs)
return (res)
# getPostProb <- function (
#
# contr_j, # j: dose level
# post_combs_i # i: simulation outcome
#
# ) {
#
# ## Fleischer, Frank, et al. "Bayesian MCPMod."
# ## Pharmaceutical Statistics 21.3 (2022): 654-670. doi:10.1002/pst.2193
# ## Sections 2.2 and 2.3
#
# ## Test statistic = sum over all components of
# ## posterior weight * normal probability distribution of
# ## critical values for doses * estimated mean / sqrt(product of critical values for doses)
#
# ## Calculation for each component of the posterior
# contr_theta <- apply(post_combs_i$means, 1, `%*%`, contr_j)
# contr_var <- apply(post_combs_i$vars, 1, `%*%`, contr_j^2)
# contr_weights <- post_combs_i$weights
#
# ## P(c_m * theta > 0 | Y = y) for a shape m (and dose j)
# post_probs <- sum(contr_weights * stats::pnorm(contr_theta / sqrt(contr_var)))
#
# return (post_probs)
#
# }
#
# post_combs_i <- getPostCompsI(posterior_i)
# post_probs <- apply(contr$contMat, 2, getPostProb, post_combs_i)
#
# res <- c(sign = ifelse(max(post_probs) > crit_prob_adj, 1, 0),
# crit_prob_adj = crit_prob_adj,
# max_post_prob = max(post_probs),
# post_probs = post_probs)
#
# return (res)
}
#' @title getContr
#'
#' @description This function calculates contrast vectors that are optimal for detecting certain alternatives via applying the function optContr() of the DoseFinding package.
#' Hereby, 4 different options can be distinguished that are automatically executed based on the input that is provided
#' 1) Bayesian approach: If dose_weights and a prior_list are provided an optimized contrasts for the posterior sample size is calculated.
#' In detail, in a first step the dose_weights (typically the number of patients per dose group) and the prior information is combined by calculating for
#' each dose group a posterior effective sample. Based on this posterior effective sample sizes the allocation ratio is derived, which allows for a calculation on
#' pseudo-optimal contrasts via regular MCPMod are calculated from the
#' regular MCPMod for these specific weights
#' 2) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the
#' regular MCPMod for these specific weights
#' 3) Bayesian approach + re-estimation: If only a cov_posterior (i.e. variability of the posterior distribution) is provided, pseudo-optimal contrasts based on these posterior weights will be calculated
#' 4) Frequentist approach+re-estimation: If only a cov_new_trial (i.e. the estimated variability of a new trial) is provided, optimal contrast vectors are calculated from the
#' regular MCPMod for this specific covariance matrix.
#'
#' @param mods An object of class 'Mods' as created by the function 'DoseFinding::Mods()'
#' @param dose_levels Vector containing the different dosage levels.
#' @param dose_weights Vector specifying weights for the different doses. Please note that in case this information is provided together with a prior (i.e. Option 1) is planned these two inputs should be provided on the same scale (e.g. patient numbers). Default NULL
#' @param prior_list A list of objects of class 'normMix' as created with 'RBesT::mixnorm()'. Only required as input for Option 1. Default NULL
#' @param cov_posterior A covariance matrix with information about the variability of the posterior distribution, only required for Option 3. Default NULL
#' @param cov_new_trial A covariance matrix with information about the observed variability, only required for Option 4. Default NULL
#'
#' @examples
#' dose_levels <- c(0, 0.5, 2, 4, 8)
#' mods <- DoseFinding::Mods(
#' linear = NULL,
#' emax = c(0.5, 1.2),
#' exponential = 2,
#' doses = dose_levels,
#' maxEff = 6)
#' cov_posterior <- diag(c(2.8, 3, 2.5, 3.5, 4)^2)
#'
#' contr_mat <- getContr(
#' mods = mods,
#' dose_levels = dose_levels,
#' cov_posterior = cov_posterior)
#'
#' @return An object of class 'optContr' as provided by the function 'DoseFinding::optContr()'.
#'
#' @export
getContr <- function (
mods,
dose_levels,
dose_weights = NULL,
prior_list = NULL,
cov_posterior = NULL,
cov_new_trial = NULL
) {
checkmate::assert_class(mods, "Mods")
checkmate::assert_numeric(dose_levels, lower = 0, any.missing = FALSE, len = length(attr(mods, "doses")))
if (!is.null(dose_weights)) {
checkmate::assert_numeric(dose_weights, any.missing = FALSE, len = length(attr(mods, "doses")))
}
if (!is.null(prior_list)) {
checkmate::assert_list(prior_list, names = "named", len = length(attr(mods, "doses")), any.missing = FALSE)
}
# Determine scenario
if (!is.null(cov_new_trial) &&
is.null(dose_weights) && is.null(prior_list) && is.null(cov_posterior)) {
w <- NULL
S <- cov_new_trial
} else if (!is.null(dose_weights) &&
is.null(cov_new_trial) && is.null(prior_list) && is.null(cov_posterior)) {
w <- dose_weights
S <- NULL
} else if (!is.null(cov_posterior) &&
is.null(cov_new_trial) && is.null(prior_list) && is.null(dose_weights)) {
w <- NULL
S <- cov_posterior
} else if (!is.null(dose_weights) && !is.null(prior_list) &&
is.null(cov_new_trial) && is.null(cov_posterior)) {
w <- dose_weights +
suppressMessages(round(unlist(lapply(prior_list, RBesT::ess))))
S <- NULL
} else {
stop("Invalid combination of inputs. Valid combinations are:
1. cov_new_trial only (frequentist with re-estimation)
2. dose_weights only (frequentist without re-estimation)
3. cov_posterior only (Bayesian with re-estimation)
4. dose_weights + prior_list (Bayesian without re-estimation)")
}
# Compute contrasts
if (is.null(w)) {
contr <- DoseFinding::optContr(
models = mods,
doses = dose_levels,
S = S
)
} else {
contr <- DoseFinding::optContr(
models = mods,
doses = dose_levels,
w = w
)
}
return(contr)
}
#' @title getCritProb
#'
#' @description This function calculates multiplicity adjusted critical values. The critical values are calculated in such a way that
#' when using non-informative priors the actual error level for falsely declaring a significant trial in the Bayesian MCPMod is controlled (by the specified alpha level).
#' Hereby optimal contrasts of the frequentist MCPMod are applied and two options can be distinguished
#' 1) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the
#' regular MCPMod for these specific weights and the corresponding critical value for this set of contrasts is calculated via the critVal() function of the DoseFinding package.
#' 2) Frequentist approach + re-estimation: If only a cov_new_trial (i.e. the covariance matrix of a new trial) is provided, optimal contrast vectors are calculated from the
#' regular MCPMod for this specific matrix. Here as well the critical value for this set of contrasts is calculated via the critVal() function of the DoseFinding package.
#'
#' @param mods An object of class "Mods" as specified in the DoseFinding package.
#' @param dose_levels Vector containing the different dosage levels.
#' @param dose_weights Vector specifying weights for the different doses, only required for Option i). Default NULL
#' @param cov_new_trial A covariance matrix, only required for Option ii). Default NULL
#' @param alpha_crit_val Significance level. Default set to 0.025.
#'
#' @examples
#' mods <- DoseFinding::Mods(linear = NULL,
#' emax = c(0.5, 1.2),
#' exponential = 2,
#' doses = c(0, 0.5, 2,4, 8))
#' dose_levels <- c(0, 0.5, 2, 4, 8)
#' critVal <- getCritProb(
#' mods = mods,
#' dose_weights = c(50,50,50,50,50), #reflecting the planned sample size
#' dose_levels = dose_levels,
#' alpha_crit_val = 0.05)
#' @return Multiplicity adjusted critical value on the probability scale.
#'
#' @export
getCritProb <- function (
mods,
dose_levels,
dose_weights = NULL,
cov_new_trial = NULL,
alpha_crit_val = 0.025
) {
checkmate::assert_class(mods, classes = "Mods")
checkmate::assert_double(dose_levels, lower = 0, any.missing = FALSE)
checkmate::assert_double(dose_weights, any.missing = FALSE, len = length(dose_levels), null.ok = TRUE)
checkmate::assert_double(alpha_crit_val, lower = 0, upper = 1)
# Get contrast using updated getContr
contr <- getContr(
mods = mods,
dose_levels = dose_levels,
dose_weights = dose_weights,
cov_new_trial = cov_new_trial
)
# Calculate critical probability
crit_prob <- stats::pnorm(DoseFinding::critVal(
corMat = contr$corMat,
alpha = alpha_crit_val,
df = 0,
alternative = "one.sided"
))
return(crit_prob)
}
getModelSuccesses <- function (b_mcp) {
stopifnot(inherits(b_mcp, "BayesianMCP"))
model_indices <- grepl("post_probs.", colnames(b_mcp))
model_names <- colnames(b_mcp)[model_indices] |>
sub(pattern = "post_probs.", replacement = "", x = _)
model_successes <- colMeans(b_mcp[, model_indices] > b_mcp[, "crit_prob_adj"])
names(model_successes) <- model_names
return (model_successes)
}
#' @title performBayesianMCPMod
#'
#' @description Performs Bayesian MCP Test step and modeling in a combined fashion. See performBayesianMCP() function for MCP Test step and getModelFits() for the modeling step
#'
#' @param posterior_list An object of class 'postList' or a list of 'postList' objects as created by getPosterior() containing information about the (mixture) posterior distribution per dose group
#' @param contr An object of class 'optContr' as created by the getContr() function. It contains the contrast matrix to be used for the testing step.
#' @param crit_prob_adj A getCritProb object, specifying the critical value to be used for the testing (on the probability scale).
#' @param simple Boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits() function. Default FALSE.
#' @param avg_fit Boolean variable, defining whether an average fit (based on generalized AIC weights) should be performed in addition to the individual models. Default TRUE.
#' @param delta A numeric value for the threshold Delta for the MED assessment. If NULL, no MED assessment is performed. Default NULL.
#' @param evidence_level A numeric value between 0 and 1 for the evidence level gamma for the MED assessment. Only required for Bayesian MED assessment, see ?getMED for details. Default NULL.
#' @param med_selection A string, either "avgFit" or "bestFit" based on the lowest gAIC, for the method of MED selection. Default "avgFit".
#' @param n_samples A numerical for the number of bootstrapped samples in case the Bayesian MED assessment is performed. Default 1e3.
#' @examples
#' mods <- DoseFinding::Mods(linear = NULL,
#' emax = c(0.5, 1.2),
#' exponential = 2,
#' doses = c(0, 0.5, 2,4, 8))
#' dose_levels <- c(0, 0.5, 2, 4, 8)
#' sd_posterior <- c(2.8, 3, 2.5, 3.5, 4)
#' contr_mat <- getContr(
#' mods = mods,
#' dose_levels = dose_levels,
#' cov_posterior = diag(sd_posterior)^2)
#' critVal <- getCritProb(
#' mods = mods,
#' dose_weights = c(50, 50, 50, 50, 50), #reflecting the planned sample size
#' dose_levels = dose_levels,
#' alpha_crit_val = 0.6) # unreasonable alpha for this example, rather choose 0.05
#' prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2),
#' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,
#' DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
#' DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2))
#' mu <- c(0, 1, 1.5, 2, 2.5)
#' S_hat <- diag(c(5, 4, 6, 7, 8)^2)
#' posterior_list <- getPosterior(
#' prior_list = prior_list,
#' mu_hat = mu,
#' S_hat = S_hat,
#' calc_ess = TRUE)
#'
#' performBayesianMCPMod(posterior_list = posterior_list,
#' contr = contr_mat,
#' crit_prob_adj = critVal,
#' simple = FALSE,
#' delta = 1)
#'
#' @return Bayesian MCP test result as well as modeling result.
#'
#' @export
performBayesianMCPMod <- function (
posterior_list,
contr,
crit_prob_adj,
simple = FALSE,
avg_fit = TRUE,
delta = NULL,
evidence_level = NULL,
med_selection = c("avgFit", "bestFit"),
n_samples = 1e3
) {
if (inherits(posterior_list, "postList")) {
posterior_list <- list(posterior_list)
}
checkmate::assert_list(posterior_list, types = "postList")
checkmate::assert_class(contr, "optContr")
checkmate::assert_class(crit_prob_adj, "numeric")
checkmate::assert_logical(simple)
checkmate::assert_logical(avg_fit)
checkmate::assert_double(delta, null.ok = TRUE)
checkmate::assert_double(evidence_level, lower = 0, upper = 1, null.ok = TRUE)
checkmate::assert_character(med_selection)
checkmate::assert_double(n_samples, lower = 1)
med_selection <- match.arg(med_selection, choices = c("avgFit", "bestFit"))
get_med <- ifelse(is.null(delta), FALSE, TRUE)
if (!is.null(evidence_level)) {
stopifnot("delta must not be NULL if evidence_level is not NULL" = !is.null(delta))
}
if (inherits(contr, "optContr")) {
model_names <- colnames(contr$contMat)
dose_levels <- as.numeric(rownames(contr$contMat))
} else if (length(contr) == length(posterior_list)) {
model_names <- colnames(contr[[1]]$contMat)
dose_levels <- as.numeric(rownames(contr[[1]]$contMat))
} else {
stop ("Argument 'contr' must be of type 'optContr'")
}
b_mcp <- performBayesianMCP(
posterior_list = posterior_list,
contr = contr,
crit_prob_adj = crit_prob_adj)
b_mod <- performBayesianMod(
b_mcp = b_mcp,
posterior_list = posterior_list,
model_names = model_names,
dose_levels = dose_levels,
simple = simple,
avg_fit = avg_fit)
b_mcp_mod <- list(BayesianMCP = b_mcp, Mod = b_mod)
if (get_med) {
med_info <- t(sapply(b_mcp_mod$Mod, function (model_fits) {
if (!is.null(model_fits)) {
if (!is.null(evidence_level)) {
bs_quantiles <- getBootstrapQuantiles(
model_fits = model_fits,
quantiles = 1 - evidence_level,
n_samples = n_samples)
med_info <- getMED(
delta = delta,
evidence_level = evidence_level,
bs_quantiles = bs_quantiles)
} else {
med_info <- getMED(
delta = delta,
model_fits = model_fits)
}
if (med_selection == "avgFit") {
return (med_info[, "avgFit"])
} else if (med_selection == "bestFit") {
model_min <- which.min(sapply(model_fits, function (x) x$gAIC))
return (med_info[, model_min])
}
} else {
return (c(med_reached = 0, med = NA))
}
}))
attr(b_mcp_mod, "MED") <- med_info
attr(b_mcp_mod, "MEDSelection") <- med_selection
}
class(b_mcp_mod) <- "BayesianMCPMod"
return (b_mcp_mod)
}
#' @title performBayesianMCP
#'
#' @description Performs Bayesian MCP Test step, as described in Fleischer et al. (2022).
#' Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate the posterior probability is calculated that the contrast is bigger than 0 (based on the posterior distribution of the dose groups).
#' In order to obtain significant test decision we consider the maximum of the posterior probabilities across the different models. This maximum is compared with a (multiplicity adjusted) critical value (on the probability scale).
#' @references Fleischer F, Bossert S, Deng Q, Loley C, Gierse J. 2022. Bayesian MCPMod. Pharmaceutical Statistics. 21(3): 654-670. doi:10.1002/pst.2193
#' @param posterior_list An object derived with getPosterior with information about the (mixture) posterior distribution per dose group
#' @param contr An object of class 'optContr' as created by the getContr() function. It contains the contrast matrix to be used for the testing step.
#' @param crit_prob_adj A getCritProb object, specifying the critical value to be used for the testing (on the probability scale)
#'
#' @examples
#' mods <- DoseFinding::Mods(linear = NULL,
#' emax = c(0.5, 1.2),
#' exponential = 2,
#' doses = c(0, 0.5, 2,4, 8))
#' dose_levels <- c(0, 0.5, 2, 4, 8)
#' sd_posterior <- c(2.8,3,2.5,3.5,4)
#' contr_mat <- getContr(
#' mods = mods,
#' dose_levels = dose_levels,
#' cov_posterior = diag(sd_posterior)^2)
#' critVal <- getCritProb(
#' mods = mods,
#' dose_weights = c(50, 50, 50, 50, 50), #reflecting the planned sample size
#' dose_levels = dose_levels,
#' alpha_crit_val = 0.05)
#' prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2),
#' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,
#' DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
#' DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2))
#' mu <- c(0, 1, 1.5, 2, 2.5)
#' S_hat <- diag(c(5, 4, 6, 7, 8)^2)
#' posterior_list <- getPosterior(
#' prior_list = prior_list,
#' mu_hat = mu,
#' S_hat = S_hat,
#' calc_ess = TRUE)
#'
#' performBayesianMCP(posterior_list = posterior_list,
#' contr = contr_mat,
#' crit_prob_adj = critVal)
#'
#' @return Bayesian MCP test result, with information about p-values for the individual dose-response shapes and overall significance
#'
#' @export
performBayesianMCP <- function(
posterior_list,
contr,
crit_prob_adj
) {
if (inherits(posterior_list, "postList")) {
posterior_list <- list(posterior_list)
}
checkmate::assert_list(posterior_list, types = "postList")
checkmate::assert(
checkmate::check_class(contr, "optContr"),
checkmate::check_list(contr, types = "optContr"),
combine = "or")
checkmate::assert_class(crit_prob_adj, "numeric")
checkmate::assert_numeric(crit_prob_adj, lower = 0, upper = Inf)
## TODO: What case is covered by the IF clause?
## Optionally providing custom contrast?
## This was written before the checkmate assertions were included
## -> if contrasts are re-evaluated, contr becomes a list of contrasts
if (inherits(contr, "optContr")) {
b_mcp <- t(sapply(posterior_list, BayesMCPi, contr, crit_prob_adj))
} else {
b_mcp <- t(mapply(BayesMCPi, posterior_list, contr, crit_prob_adj))
}
class(b_mcp) <- "BayesianMCP"
attr(b_mcp, "critProbAdj") <- crit_prob_adj
attr(b_mcp, "successRate") <- mean(b_mcp[, 1])
attr(b_mcp, "essAvg") <- ifelse(
test = is.na(attr(posterior_list[[1]], "ess")),
yes = numeric(0),
no = rowMeans(sapply(posterior_list,
function (posteriors) attr(posteriors, "ess"))))
return (b_mcp)
}
performBayesianMod <- function (
b_mcp,
posterior_list,
model_names,
dose_levels,
simple,
avg_fit
) {
## Make parallel processing optional
if (requireNamespace("future.apply", quietly = TRUE)) {
optPar_lapply <- future.apply::future_lapply
} else {
optPar_lapply <- lapply
}
fits_list <- optPar_lapply(seq_along(posterior_list), function (i) {
if (b_mcp[i, 1]) {
model_fits <- getModelFits(
models = model_names,
dose_levels = dose_levels,
posterior = posterior_list[[i]],
avg_fit = avg_fit,
simple = simple)
sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "critProbAdj")
model_fits <- addSignificance(model_fits, sign_models)
} else {
NULL
}
})
}
Try the BayesianMCPMod package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.