## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse   = TRUE,
  comment    = "#>",
  message    = FALSE,
  warning    = FALSE,
  fig.width  = 7,
  fig.height = 4.5,
  dpi        = 96
)

## ----libraries----------------------------------------------------------------
library(shrinkr)
library(beastt)
library(distributional)
library(tibble)
library(ggplot2)

set.seed(1104)
sigma_known <- 2   # known within-arm response SD (same in all arms)

## ----hist-data----------------------------------------------------------------
n_hist <- c(58, 43, 62, 51, 47, 55)
hist <- tibble(
  study = paste0("H", 1:6),
  n     = n_hist,
  ybar  = c(9.9, 10.3, 9.7, 10.1, 10.4, 9.8),
  se    = sigma_known / sqrt(n_hist)
)
hist

## ----forest, fig.height = 3.4-------------------------------------------------
ggplot(hist, aes(study, ybar)) +
  geom_hline(yintercept = mean(hist$ybar), linetype = "dashed", color = "grey50") +
  geom_pointrange(aes(ymin = ybar - 1.96 * se, ymax = ybar + 1.96 * se),
                  color = "steelblue", linewidth = 0.8) +
  labs(x = NULL, y = "Control mean", title = "Historical control arms") +
  theme_minimal(base_size = 12)

## ----internal-data------------------------------------------------------------
n_int    <- 70
int_ctrl <- tibble(y = rnorm(n_int, mean = 10, sd = sigma_known))

## ----priors-------------------------------------------------------------------
hierarchical_priors <- list(
  mu  = dist_normal(0, 100),
  tau = dist_truncated(dist_normal(0, sigma_known / 4), lower = 0)
)

## ----prior-predictive---------------------------------------------------------
prior_pred <- sample_prior_predictive(hierarchical_priors,
                                      n_groups = nrow(hist), n_draws = 2000)
plot(prior_pairwise_differences(prior_pred))

## ----fit, results = "hide"----------------------------------------------------
fit_map <- shrink(
  mle                 = hist$ybar,
  var_matrix          = hist$se^2,
  hierarchical_priors = hierarchical_priors,
  chains = 4, iter = 4000, warmup = 1000,
  seed = 2026, refresh = 0, verbose = FALSE
)

## ----fit-summary--------------------------------------------------------------
summarize_mu_tau(fit_map)

## ----make-map-----------------------------------------------------------------
make_map <- function(fit) {
  d <- extract_mu_tau(fit)
  dist_normal(mean(d$mu), sqrt(var(d$mu) + mean(d$tau_squared)))
}

map_prior <- make_map(fit_map)
map_prior

## ----map-ess------------------------------------------------------------------
map_ess <- sigma_known^2 / variance(map_prior)
map_ess

## ----robustify----------------------------------------------------------------
rmp         <- robustify_norm(map_prior, n = map_ess, weights = c(0.5, 0.5))
vague_prior <- dist_normal(mix_means(rmp)[["vague"]], mix_sigmas(rmp)[["vague"]])

plot_dist("MAP (informative)" = map_prior,
          "Vague component"   = vague_prior,
          "Robust mixture"    = rmp)

## ----posterior----------------------------------------------------------------
post_borrow <- calc_post_norm(int_ctrl, response = y,
                              prior = rmp, internal_sd = sigma_known)
post_nobrrw <- calc_post_norm(int_ctrl, response = y,
                              prior = vague_prior, internal_sd = sigma_known)

plot_dist("No borrowing"           = post_nobrrw,
          "Borrowing (robust MAP)" = post_borrow)

## ----ess-table----------------------------------------------------------------
ess_post <- n_int * variance(post_nobrrw) / variance(post_borrow)

tibble(
  quantity     = c("Posterior mean", "Posterior SD", "Effective sample size"),
  no_borrowing = round(c(mean(post_nobrrw), sqrt(variance(post_nobrrw)), n_int), 2),
  robust_map   = round(c(mean(post_borrow), sqrt(variance(post_borrow)), ess_post), 2)
)