## load packages
library("RBesT")
library("ggplot2")
theme_set(theme_bw(10))

## 4.1. Prior derivation

AS
colSums(AS[, c("n", "r")])

## gMAP example for AS data set
set.seed(35667)
map_mcmc <- gMAP(cbind(r, n-r) ~ 1 | study, family = binomial, data = AS,
                 tau.dist = "HalfNormal", tau.prior = 1, beta.prior = cbind(0, 2))
map_mcmc

## diagnostic plot of model estimate enhanced forest plot
model_plots <- plot(map_mcmc, size = 0.5)
names(model_plots)

## Figure 2
model_plots$forest_model

## 4.2. Parametric mixture prior derivation

## turn MCMC MAP prior into parametric mixture
map <- automixfit(map_mcmc)
map

## diagnostic plot to check how well the EM fitted the MCMC sample
em_diagnostic <- plot(map)

## Figure 3
em_diagnostic$mix + ggtitle(NULL)

em_diagnostic$mixdens
em_diagnostic_link <- plot(map, link = "logit")
em_diagnostic_link$mix

ess(map)

## robustify MAP prior
map_robust <- robustify(map, weight = 0.2, mean = 0.5)
ess(map_robust)

## 4.3. Trial design evaluation

## now evaluate the trial design properties when using the informative MAP prior

## decision function as used in the trial
success <- decision2S(0.95, 0, lower.tail = FALSE)
## an alternative which demands in addition to see at least
## a median difference of at least 10%
success_dual <- decision2S(c(0.95, 0.5), c(0, 0.1), lower.tail = FALSE)

success
success_dual

power.prop.test(p1 = 0.25, p2 = 0.6, power = 0.8, sig.level = 0.05,
                alternative = "one.sided")

## evaluate operating characteristics (the conditional power)

## prior used for treatment in trial
treat_prior <- mixbeta(c(1, 1/2, 1))

## explore different trial designs
design_nonrobust <- oc2S(treat_prior, map, 24, 6, success)
design_robust <- oc2S(treat_prior, map_robust, 24, 6, success)

## calculate point-wise CP
theta  <- c(0.25, 0.5, 0.75)
round(design_nonrobust(theta, theta), 3)
round(design_robust(theta, theta), 3)

round(design_nonrobust(0.25 + 0.35, 0.25), 3)
round(design_robust(0.25 + 0.35, 0.25), 3)

## or better: explore large parameter spaces using graphical plots
curve(design_nonrobust(0.25, 0.25 + x), 0, 0.5, xlab = "Delta",
      ylab = "Power", main = "Power using non-robust (solid),\nrobust (dashed) MAP prior\nand no prior (dotted) for control")
curve(design_robust(0.25, 0.25 + x), add = TRUE, lty = 2)

## 4.4. Probability of success for interim trial analysis

interim_treat <- postmix(treat_prior, r = 5, n = 12)
interim_control <- postmix(map_robust, r = 1, n = 3)

pos_robust <- pos2S(interim_treat, interim_control, 24 - 12, 6 - 3,
                    success)
pos_robust(interim_treat, interim_control)

interim_control_nonrobust <- postmix(map, r = 1, n = 3)
pos_robust(interim_treat, interim_control_nonrobust)

## 4.5. Trial analysis

final_treat <- postmix(treat_prior, r = 15, n = 24)
final_control <- postmix(map_robust, r = 1, n = 6)
success(final_treat, final_control)

pmixdiff(final_treat, final_control, 0, lower.tail = FALSE)

## for more examples and introductory material, please refer to the
## vignettes of the R package https://CRAN.R-project.org/package=RBesT