## ======================================================================
## GBT-PP: Generalized Bayesian TOST with Power-Discounted Priors
## Script faithful to Monleon-Getino (2026), Sections 4-7 and
## Appendices A-C, plus illustrative plots (Section 6 sensitivity
## analysis and TOST-style posterior visualization).
##
## Two gaps in the Appendix B listing are completed here because the
## article's text (Section 7.2) explicitly requires them to reproduce
## Table 2, even though the printed code did not include them:
##   (1) the simple Normal approximation to Peq (gamma_theta = 0) in the
##       frequency case, alongside Monte Carlo and Edgeworth;
##   (2) the dual decision rule (Eq. 11) and the sample-size formula
##       (Corollary 5, Eqs. 17-18), which the text formalizes but does
##       not render as executable code.
## Everything else follows the Appendix B listing as written.
## ======================================================================

## ----------------------------------------------------------------------
## SECTION 4.2 / 4.4 / APPENDIX A.1
## Continuous case: Normal-Normal prior with power discounting
## Eqs. (4)-(9)
## ----------------------------------------------------------------------
gbtpp_continuous <- function(x0_bar, s0, n0, x_bar, s, n, w, delta,
                             sigma0 = NULL, sigma0_max_mult = 10) {
  
  if (is.null(sigma0)) {
    sigma0 <- sigma0_max_mult * s0 / sqrt(n0)
  }
  
  n0_eff    <- w * n0
  prec_hist <- ifelse(n0_eff > 0, n0_eff / s0^2, 0)
  prec_cur  <- n / s^2
  
  tau2    <- 1 / (prec_hist + prec_cur)
  theta_n <- tau2 * (prec_hist * x0_bar + prec_cur * x_bar)
  tau_n   <- sqrt(tau2)
  
  Peq <- pnorm((delta - theta_n) / tau_n) - pnorm((-delta - theta_n) / tau_n)
  
  pi_eq <- pnorm(delta / sigma0) - pnorm(-delta / sigma0)
  
  post_odds  <- Peq / (1 - Peq)
  prior_odds <- pi_eq / (1 - pi_eq)
  BF01 <- post_odds / prior_odds
  
  list(theta_n = theta_n, tau_n = tau_n, Peq = Peq,
       pi_eq = pi_eq, BF01 = BF01, sigma0 = sigma0)
}

## ----------------------------------------------------------------------
## SECTION 4.3 / APPENDIX A.2
## Frequency case: Beta-Binomial prior with power discounting
## Eqs. (6), (10), (20)-(21), (26)
## ----------------------------------------------------------------------
gbtpp_frequency <- function(s0, n0, s, n, w, a0 = 1, b0 = 1,
                            delta, a2 = NULL, b2 = NULL,
                            M = 2e5, seed = 123, return_draws = FALSE) {
  
  f0 <- n0 - s0
  f  <- n  - s
  
  a_ref <- a0 + w * s0 + s
  b_ref <- b0 + w * f0 + f
  
  if (is.null(a2) || is.null(b2)) {
    stop("Provide control-arm posterior parameters a2, b2 for the ",
         "two-proportion contrast theta = p1 - p2 (Section 4.3).")
  }
  
  beta_moments <- function(a, b) {
    mu   <- a / (a + b)
    sig2 <- a * b / ((a + b)^2 * (a + b + 1))
    gam1 <- 2 * (b - a) * sqrt(a + b + 1) / ((a + b + 2) * sqrt(a * b))
    c(mu = mu, sig2 = sig2, gam1 = gam1)
  }
  m1 <- beta_moments(a_ref, b_ref)
  m2 <- beta_moments(a2, b2)
  
  mu_theta  <- unname(m1["mu"] - m2["mu"])
  sig_theta <- unname(sqrt(m1["sig2"] + m2["sig2"]))
  gam_theta <- unname((m1["gam1"] * m1["sig2"]^1.5 -
                         m2["gam1"] * m2["sig2"]^1.5) / sig_theta^3)
  
  zplus  <- (delta  - mu_theta) / sig_theta
  zminus <- (-delta - mu_theta) / sig_theta
  
  set.seed(seed)
  p1 <- rbeta(M, a_ref, b_ref)
  p2 <- rbeta(M, a2, b2)
  theta  <- p1 - p2
  Peq_MC <- mean(theta > -delta & theta < delta)
  SE_MC  <- sqrt(Peq_MC * (1 - Peq_MC) / M)
  
  edge_cdf <- function(z) pnorm(z) - (gam_theta / 6) * dnorm(z) * (z^2 - 1)
  Peq_Edge <- edge_cdf(zplus) - edge_cdf(zminus)
  
  Peq_Normal <- pnorm(zplus) - pnorm(zminus)
  
  out <- list(a_ref = a_ref, b_ref = b_ref,
              mu_theta = mu_theta, sig_theta = sig_theta, gam_theta = gam_theta,
              Peq_MC = Peq_MC, SE_MC = SE_MC,
              Peq_Edge = Peq_Edge, Peq_Normal = Peq_Normal)
  if (return_draws) out$theta_draws <- theta
  out
}

