inst/doc/analysis_normal.R

In BayesianMCPMod: Simulate, Evaluate, and Analyze Dose Finding Trials with Bayesian MCPMod

## ----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)