---
title: "Robust, Sceptical, and Power Priors"
author: "Ndoh Penn"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
    number_sections: true
vignette: >
  %\VignetteIndexEntry{Robust, Sceptical, and Power Priors}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r 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)
```

# Overview

Beyond the elicited informative prior, regulatory submissions typically require
one or more alternative prior specifications to demonstrate robustness. bayprior
provides three well-established alternatives:

| Prior type | Function | Method | Reference |
|---|---|---|---|
| **MAP (from historical trials)** | `map_prior()` | Random-effects meta-analysis with explicit tau prior | Roever et al. (2021) |
| **Robust mixture** | `robust_prior()` | Mixes informative + vague | Schmidli et al. (2014) |
| **Sceptical** | `sceptical_prior()` | Centred at null effect | Spiegelhalter et al. (1994) |
| **Power prior** | `calibrate_power_prior()` | Down-weights historical data | Ibrahim & Chen (2000) |

---

# Robust Mixture Prior

## Concept

A robust prior protects against prior misspecification by mixing the informative
elicited prior with a vague (diffuse) component:

$$\pi_{\text{robust}}(\theta) = (1 - w) \cdot \pi_{\text{informative}}(\theta)
+ w \cdot \pi_{\text{vague}}(\theta)$$

The vague component — a wide Normal centred at the informative prior mean —
ensures the posterior is never completely dominated by a conflicting prior.
The default weight is $w = 0.20$ (80% informative, 20% vague).

## Usage

```{r 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)
```

```{r 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")
```

## Varying the Vague Weight

A higher vague weight makes the prior more diffuse and reduces its influence
on the posterior:

```{r 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)       
```

## Vague SD Multiplier

The vague component's SD defaults to 10× the informative prior's SD. Adjust
with `vague_sd`:

```{r 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")
```

---

# Deriving the Informative Prior from Historical Trials (MAP)

## Concept

The examples above assume the informative prior fed into `robust_prior()`
already exists (from elicitation, or otherwise). In practice, that informative
component is often itself derived from **multiple historical trials** via a
random-effects meta-analysis. `map_prior()` performs this step directly:
given trial-level effect estimates and standard errors, it returns a
meta-analytic-predictive (MAP) prior suitable for passing straight into
`robust_prior()`.

Unlike treating the between-trial heterogeneity (the parameter $\tau$) as
fixed, `map_prior()` always places an explicit **prior on $\tau$**, following
the weakly-informative recommendations of Roever et al. (2021). This matters
because how much heterogeneity is assumed among historical trials directly
affects how much borrowing is ultimately justified — an assumption that
should be stated and justified explicitly, not left as an unstated default.
Sensible outcome-specific defaults are built in; see `?resolve_tau_prior`.

`map_prior()` fits this model itself, entirely in base R: the posterior of
$\mu$ given $\tau$ has a closed form, so only the marginal posterior of
$\tau$ needs numerical integration (`stats::integrate()`). `bayprior` has no
dependency, direct or optional, on any external meta-analysis package for
this. See `?map_prior` for the underlying model.

## Usage

```{r map-prior}
# Five historical control-arm trials -- i.e. a single-arm log-odds scale,
# not a two-arm log-odds RATIO. This distinction matters: see below.
y  <- c(0.10, -0.05, 0.22, 0.05, 0.15)
se <- c(0.12, 0.15, 0.10, 0.18, 0.14)

map <- map_prior(
  y, se,
  outcome_type = "single_arm_log_odds",
  label        = "Historical control (log-odds)"
)
print(map)
```

`map_prior()` returns an ordinary `bayprior` object (a Normal approximation
to the marginal posterior of the meta-analytic mean, integrating over $\tau$),
so it can be used exactly like any elicited prior:

```{r map-into-robust}
rob_map <- robust_prior(map, vague_weight = 0.20)
plot(rob_map)
```

## Getting `y` and `se` from Raw Trial Data

The example above assumes `y`/`se` are already computed. In practice
they're usually derived from each historical trial's own reported summary
statistics (event counts, arm means/SDs, ...), and
[`metafor::escalc()`](https://cran.r-project.org/package=metafor) is the
standard tool for that conversion. `historical_effect_sizes()` is a thin
wrapper around it that maps `map_prior()`'s `outcome_type` vocabulary onto
the matching `metafor` measure, so raw trial data can be turned into
`y`/`se` in one step:

```{r historical-effect-sizes, eval=requireNamespace("metafor", quietly = TRUE)}
# The same five trials as above, but as originally reported -- event
# counts (xi out of ni) rather than pre-computed log-odds estimates.
trials <- data.frame(
  xi = c(11, 9, 14, 10, 12),
  ni = c(45, 40, 55, 48, 42)
)
es <- historical_effect_sizes("single_arm_log_odds", trials)
es[, c("xi", "ni", "y", "se")]

map_from_counts <- map_prior(
  y = es$y, se = es$se,
  outcome_type = "single_arm_log_odds",
  label        = "Historical control (from raw counts)"
)
print(map_from_counts)
```

This is not the only way to get `y`/`se` -- they can equally be read by
hand off a published point estimate and 95% CI (e.g.
`se = (log(upper) - log(lower)) / (2 * 1.96)` on whatever scale the CI was
reported) -- and `historical_effect_sizes()` is only usable if `metafor` is
installed (it is a `Suggests` dependency, not `Imports`, so it is not
required just to use `map_prior()` itself). See `?historical_effect_sizes`
for the full column-name mapping across outcome types -- two-arm odds
ratios, standardised/raw mean differences, incidence-rate ratios, and
correlations, in addition to the single-arm case shown here.

## Choosing the Heterogeneity Prior Scale

`outcome_type` is not just a label: unless `tau_prior` is supplied
explicitly, it selects the default $\tau$ prior from Roever et al. (2021)'s
outcome-specific recommendations, via `resolve_tau_prior()`:

```{r tau-presets}
resolve_tau_prior("log_or")               # two-arm log-odds ratio
resolve_tau_prior("single_arm_log_odds")  # e.g. a historical control rate
```

Note that these give **different scales (0.5 vs 1.0)**, despite both being
"log-odds" in casual speech -- they are different quantities with different
heterogeneity properties, and conflating them silently applies the wrong
default. `"mean_difference"` (raw, unstandardised units) has no preset at
all -- its scale is entirely endpoint-specific (mmHg, cm, days, ...), so
`map_prior()` requires `tau_prior` to be supplied explicitly for it:

```{r map-tau-sensitivity}
# An explicit tau_prior always overrides the outcome_type preset, and is
# required (map_prior() errors otherwise) for "mean_difference".
map_tight <- map_prior(y, se, tau_prior = list(family = "half_normal", scale = 0.1))
map_wide  <- map_prior(y, se, tau_prior = list(family = "half_normal", scale = 1.0))

cat("Tight tau prior -- posterior SD:", round(map_tight$fit_summary$sd, 4), "\n")
cat("Wide tau prior  -- posterior SD:", round(map_wide$fit_summary$sd, 4), "\n")
```

Allowing more heterogeneity (a wider $\tau$ prior) generally widens the
resulting informative prior, which in turn results in a more conservative
(less confident) prior being carried forward into `robust_prior()`.

## Inspecting the Heterogeneity Posterior

The full posterior distribution of $\tau$ (evaluated on a plotting grid
whose range is chosen automatically from the fitted posterior, not a fixed
guess) is available via `$tau_posterior`, useful for checking whether the
data are actually informative about heterogeneity or whether the posterior
for $\tau$ still closely tracks the prior:

```{r map-tau-posterior}
plot(map$tau_posterior$tau, map$tau_posterior$density, type = "l",
     xlab = expression(tau), ylab = "Posterior density",
     main = "Posterior distribution of the heterogeneity SD")
```

---

# Sceptical Prior

## Concept

The sceptical prior (Spiegelhalter & Freedman, 1994) represents the view of a
conservative regulator who is sceptical of a treatment effect. It is centred at
the null value of the treatment effect with width calibrated to a chosen
strength of scepticism.

This sensitivity prior follows Spiegelhalter and Freedman's (1994) approach
for trials using informative priors: the trial conclusions should hold even
under a prior that places most mass at "no effect."

## Normal Family (Mean Differences, Log Odds Ratios)

For continuous or log-scale quantities where the null is typically 0:

```{r 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")
```

```{r sceptical-normal-plot}
plot(sc_moderate)
```

The SD mapping by strength:

| Strength | SD | Interpretation |
|---|---|---|
| `weak` | 1.0 | Vague scepticism — wide prior around null |
| `moderate` | 0.5 | 2-SD departure from null has ~5% prior probability |
| `strong` | 0.25 | Very concentrated at null — very sceptical |

## Beta Family (Response Rates)

For binary endpoints, `null_value` must be in $(0, 1)$ — it represents the null
response rate, not a difference:

```{r 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)
```

## Log-Normal Family (Hazard Ratios)

For hazard ratios, the null is HR = 1, which corresponds to `null_value = 0`
on the log scale:

```{r 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 vs Sceptical Pair

A common regulatory sensitivity practice is to present conclusions under both
an enthusiastic prior (favouring treatment benefit) and a sceptical prior,
following Spiegelhalter and Freedman (1994):

```{r 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

## Concept

The power prior (Ibrahim & Chen, 2000) provides a principled method for
incorporating historical data by down-weighting it by a factor
$\delta \in (0, 1]$:

$$\pi(\theta | D_0, \delta) \propto L(\theta | D_0)^\delta \cdot \pi_0(\theta)$$

where $D_0$ is the historical data and $\delta$ controls how much weight it
receives. $\delta = 1$ fully incorporates the historical data (standard
Bayesian updating); $\delta \to 0$ ignores it entirely.

## Calibrating $\delta$

`calibrate_power_prior()` selects $\delta$ to achieve a target Bayes Factor
between the historical-data-informed prior and the current likelihood —
ensuring the historical data is incorporated only to the extent it is
compatible with current data:

```{r 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)
```

```{r power-plot}
plot(calib)
```

The calibration curves show:

- **Top panel:** Bayes Factor vs $\delta$. The dashed red line is the target BF;
  the dotted green line marks the optimal $\delta$.
- **Bottom panel:** Box p-value vs $\delta$. Values below 0.05 indicate
  conflict at that weight.

## Compatibility Method

Alternatively, select $\delta$ to be the largest value for which the
historical-data-informed prior shows no conflict with the current data:

```{r 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 Prior for Other Families

Power prior updating is supported for Beta, Normal, Gamma, Log-Normal, and
Mixture priors. For a Normal prior with continuous historical data:

```{r 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)
```

---

# Choosing the Right Alternative Prior

```{r 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")
```

---

# Including Robust Priors in the Regulatory Report

All three robust prior types — robust mixture, sceptical, and power prior —
are automatically included in the downloaded prior justification report when
they have been computed in the session. Pass them directly to `prior_report()`:

```{r 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"
)
```

Each section in the report includes a parameter summary table and the
corresponding density or calibration plot. The compliance checklist in the
report automatically marks "Robust / sceptical prior computed" as Complete
when any of the three types is supplied.

When using the **Shiny app**, the robust priors flow into the report
automatically — simply run the analyses in the Robust Priors panel before
clicking Download Report.

> **Note:** `prior_report()` requires `devtools::install()`, not just
> `devtools::load_all()`. Quarto spawns a fresh R session that requires the
> package to be properly installed.

---

# References

Ibrahim, J. G. & Chen, M.-H. (2000). Power prior distributions for regression
models. *Statistical Science*, 15, 46–60.

Schmidli, H., Gsteiger, S., Roychoudhury, S., O'Hagan, A., Spiegelhalter, D.,
& Neuenschwander, B. (2014). Robust meta-analytic-predictive priors in clinical
trials with historical control information. *Biometrics*, 70, 1023–1032.

Spiegelhalter, D. J., Freedman, L. S., & Parmar, M. K. B. (1994). Bayesian
approaches to randomized trials. *Journal of the Royal Statistical Society A*,
157, 357–416.

Gravestock, I. & Held, L. (2017). Adaptive power priors with empirical Bayes
for clinical trials. *Pharmaceutical Statistics*, 16, 349–360.

Roever, C., Bender, R., Dias, S., Schmid, C. H., Schmidli, H., Sturtz, S.,
Weber, S., & Friede, T. (2021). On weakly informative prior distributions for
the heterogeneity parameter in Bayesian random-effects meta-analysis.
*Research Synthesis Methods*, 12(4), 448–474.