---
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 distribution in the same family as the
informative prior (a uniform `Beta(1, 1)` for a Beta prior, a wide Gamma for a
Gamma prior, a wide Normal for a Normal prior) —
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 Component and its SD

The vague component is chosen in the family of the informative prior, so it
never puts mass outside the parameter space and the mixture stays conjugate.
For a Beta prior the default is the uniform `Beta(1, 1)` (SD 0.289). For a
Normal prior it is a Normal with 10× the informative SD. Adjust the spread
with `vague_sd` (for a Beta prior it must stay below
`sqrt(mean * (1 - mean))`):

```{r robust-sd}
rob_narrow <- robust_prior(informative, vague_weight = 0.20, vague_sd = 0.15)
rob_wide   <- robust_prior(informative, vague_weight = 0.20, vague_sd = 0.35)

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

The 95% credible interval of a mixture is computed from its distribution
function, not as mean ± 1.96 SD:

```{r robust-cri}
round(unlist(rob$fit_summary[c("q025", "q975")]), 3)
```

---

# 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, so it can be used exactly
like any elicited prior. It is the **predictive** distribution of the parameter
of a *new* trial, $\theta^* \sim \mathrm{Normal}(\mu, \tau^2)$, given the
historical data, integrating over both $\mu$ and $\tau$, and matched by moments
to a Normal. Its SD therefore includes the between-trial variance $\tau^2$ as
well as the uncertainty about $\mu$, and is wider than the posterior SD of $\mu$
alone (kept in `map$mu_summary`). The predictive distribution is a scale
mixture of Normals and heavier-tailed than the Normal returned; its exact 2.5%,
50% and 97.5% quantiles are in `map$predictive`.

```{r map-predictive}
c(prior_sd = map$fit_summary$sd, mu_only_sd = map$mu_summary$sd)
unlist(map$predictive[c("q025", "q500", "q975")])
```

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