## ----------------------------------------------------------------------
## SECTION 5 / Eq. (11): dual decision rule
## ----------------------------------------------------------------------
gbtpp_decision <- function(Peq, BF01, gamma = 0.95, BF_star = 3) {
  decision <- (Peq >= gamma) && (BF01 >= BF_star)
  list(decision = decision, Peq_pass = Peq >= gamma,
       BF_pass = BF01 >= BF_star)
}

## ----------------------------------------------------------------------
## SECTION 5 / Corollary 5, Eqs. (17)-(18): approximate sample size
## ----------------------------------------------------------------------
gbtpp_sample_size <- function(delta, sigma, w, n0, gamma = 0.95, s0 = sigma) {
  z_gamma <- qnorm((1 + gamma) / 2)
  prec_hist <- (w * n0) / s0^2
  prec_needed <- (z_gamma / delta)^2
  prec_cur_needed <- prec_needed - prec_hist
  if (prec_cur_needed <= 0) {
    return(list(n_required = 0,
                note = "Historical data (discounted by w) already suffice."))
  }
  n_required <- ceiling(prec_cur_needed * sigma^2)
  list(n_required = n_required, prec_hist = prec_hist,
       prec_needed = prec_needed)
}

## ----------------------------------------------------------------------
## SECTION 6: sensitivity analysis over a grid of discount weights w
## ----------------------------------------------------------------------
gbtpp_sensitivity_continuous <- function(w_grid = c(0, 0.25, 0.5, 0.75, 1), ...) {
  res <- lapply(w_grid, function(w) {
    out <- gbtpp_continuous(w = w, ...)
    data.frame(w = w, theta_n = out$theta_n, tau_n = out$tau_n,
               Peq = out$Peq, pi_eq = out$pi_eq, BF01 = out$BF01)
  })
  do.call(rbind, res)
}

gbtpp_sensitivity_frequency <- function(w_grid = c(0, 0.25, 0.5, 0.75, 1), ...) {
  res <- lapply(w_grid, function(w) {
    out <- gbtpp_frequency(w = w, ...)
    data.frame(w = w, a_ref = out$a_ref, b_ref = out$b_ref,
               Peq_MC = out$Peq_MC, SE_MC = out$SE_MC,
               Peq_Edge = out$Peq_Edge, Peq_Normal = out$Peq_Normal)
  })
  do.call(rbind, res)
}

## ======================================================================
## SECTION 7.1 - Worked example: continuous case (bioequivalence)
## Reproduces Table 1 of the article.
## ======================================================================
sigma0_fixed <- 10 * 0.15 / sqrt(24)

table1 <- gbtpp_sensitivity_continuous(
  w_grid  = c(0, 0.25, 0.5, 0.75, 1),
  x0_bar  = 0.02, s0 = 0.15, n0 = 24,
  x_bar   = 0.03, s  = 0.12, n  = 30,
  delta   = 0.05,
  sigma0  = sigma0_fixed
)
cat("\n--- Table 1: continuous case (Section 7.1) ---\n")
print(table1, digits = 4)

