## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
suppressPackageStartupMessages({
  library(BayesianMCPMod)
  library(RBesT)
  library(DoseFinding)
  library(dplyr)
})

set.seed(7015)

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)

## -----------------------------------------------------------------------------
trial <- c("trial_1", "trial_2", "trial_3")
n     <- c(70,  115, 147) # sample size per trial
r     <- c( 6,   16,  16) # n responders per trial

## -----------------------------------------------------------------------------
dose_levels <- c(0, 2.5, 5, 10, 20, 50, 100, 200)

# 1) Establish MAP prior (beta mixture distribution)
set.seed(7015) # re-set seed only for this example; remove in your analysis script
map <- gMAP(
  cbind(r, n - r) ~ 1 | trial,
  family     = binomial,
  tau.dist   = "HalfNormal",
  tau.prior  = 0.5,
  beta.prior = (1 / sqrt(0.1 * 0.9)),
  warmup     = 1000,
  iter       = 10000,
  chains     = 2,
  thin       = 1
)
map

prior <- automixfit(map) #fits mixture distribution from MCMC samples from above
p     <- summary(prior)[1]

# 2) Robustify prior
prior_rob <- RBesT::robustify(priormix = prior,
                              mean     = 0.5,
                              weight   = 0.4)

# 3) Translate prior to logit scale (to approximate via normal mixture model)
r                <- rmix(prior_rob, n = 1e4)
log_r            <- RBesT::logit(r)
prior_ctr        <- automixfit(log_r, type = "norm")

# Specification of reference scale (this follows the idea of [@Neuenschwander2016]).
sigma(prior_ctr) <- sqrt(1 / (p * (1 - p)))

# Specify a prior list
prior_trt <- RBesT::mixnorm(
  comp1 = c(
    w = 1,
    m = logit(summary(prior)[1]),
    n = 1
  ),
  sigma = sqrt(1 / (p * (1 - p))),
  param = "mn"
)

prior_list <- c(list(prior_ctr),
                rep(x     = list(prior_trt),
                    times = length(dose_levels[-1])))

dose_names        <- c("Ctr", paste0("DG_", seq_along(dose_levels[-1])))
names(prior_list) <- dose_names

## -----------------------------------------------------------------------------
models <- Mods(
  linear      = NULL,
  sigEmax     = c(50, 3),
  quadratic   = -1 / 250,
  logistic    = c(110, 15),
  exponential = 80,
  emax        = 10,
  doses       = dose_levels,
  placEff     = RBesT::logit(0.118),
  maxEff      = RBesT::logit(0.3) - RBesT::logit(0.118)
)

plot(models)

## -----------------------------------------------------------------------------
data("migraine") # data set from the DoseFinding package

doses_fact <- as.factor(dose_levels)
n_patients <- migraine$ntrt
resp_rate  <- migraine$painfree/n_patients

## Execution of logistic regression and readout of parameters 
## Note that estimates are automatically on the logit scale.
log_fit <- glm(resp_rate ~ doses_fact - 1, family = binomial, weights = n_patients)
mu_hat  <- coef(log_fit)
S_hat   <- vcov(log_fit)

## -----------------------------------------------------------------------------
post_logit <- getPosterior(prior_list, mu_hat = mu_hat, S_hat  = S_hat)

## -----------------------------------------------------------------------------
summary(post_logit, probability_scale = TRUE)

## -----------------------------------------------------------------------------
contr_mat_prior <- getContr(
  mods           = models,
  dose_levels    = dose_levels,
  dose_weights   = n_patients)

set.seed(7015) # re-sets seed only for this example; remove in your analysis script
crit_pval <- getCritProb(
  mods           = models,
  dose_levels    = dose_levels,
  cov_new_trial  = S_hat,
  alpha_crit_val = 0.05
)

## -----------------------------------------------------------------------------
BMCP_result <- performBayesianMCP(
  posterior_list = post_logit,
  contr          = contr_mat_prior, 
  crit_prob_adj  = crit_pval)

## -----------------------------------------------------------------------------
BMCP_result

## -----------------------------------------------------------------------------
model_fits <- getModelFits(
  models            = models,
  dose_levels       = dose_levels,
  posterior         = post_logit,
  simple            = TRUE,
  probability_scale = TRUE)

## -----------------------------------------------------------------------------
plot(model_fits, cr_bands = TRUE)

## -----------------------------------------------------------------------------
plot(model_fits, probability_scale = FALSE)

## -----------------------------------------------------------------------------
display_params_table(stats::predict(model_fits, doses = c(0, 2.5, 10,150, 200)))

## -----------------------------------------------------------------------------
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      = dose_levels,
  n_samples  = 10)

## -----------------------------------------------------------------------------
getMED(
  delta       = 0.16, # on probability scale
  model_fits  = model_fits,
  dose_levels = seq(min(dose_levels), max(dose_levels), by = 1))

## ----eval = FALSE-------------------------------------------------------------
# BMCPMod_result <- performBayesianMCPMod(
#   posterior_list    = post_logit,
#   contr             = contr_mat_prior,
#   crit_prob_adj     = crit_pval,
#   simple            = TRUE,
#   delta             = 0.16,
#   probability_scale = TRUE
# )