inst/doc/robust-priors.R

In bayprior: Bayesian Prior Elicitation and Diagnostics for Clinical Trials

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse   = TRUE,
  comment    = "#>",
  fig.width  = 6,
  fig.height = 3.8,
  out.width  = "95%",
  dpi        = 120,
  warning    = FALSE,
  message    = FALSE
)
library(bayprior)

## ----robust-------------------------------------------------------------------
informative <- elicit_beta(
  mean      = 0.30,
  sd        = 0.08,
  method    = "moments",
  label     = "Response rate"
)

rob <- robust_prior(
  informative  = informative,
  vague_weight = 0.20,
  label        = "Robust mixture prior"
)
plot(rob)

## ----robust-summary-----------------------------------------------------------
cat("Informative component weight:", 1 - rob$vague_weight, "\n")
cat("Vague component weight:      ", rob$vague_weight, "\n")
cat("Mixture mean:", round(rob$fit_summary$mean, 4), "\n")
cat("Mixture SD:  ", round(rob$fit_summary$sd,   4), "\n")

## ----robust-weights, fig.height = 4-------------------------------------------
oldpar <- par(mfrow = c(1, 1))
weights <- c(0.10, 0.20, 0.30, 0.50)
cols    <- c("#185FA5", "#1D9E75", "#D85A30", "#888780")

x <- seq(0, 0.8, length.out = 300)
plot(x, bayprior:::.eval_density_vec(informative, x),
     type = "l", lwd = 2, col = "#185FA5",
     xlab = "Response rate", ylab = "Density",
     main = "Effect of vague weight on robust prior",
     ylim = c(0, 6))
for (i in seq_along(weights)) {
  r <- robust_prior(informative, vague_weight = weights[i])
  lines(x, bayprior:::.eval_density_vec(r, x),
        col = cols[i], lwd = 1.5, lty = i + 1)
}
legend("topright",
       legend = c("Informative",
                  paste0("w = ", weights)),
       col    = c("#185FA5", cols),
       lwd    = 2,
       lty    = c(1, 2, 3, 4, 5),
       bty    = "n", cex = 0.8)
par(oldpar)       

## ----robust-sd----------------------------------------------------------------
rob_narrow <- robust_prior(informative, vague_weight = 0.20,
                           vague_sd = 2 * informative$fit_summary$sd)
rob_wide   <- robust_prior(informative, vague_weight = 0.20,
                           vague_sd = 20 * informative$fit_summary$sd)

cat("Narrow vague SD: ", round(rob_narrow$components$vague$fit_summary$sd, 3), "\n")
cat("Default vague SD:", round(rob$components$vague$fit_summary$sd, 3), "\n")
cat("Wide vague SD:   ", round(rob_wide$components$vague$fit_summary$sd, 3), "\n")

## ----sceptical-normal---------------------------------------------------------
sc_weak <- sceptical_prior(
  null_value = 0, family = "normal", strength = "weak",
  label = "Log OR (weak sceptic)"
)
sc_moderate <- sceptical_prior(
  null_value = 0, family = "normal", strength = "moderate",
  label = "Log OR (moderate sceptic)"
)
sc_strong <- sceptical_prior(
  null_value = 0, family = "normal", strength = "strong",
  label = "Log OR (strong sceptic)"
)

cat("Weak SD:    ", sc_weak$fit_summary$sd, "\n")
cat("Moderate SD:", sc_moderate$fit_summary$sd, "\n")
cat("Strong SD:  ", sc_strong$fit_summary$sd, "\n")

## ----sceptical-normal-plot----------------------------------------------------
plot(sc_moderate)

## ----sceptical-beta-----------------------------------------------------------
# Null response rate of 20%: sceptic believes treatment is no better than 20%
sc_beta <- sceptical_prior(
  null_value = 0.20,
  family     = "beta",
  strength   = "moderate",
  label      = "Response rate (sceptical)"
)
plot(sc_beta)

## ----sceptical-lognormal------------------------------------------------------
sc_hr <- sceptical_prior(
  null_value = 0,         # log(1) = 0, i.e. HR = 1
  family     = "lognormal",
  strength   = "moderate",
  label      = "Hazard ratio (sceptical)"
)
plot(sc_hr)