## ======================================================================
## SECTION 7.2 - Worked example: frequency case (cure rates)
## Reproduces Table 2 of the article.
## ======================================================================
table2 <- gbtpp_sensitivity_frequency(
  w_grid = c(0, 0.25, 0.5, 0.75, 1),
  s0 = 45, n0 = 50,
  s  = 42, n  = 50,
  a0 = 1,  b0 = 1,
  delta = 0.10,
  a2 = 38, b2 = 14,
  M  = 2e5, seed = 123
)
cat("\n--- Table 2: frequency case (Section 7.2) ---\n")
print(table2, digits = 4)

## ======================================================================
## Dual decision rule example (w = 1, continuous case)
## ======================================================================
res_w1 <- gbtpp_continuous(x0_bar = 0.02, s0 = 0.15, n0 = 24,
                           x_bar = 0.03, s = 0.12, n = 30,
                           w = 1, delta = 0.05, sigma0 = sigma0_fixed)
decision_w1 <- gbtpp_decision(res_w1$Peq, res_w1$BF01)
cat("\nDual decision at w = 1 (continuous case):\n")
print(decision_w1)

## ======================================================================
## Sample-size example (Corollary 5) at w = 0.5
## ======================================================================
n_needed <- gbtpp_sample_size(delta = 0.05, sigma = 0.12, w = 0.5,
                              n0 = 24, gamma = 0.95, s0 = 0.15)
cat("\nRequired current-study n at w = 0.5, gamma = 0.95:\n")
print(n_needed)


## ======================================================================
## ======================================================================
##                              PLOTS
## Four figures, each illustrating a key result of the article.
## All use base R graphics only (no extra packages required).
## ======================================================================
## ======================================================================

## Colorblind-friendly palette (Wong, 2011)
col_blue   <- "#0072B2"
col_orange <- "#D55E00"
col_green  <- "#009E73"
col_grey   <- "grey50"

## ------------------------------------------------------------------
## FIGURE 1. Sensitivity curves for the continuous case (Section 6)
##
## WHAT IT SHOWS: how the posterior equivalence probability Peq and the
## Bayes factor BF01 change as the discount weight w increases from 0
## (fully discounted / ignore historical data) to 1 (fully pooled /
## historical data treated as equally informative as current data).
##
## WHY IT MATTERS: this is the core robustness check proposed in
## Section 6 of the article. A flat or smoothly increasing Peq curve
## with w indicates that the equivalence conclusion is NOT an artifact
## of over-borrowing from historical data (Remark 3). The horizontal
## reference lines mark the decision thresholds gamma = 0.95 (for Peq)
## and BF* = 3 (for BF01), as used in the dual decision rule, Eq. (11).
## ------------------------------------------------------------------
plot_sensitivity_continuous <- function(table1, gamma = 0.95, BF_star = 3) {
  op <- par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3, 1))
  
  ## Panel A: Peq vs w
  plot(table1$w, table1$Peq, type = "b", pch = 19, lwd = 2,
       col = col_blue, ylim = c(0, 1),
       xlab = "Discount weight (w)",
       ylab = expression(P[eq]~"(posterior equivalence probability)"),
       main = "A. Equivalence probability vs. w")
  abline(h = gamma, col = col_grey, lty = 2, lwd = 1.5)
  text(0.02, gamma + 0.03, paste0("gamma = ", gamma), col = col_grey,
       adj = 0, cex = 0.85)
  grid(col = "grey90")
  
  ## Panel B: BF01 vs w (log scale, since BF can span orders of magnitude)
  plot(table1$w, table1$BF01, type = "b", pch = 19, lwd = 2,
       col = col_orange, log = "y",
       xlab = "Discount weight (w)",
       ylab = expression(BF[0][1]~"(Bayes factor for H0: equivalence)"),
       main = "B. Bayes factor vs. w")
  abline(h = BF_star, col = col_grey, lty = 2, lwd = 1.5)
  text(0.02, BF_star * 1.15, paste0("BF* = ", BF_star), col = col_grey,
       adj = 0, cex = 0.85)
  grid(col = "grey90")
  
  par(op)
  mtext("Figure 1: Sensitivity analysis, continuous case (Section 6, Table 1)",
        side = 3, line = -1.5, outer = TRUE, cex = 0.9, font = 2)
}

