## ----setup-------------------------------------------------------------------- library(BayesianMCPMod) library(RBesT) library(clinDR) library(dplyr) library(tibble) library(reactable) set.seed(7015) ## ----table_function, collapse=TRUE-------------------------------------------- display_params_table <- function(named_list) { round_numeric <- function(x, digits = 3) if (is.numeric(x)) round(x, digits) else x tbl <- data.frame( Name = names(named_list), Value = I(lapply(named_list, function(v) { if (inherits(v, "Date")) v <- as.character(v) if (!is.null(names(v))) paste0("{", paste(names(v), v, sep="=", collapse=", "), "}") else v })) ) tbl$Value <- lapply(tbl$Value, round_numeric) knitr::kable(tbl) } ## ----eval = FALSE------------------------------------------------------------- # future::plan(future::multisession, workers = 4L) ## ----Historical Data for Control Arm------------------------------------------ data("metaData") dataset <- filter(as.data.frame(metaData), bname == "BRINTELLIX") histcontrol <- filter( dataset, dose == 0, primtime == 8, indication == "MAJOR DEPRESSIVE DISORDER", protid != 5) hist_data <- data.frame( trial = histcontrol$nctno, est = histcontrol$rslt, se = histcontrol$se, sd = histcontrol$sd, n = histcontrol$sampsize) ## ----Defining MAP prior function---------------------------------------------- getPriorList <- function ( hist_data, dose_levels, dose_names = NULL, robust_weight = 0.5 ) { sd_tot <- with(hist_data, sum(sd * n) / sum(n)) gmap <- RBesT::gMAP( formula = cbind(est, se) ~ 1 | trial, weights = hist_data$n, data = hist_data, family = gaussian, beta.prior = cbind(0, 100 * sd_tot), tau.dist = "HalfNormal", tau.prior = cbind(0, sd_tot / 4)) prior_ctr <- RBesT::automixfit(gmap) if (!is.null(robust_weight)) { prior_ctr <- suppressMessages(RBesT::robustify( priormix = prior_ctr, weight = robust_weight, sigma = sd_tot)) } prior_trt <- RBesT::mixnorm( comp1 = c(w = 1, m = summary(prior_ctr)[1], n = 1), sigma = sd_tot, param = "mn") prior_list <- c(list(prior_ctr), rep(x = list(prior_trt), times = length(dose_levels[-1]))) if (is.null(dose_names)) { dose_names <- c("Ctr", paste0("DG_", seq_along(dose_levels[-1]))) } names(prior_list) <- dose_names return (prior_list) } ## ----Getting the MAP prior---------------------------------------------------- dose_levels <- c(0, 2.5, 5, 10) set.seed(7015) # re-sets seed only for this example; remove in your analysis script prior_list <- getPriorList( hist_data = hist_data, dose_levels = dose_levels, robust_weight = 0.3) getESS(prior_list) ## ----------------------------------------------------------------------------- # Guesstimate estimation exp_guesst <- DoseFinding::guesst( model = "exponential", d = 5, p = 0.2, Maxd = max(dose_levels) ) emax_guesst <- DoseFinding::guesst( model = "emax", d = 2.5, p = 0.9 ) sigEmax_guesst <- DoseFinding::guesst( model = "sigEmax", d = c(2.5, 5), p = c(0.5, 0.95) ) logistic_guesst <- DoseFinding::guesst( model = "logistic", d = c(5, 10), p = c(0.1, 0.85) ) ## ----------------------------------------------------------------------------- betaMod_params <- c(delta1 = 1, delta2 = 1) quadratic_params <- c(delta2 = -0.1) ## ----------------------------------------------------------------------------- mods <- DoseFinding::Mods( linear = NULL, # guesstimate scale exponential = exp_guesst, emax = emax_guesst, sigEmax = sigEmax_guesst, logistic = logistic_guesst, # parameter scale betaMod = betaMod_params, quadratic = quadratic_params, # Options for all models doses = dose_levels, maxEff = -1, placEff = -12.8 ) plot(mods) ## ----------------------------------------------------------------------------- display_params_table(mods) ## ----------------------------------------------------------------------------- knitr::kable(DoseFinding::getResp(mods, doses = dose_levels)) ## ----------------------------------------------------------------------------- data("metaData") trial_data <- dplyr::filter( dplyr::filter(tibble::tibble(metaData), bname == "BRINTELLIX"), primtime == 8, indication == "MAJOR DEPRESSIVE DISORDER", protid == 5 ) n_patients <- c(128, 124, 129, 122) ## ----------------------------------------------------------------------------- posterior <- getPosterior( prior_list = prior_list, mu_hat = trial_data$rslt, S_hat = diag(trial_data$se^2), calc_ess = TRUE ) knitr::kable(summary(posterior)) ## ----------------------------------------------------------------------------- set.seed(7015) # re-sets seed only for this example; remove in your analysis script crit_pval <- getCritProb( mods = mods, dose_levels = dose_levels, cov_new_trial = diag(trial_data$se^2), alpha_crit_val = 0.05 ) contr_mat <- getContr( mods = mods, dose_levels = dose_levels, cov_posterior = diag(summary(posterior)[, 2]^2) ) ## ----------------------------------------------------------------------------- # # i) the frequentist contrast # contr_mat_prior <- getContr( # mods = mods, # dose_levels = dose_levels, # dose_weights = n_patients, # prior_list = prior_list) # # ii) re-estimated frequentist contrasts # contr_mat_prior <- getContr( # mods = mods, # dose_levels = dose_levels, # cov_new_trial = diag(trial_data$se^2)) # # iii) Bayesian approach using number of patients for new trial and prior distribution # contr_mat_prior <- getContr( # mods = mods, # dose_levels = dose_levels, # dose_weights = n_patients, # prior_list = prior_list) ## ----------------------------------------------------------------------------- BMCP_result <- performBayesianMCP( posterior_list = posterior, contr = contr_mat, crit_prob_adj = crit_pval) ## ----------------------------------------------------------------------------- BMCP_result ## ----------------------------------------------------------------------------- # If simple = TRUE, uses approx posterior # Here we use complete posterior distribution model_fits <- getModelFits( models = mods, dose_levels = dose_levels, posterior = posterior, simple = FALSE) ## ----------------------------------------------------------------------------- display_params_table(stats::predict(model_fits, doses = c(0, 2.5, 4, 5, 7, 10))) ## ----------------------------------------------------------------------------- plot(model_fits) ## ----------------------------------------------------------------------------- plot(model_fits, cr_bands = TRUE) ## ----------------------------------------------------------------------------- set.seed(7015) # re-sets seed only for this example; remove in your analysis script bootstrap_quantiles <- getBootstrapQuantiles( model_fits = model_fits, quantiles = c(0.025, 0.5, 0.975), doses = c(0, 2.5, 4, 5, 7, 10), n_samples = 10) ## ----collapse = TRUE---------------------------------------------------------- reactable::reactable( data = bootstrap_quantiles |> tidyr::pivot_wider(names_from = q_prob, values_from = q_val), groupBy = "model", columns = list( dose = colDef(format = list(aggregated = colFormat(suffix = " dose"))), "0.025" = colDef(format = list(cell = colFormat(digits = 4))), "0.5" = colDef(format = list(cell = colFormat(digits = 4))), "0.975" = colDef(format = list(cell = colFormat(digits = 4))) ) ) ## ----------------------------------------------------------------------------- getMED( delta = 4, model_fits = model_fits, dose_levels = seq(min(dose_levels), max(dose_levels), by = 0.01)) ## ----------------------------------------------------------------------------- # performBayesianMCPMod( # posterior_list = posterior, # contr = contr_mat, # crit_prob_adj = crit_pval, # delta = 4, # simple = FALSE)