## ----enthusiastic-sceptical---------------------------------------------------
enthusiastic <- elicit_beta(
  mean = 0.45, sd = 0.08,
  method = "moments", label = "Response rate (enthusiastic)"
)
sceptical <- sceptical_prior(
  null_value = 0.20, family = "beta", strength = "moderate",
  label = "Response rate (sceptical)"
)

data_obs <- list(type = "binary", x = 18, n = 40)

post_enth <- bayprior:::.conjugate_update(enthusiastic, data_obs)
post_scep <- bayprior:::.conjugate_update(sceptical,    data_obs)

cat("Posterior mean (enthusiastic):", round(post_enth$fit_summary$mean, 3), "\n")
cat("Posterior mean (sceptical):   ", round(post_scep$fit_summary$mean, 3), "\n")
cat("Posterior SD (enthusiastic):  ", round(post_enth$fit_summary$sd,   3), "\n")
cat("Posterior SD (sceptical):     ", round(post_scep$fit_summary$sd,   3), "\n")

## ----power-prior, cache = TRUE------------------------------------------------
base <- elicit_beta(
  mean   = 0.50,
  sd     = 0.20,
  method = "moments",
  label  = "Response rate"
)

calib <- calibrate_power_prior(
  historical_data = list(type = "binary", x = 12, n = 40),
  current_data    = list(type = "binary", x = 18, n = 50),
  base_prior      = base,
  target_bf       = 3,
  delta_grid      = seq(0.05, 1.0, by = 0.05),
  method          = "bayes_factor"
)
print(calib)

## ----power-plot---------------------------------------------------------------
plot(calib)

## ----power-compat, cache = TRUE-----------------------------------------------
calib_compat <- calibrate_power_prior(
  historical_data = list(type = "binary", x = 12, n = 40),
  current_data    = list(type = "binary", x = 18, n = 50),
  base_prior      = base,
  method          = "compatibility",
  delta_grid      = seq(0.05, 1.0, by = 0.05)
)
cat("Optimal delta (BF method):           ", calib$delta_opt, "\n")
cat("Optimal delta (compatibility method):", calib_compat$delta_opt, "\n")

## ----power-normal, cache = TRUE-----------------------------------------------
base_norm <- elicit_normal(
  mean = 0.0, sd = 0.5,
  method = "moments", label = "Mean difference"
)

calib_norm <- calibrate_power_prior(
  historical_data = list(type = "continuous", x = 0.35, sd = 0.3, n = 60),
  current_data    = list(type = "continuous", x = 0.42, sd = 0.3, n = 80),
  base_prior      = base_norm,
  target_bf       = 3,
  delta_grid      = seq(0.05, 1.0, by = 0.10),
  method          = "bayes_factor"
)
print(calib_norm)

## ----choice-table-------------------------------------------------------------
library(knitr)
kable(data.frame(
  Situation = c(
    "No conflict, regulatory requirement",
    "Mild conflict detected",
    "Severe conflict detected",
    "Historical data available",
    "Enthusiastic/sceptical pair for sensitivity analysis"
  ),
  `Recommended prior` = c(
    "Robust mixture (w = 0.20)",
    "Robust mixture (w = 0.30-0.40)",
    "Sceptical prior (moderate-strong)",
    "Power prior (calibrated)",
    "Sceptical prior as second arm"
  ),
  check.names = FALSE
), align = "ll")

## ----report-integration, eval = FALSE-----------------------------------------
# # After running the analyses above...
# prior_report(
#   prior           = prior,
#   conflict        = cd,
#   sensitivity     = sa,
#   robust_prior    = rob,        # adds "Robust Mixture" section to report
#   sceptical_prior = scep,       # adds "Sceptical Prior" section to report
#   power_prior     = calib,      # adds "Power Prior" section with calibration table
#   output_format   = "html",
#   output_file     = "prior_justification_report",
#   trial_name      = "TRIAL-001",
#   sponsor         = "BioPharma Ltd",
#   author          = "J. Smith, Biostatistician"
# )