## ------------------------------------------------------------------
## FIGURE 2. TOST-style posterior visualization (continuous case)
##
## WHAT IT SHOWS: the posterior distribution of theta (the parameter
## of interest, e.g. a mean difference) under a chosen w, together with
## the equivalence region E = (-delta, delta). The shaded area under
## the curve within E is exactly Peq.
##
## WHY IT MATTERS: this is the direct graphical counterpart of the
## generalized Bayesian TOST (Eq. 1 and Section 3). It makes visible,
## in a single picture, why Peq is high (most posterior mass inside E)
## or low (posterior concentrated outside E, or straddling the
## boundary), which a table of numbers alone does not convey.
## ------------------------------------------------------------------
plot_posterior_continuous <- function(theta_n, tau_n, delta, Peq, w,
                                      xlim_mult = 5) {
  xs <- seq(theta_n - xlim_mult * tau_n, theta_n + xlim_mult * tau_n,
            length.out = 1000)
  ys <- dnorm(xs, theta_n, tau_n)
  
  plot(xs, ys, type = "n",
       xlab = expression(theta~"(parameter of interest)"),
       ylab = "Posterior density",
       main = paste0("Figure 2: Posterior distribution and equivalence ",
                     "region (w = ", w, ")"))
  
  ## Shade the equivalence region E = (-delta, delta)
  x_in <- xs[xs > -delta & xs < delta]
  y_in <- dnorm(x_in, theta_n, tau_n)
  polygon(c(-delta, x_in, delta), c(0, y_in, 0),
          col = adjustcolor(col_blue, alpha.f = 0.35), border = NA)
  
  ## Full posterior curve
  lines(xs, ys, col = col_blue, lwd = 2)
  
  ## Equivalence margins and posterior mean
  abline(v = c(-delta, delta), col = col_orange, lty = 2, lwd = 1.5)
  abline(v = theta_n, col = col_grey, lty = 3, lwd = 1.5)
  
  legend("topright", bty = "n", cex = 0.85,
         legend = c(expression("Posterior of"~theta),
                    expression("Equivalence margins"~-delta*", "*delta),
                    "Posterior mean",
                    paste0("Peq = ", round(Peq, 3))),
         col = c(col_blue, col_orange, col_grey, NA),
         lty = c(1, 2, 3, NA), lwd = c(2, 1.5, 1.5, NA))
}

## ------------------------------------------------------------------
## FIGURE 3. Sensitivity curves for the frequency case, comparing the
## three approximation methods to Peq (Section 7.2 / Proposition 6)
##
## WHAT IT SHOWS: Peq computed via Monte Carlo (gold standard), the
## Edgeworth approximation, and the simple Normal approximation, all
## plotted together across the discount-weight grid.
##
## WHY IT MATTERS: this directly validates Proposition 6 of the
## article: the Edgeworth correction should track the Monte Carlo curve
## much more closely than the simple Normal approximation, especially
## when the Beta posteriors are skewed (small/large a_ref, b_ref). If
## the three lines are visually indistinguishable, the Normal
## approximation is adequate for this dataset; if they diverge, the
## skewness correction (gamma_theta) matters.
## ------------------------------------------------------------------
plot_sensitivity_frequency <- function(table2) {
  plot(table2$w, table2$Peq_MC, type = "b", pch = 19, lwd = 2,
       col = col_blue, ylim = range(c(table2$Peq_MC, table2$Peq_Edge,
                                      table2$Peq_Normal)),
       xlab = "Discount weight (w)",
       ylab = expression(P[eq]),
       main = "Figure 3: Peq by approximation method, frequency case")
  lines(table2$w, table2$Peq_Edge, type = "b", pch = 17, lwd = 2,
        col = col_orange, lty = 2)
  lines(table2$w, table2$Peq_Normal, type = "b", pch = 15, lwd = 2,
        col = col_green, lty = 3)
  grid(col = "grey90")
  legend("bottomleft", bty = "n", cex = 0.85,
         legend = c("Monte Carlo (gold standard, Eq. 26)",
                    "Edgeworth approximation (Prop. 6)",
                    "Simple Normal approximation"),
         col = c(col_blue, col_orange, col_green),
         pch = c(19, 17, 15), lty = c(1, 2, 3), lwd = 2)
}

## ------------------------------------------------------------------
## FIGURE 4. TOST-style posterior visualization (frequency case)
##
## WHAT IT SHOWS: the Monte Carlo-simulated posterior distribution of
## theta = p1 - p2 (risk difference) at a chosen w, with the
## equivalence region E = (-delta, delta) shaded.
##
## WHY IT MATTERS: analogous to Figure 2 but for the Beta-Binomial
## model; it is the empirical (simulation-based) counterpart used to
## sanity-check the analytic Edgeworth/Normal approximations plotted
## in Figure 3.
## ------------------------------------------------------------------
plot_posterior_frequency <- function(theta_draws, delta, Peq_MC, w) {
  dens <- density(theta_draws)
  
  plot(dens, col = col_blue, lwd = 2,
       xlab = expression(theta == p[1] - p[2]~"(risk difference)"),
       ylab = "Posterior density (Monte Carlo)",
       main = paste0("Figure 4: Simulated posterior of ",
                     "theta = p1 - p2 (w = ", w, ")"))
  
  x_in <- dens$x[dens$x > -delta & dens$x < delta]
  y_in <- dens$y[dens$x > -delta & dens$x < delta]
  polygon(c(-delta, x_in, delta), c(0, y_in, 0),
          col = adjustcolor(col_blue, alpha.f = 0.35), border = NA)
  
  abline(v = c(-delta, delta), col = col_orange, lty = 2, lwd = 1.5)
  abline(v = mean(theta_draws), col = col_grey, lty = 3, lwd = 1.5)
  
  legend("topright", bty = "n", cex = 0.85,
         legend = c("Monte Carlo posterior density",
                    expression("Equivalence margins"~-delta*", "*delta),
                    "Posterior mean",
                    paste0("Peq (MC) = ", round(Peq_MC, 3))),
         col = c(col_blue, col_orange, col_grey, NA),
         lty = c(1, 2, 3, NA), lwd = c(2, 1.5, 1.5, NA))
}

## ======================================================================
## GENERATE THE FOUR FIGURES
## ======================================================================

## Figure 1: sensitivity curves, continuous case
plot_sensitivity_continuous(table1, gamma = 0.95, BF_star = 3)

## Figure 2: posterior + equivalence region at w = 1 (continuous case)
#dev.new()  # opens a new plotting window; remove this line if running
# in RStudio / RMarkdown where each plot call is shown separately
plot_posterior_continuous(theta_n = res_w1$theta_n, tau_n = res_w1$tau_n,
                          delta = 0.05, Peq = res_w1$Peq, w = 1)

## Figure 3: sensitivity curves across approximation methods, frequency case
#dev.new()
plot_sensitivity_frequency(table2)

## Figure 4: Monte Carlo posterior + equivalence region (frequency case, w = 1)
res_freq_w1 <- gbtpp_frequency(s0 = 45, n0 = 50, s = 42, n = 50,
                               w = 1, a0 = 1, b0 = 1, delta = 0.10,
                               a2 = 38, b2 = 14, M = 2e5, seed = 123,
                               return_draws = TRUE)
#dev.new()
plot_posterior_frequency(theta_draws = res_freq_w1$theta_draws,
                         delta = 0.10, Peq_MC = res_freq_w1$Peq_MC, w = 1)