---
title: "Introduction to bayprior"
author: "Ndoh Penn"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
    number_sections: true
vignette: >
  %\VignetteIndexEntry{Introduction to bayprior}
  %\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
)
```

# Overview

**bayprior** is an R package for principled Bayesian prior elicitation,
conflict detection, and sensitivity analysis in clinical trials. It supports
the documented, systematic prior construction process called for in the
FDA's 2026 draft guidance, *Use of Bayesian Methodology in Clinical Trials
of Drug and Biological Products*, and implements the SHELF elicitation
framework (O'Hagan et al., 2006).

The core analytical workflow is:

1. **Elicit** — Specify a prior from expert knowledge or historical data
2. **Pool** — Combine priors from multiple experts into a consensus prior
3. **Diagnose** — Detect conflict between the prior and observed data
4. **Analyse** — Assess how sensitive posterior conclusions are to prior choice
5. **Robust** — Construct robust, sceptical, or power priors as alternatives
6. **Report** — Generate a regulatory-ready prior justification report

---

# Prior Elicitation

## Beta Prior (Binary Endpoints)

The Beta distribution is the natural conjugate prior for a binary response rate
$\theta \in (0, 1)$. `elicit_beta()` supports two elicitation methods.

### Quantile Matching

The expert specifies quantiles directly — for example: "I believe the true
response rate has a 5% chance of being below 10%, a median of around 30%, and
a 95% chance of being below 60%."

```{r beta-quantile}
library(bayprior)

prior_q <- elicit_beta(
  quantiles = c("0.05" = 0.10, "0.50" = 0.30, "0.95" = 0.60),
  expert_id = "Expert_1",
  label     = "Response rate (treatment arm)"
)
print(prior_q)
```

### Moment Matching

Alternatively, elicit the prior mean and standard deviation directly:

```{r beta-moments}
prior_m <- elicit_beta(
  mean      = 0.30,
  sd        = 0.10,
  method    = "moments",
  expert_id = "Expert_1",
  label     = "Response rate (treatment arm)"
)
plot(prior_m)
```

The shaded region shows the 95% credible interval and the dashed vertical line
marks the prior mean.

## Normal Prior (Continuous Endpoints)

For continuous quantities such as mean differences or log odds ratios:

```{r normal-prior}
prior_nor <- elicit_normal(
  quantiles = c("0.025" = -0.5, "0.50" = 0.20, "0.975" = 0.90),
  label     = "Log odds ratio"
)
print(prior_nor)
```

## Gamma Prior (Rate / Count Endpoints)

For positive-valued quantities such as Poisson event rates or median survival:

```{r gamma-prior}
prior_gam <- elicit_gamma(
  mean   = 5,
  sd     = 2,
  method = "moments",
  label  = "Median OS (months)"
)
plot(prior_gam)
```

## Log-Normal Prior (Hazard Ratios)

For multiplicative quantities that must remain positive:

```{r lognormal-prior}
prior_ln <- elicit_lognormal(
  quantiles = c("0.05" = 0.40, "0.50" = 0.70, "0.95" = 1.20),
  label     = "Hazard ratio (treatment vs control)"
)
plot(prior_ln)
```

## Roulette Method

The roulette method (Oakley & O'Hagan, 2010) converts chip allocations across
pre-defined bins into a parametric prior via quantile matching. It is available
interactively through the Shiny app and programmatically here:

```{r roulette}
prior_rou <- elicit_roulette(
  chips  = c(0L, 2L, 5L, 8L, 5L, 2L, 1L),
  breaks = seq(0, 0.7, by = 0.1),
  family = "beta",
  label  = "Response rate"
)
print(prior_rou)
```

---

# Expert Pooling

When multiple experts contribute, their priors are aggregated into a single
consensus prior via `aggregate_experts()`.

## Linear Pooling

Linear pooling produces a weighted mixture of the individual priors. It is the
most common approach in clinical trial settings (externally Bayesian):

```{r linear-pooling}
p1 <- elicit_beta(mean = 0.25, sd = 0.08, method = "moments",
                  expert_id = "E1", label = "Response rate")
p2 <- elicit_beta(mean = 0.35, sd = 0.10, method = "moments",
                  expert_id = "E2", label = "Response rate")
p3 <- elicit_beta(mean = 0.30, sd = 0.09, method = "moments",
                  expert_id = "E3", label = "Response rate")

consensus <- aggregate_experts(
  priors  = list(E1 = p1, E2 = p2, E3 = p3),
  weights = c(0.40, 0.30, 0.30),
  method  = "linear"
)
print(consensus)
plot(consensus)
```

`aggregate_experts()` automatically computes pairwise Bhattacharyya
coefficients and warns when expert disagreement is substantial.

## Manual Mixture Construction

You can also build a mixture prior directly from components:

```{r manual-mixture}
mix <- elicit_mixture(
  components = list(p1, p2),
  weights    = c(0.5, 0.5),
  label      = "50-50 pooled prior"
)
plot(mix)
```

---

# Prior-Data Conflict Diagnostics

Before updating the prior with trial data it is good practice to check
compatibility between the prior and the observed data. `prior_conflict()`
computes four complementary diagnostics:

| Diagnostic | Interpretation |
|---|---|
| Box p-value | Prior predictive p-value; < 0.05 flags conflict |
| Surprise index | Standardised distance; > 3 indicates high surprise |
| KL divergence | Information distance from prior to likelihood |
| Bhattacharyya overlap | Distributional overlap; < 0.3 is concerning |

```{r conflict}
prior <- elicit_beta(
  mean   = 0.30,
  sd     = 0.10,
  method = "moments",
  label  = "Response rate"
)

# Observed: 18 events in 40 patients
cd <- prior_conflict(
  prior        = prior,
  data_summary = list(type = "binary", x = 18, n = 40),  # also: "continuous", "poisson", "survival"
  alpha        = 0.05
)
print(cd)
```

Visualise the prior-likelihood-posterior overlay to assess compatibility
graphically:

```{r overlay, fig.height = 4}
plot_prior_likelihood(
  prior,
  data_summary   = list(type = "binary", x = 18, n = 40),
  show_posterior = TRUE
)
```

When the prior (blue) and likelihood (orange) overlap substantially there is no
meaningful conflict, and the posterior (green) represents a credible update.

---

# Sensitivity Analysis

Regulatory guidelines require demonstrating that trial conclusions are robust
to plausible prior variations. `sensitivity_grid()` evaluates how posterior
inferences change across a grid of hyperparameter values.

Note: the grids below are intentionally coarse for vignette build speed. In
practice, use finer grids (e.g. `seq(1, 8, 0.5)`) to obtain smooth surfaces.

```{r sensitivity, cache = TRUE}
sa <- sensitivity_grid(
  prior        = prior,
  data_summary = list(type = "binary", x = 14, n = 40),  # also: "continuous", "poisson", "survival"
  param_grid   = list(alpha = seq(1, 6, 1), beta = seq(2, 14, 2)),
  target       = c("posterior_mean", "prob_efficacy"),
  threshold    = 0.30
)
```

## Tornado Plot

The tornado plot ranks posterior quantities by their range across the
hyperparameter grid. Wider bars indicate higher prior sensitivity:

```{r tornado, fig.height = 3}
plot_tornado(sa)
```

## Influence Heatmap

The heatmap shows how the posterior mean changes as both hyperparameters vary.
The orange diamond marks the reference prior:

```{r heatmap, fig.height = 4}
plot_sensitivity(sa, target = "posterior_mean")
```

## Credible Interval Sensitivity

`sensitivity_cri()` focuses specifically on the posterior credible interval
width — a key regulatory quantity:

```{r cri-sensitivity, cache = TRUE, fig.height = 4}
cri_sa <- sensitivity_cri(
  prior        = prior,
  data_summary = list(type = "binary", x = 14, n = 40),  # also: "continuous", "poisson", "survival"
  param_grid   = list(alpha = seq(1, 6, 1), beta = seq(2, 14, 2)),
  cri_level    = 0.95,
  threshold    = 0.30
)
plot_sensitivity(cri_sa, target = "cri_width")
```

---

# Robust and Sceptical Priors

When prior-data conflict is detected, or when a conservative regulatory stance
is required, bayprior provides ready-made alternative prior constructions.

## Robust Mixture Prior

A robust prior mixes an informative component with a vague (diffuse) component
(Schmidli et al., 2014). The vague component ensures the posterior is never
dominated by a conflicting informative prior. The default `vague_weight = 0.20`
gives an 80/20 informative/vague split:

```{r robust, fig.height = 3.5}
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)
```

## Sceptical Prior

A sceptical prior is centred at the null value of the treatment effect and
represents a conservative regulatory stance (Spiegelhalter et al., 1994).

For a **Normal** family (e.g. log odds ratio, mean difference), `null_value`
is the null treatment effect (typically 0):

```{r sceptical-normal, fig.height = 3.5}
sc_norm <- sceptical_prior(
  null_value = 0,
  family     = "normal",
  strength   = "moderate",
  label      = "Log odds ratio (sceptical)"
)
plot(sc_norm)
```

For a **Beta** family (binary response rate), `null_value` must be in (0, 1):

```{r sceptical-beta, fig.height = 3.5}
sc_beta <- sceptical_prior(
  null_value = 0.20,   # null response rate of 20%
  family     = "beta",
  strength   = "moderate",
  label      = "Response rate (sceptical)"
)
plot(sc_beta)
```

---

# Power Prior Calibration

When relevant historical data exist, a power prior down-weights the historical
evidence by $\delta \in (0, 1]$ before incorporating it. `calibrate_power_prior()`
selects $\delta$ to achieve a target Bayes Factor (Ibrahim & Chen, 2000;
Gravestock & Held, 2017).

```{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.10, 1.0, by = 0.10),  # coarse grid for vignette speed
  method          = "bayes_factor"
)
print(calib)
plot(calib)
```

The calibration curves show the Bayes Factor (top) and Box p-value (bottom)
across the $\delta$ grid. The vertical dotted green line marks the optimal
weight and the horizontal dashed lines mark the target BF and $\alpha = 0.05$.

---

# Multivariate Conflict Diagnostics

For trials with co-primary endpoints, `conflict_mahalanobis()` extends the
conflict check to the multivariate setting. The Mahalanobis distance has a
known $\chi^2$ reference distribution under compatibility:

```{r mahal}
mv <- conflict_mahalanobis(
  prior_means = c(0.35, 0.60),
  prior_cov   = matrix(c(0.01, 0.003, 0.003, 0.015), 2, 2),
  obs_means   = c(0.55, 0.58),
  obs_cov     = matrix(c(0.008, 0.002, 0.002, 0.010), 2, 2) / 50,
  labels      = c("Response rate", "OS rate")
)
print(mv)
```

The `marginal_z` scores identify which endpoint contributes most to any
detected conflict.

---

# Generating a Prior Justification Report

Once the analysis is complete, `prior_report()` renders a self-contained
regulatory-ready document in HTML, PDF, or Word format. The report includes
the prior specification, conflict diagnostics, sensitivity visualisations, and
an FDA/EMA compliance checklist.

```{r report, eval = FALSE}
# prior_report() calls rmarkdown::render() internally and cannot be run
# inside a vignette build. Run interactively after loading bayprior.

prior_report(
  prior         = prior,
  conflict      = cd,
  sensitivity   = sa,
  output_format = "html",       # or "pdf" or "docx"
  trial_name    = "TRIAL-001",
  sponsor       = "BioPharma Ltd",
  author        = "J. Smith",
  notes         = paste0(
    "Prior based on Phase 2 internal data and two external expert ",
    "elicitations. Sensitivity analysis confirms robustness across a ",
    "wide range of prior hyperparameter values."
  )
)
```

---

# Complete Worked Example

The following end-to-end example mirrors a realistic binary-endpoint trial
workflow, combining three expert opinions and assessing sensitivity.

```{r full-example, cache = TRUE, fig.height = 4}
# 1. Elicit priors from three experts using different methods
e1 <- elicit_beta(mean = 0.28, sd = 0.08, method = "moments",
                  expert_id = "E1", label = "ORR")
e2 <- elicit_beta(mean = 0.35, sd = 0.10, method = "moments",
                  expert_id = "E2", label = "ORR")
e3 <- elicit_beta(
  quantiles = c("0.10" = 0.15, "0.50" = 0.30, "0.90" = 0.52),
  expert_id = "E3", label = "ORR"
)

# 2. Pool into a consensus prior
consensus <- aggregate_experts(
  priors  = list(E1 = e1, E2 = e2, E3 = e3),
  weights = c(0.40, 0.35, 0.25),
  method  = "linear"
)

# 3. Check conflict with interim data (20 responses in 55 patients)
data_obs <- list(type = "binary", x = 20, n = 55)
cd_full  <- prior_conflict(consensus, data_obs)
print(cd_full)

# 4. Visualise prior-likelihood-posterior
plot_prior_likelihood(consensus, data_obs, show_posterior = TRUE)
```

```{r full-sensitivity, cache = TRUE, fig.height = 3}
# 5. Sensitivity analysis (coarse grid for vignette speed)
sa_full <- sensitivity_grid(
  prior        = consensus,
  data_summary = data_obs,
  param_grid   = list(alpha = seq(1, 6, 1), beta = seq(2, 12, 2)),
  target       = c("posterior_mean", "prob_efficacy"),
  threshold    = 0.25
)
plot_tornado(sa_full)
```

```{r full-heatmap, cache = TRUE, fig.height = 4}
plot_sensitivity(sa_full, target = "prob_efficacy")
```

```{r full-robust, fig.height = 3.5}
# 6. Robust alternative in case conflict worsens at final analysis
rob_full <- robust_prior(consensus, vague_weight = 0.20)
plot(rob_full)
```

---

# Function Reference

| Function | Purpose |
|---|---|
| `elicit_beta()` | Beta prior for binary endpoints — support (0, 1) |
| `elicit_normal()` | Normal prior for continuous endpoints |
| `elicit_gamma()` | Gamma prior for rate/count endpoints |
| `elicit_lognormal()` | Log-normal prior for hazard ratios |
| `elicit_exponential()` | Exponential prior for constant hazard / Poisson rates |
| `elicit_weibull()` | Weibull prior for non-constant hazard survival times |
| `elicit_roulette()` | Chip-based roulette elicitation |
| `elicit_mixture()` | Manual mixture prior construction |
| `aggregate_experts()` | Pool multiple expert priors |
| `prior_conflict()` | Univariate prior-data conflict diagnostics |
| `conflict_mahalanobis()` | Multivariate conflict diagnostics |
| `sensitivity_grid()` | Grid sensitivity of posterior quantities |
| `sensitivity_cri()` | Credible interval sensitivity |
| `plot_sensitivity()` | Heatmap / line plot of sensitivity results |
| `plot_tornado()` | Tornado plot of prior influence |
| `plot_prior_likelihood()` | Prior-likelihood-posterior overlay |
| `robust_prior()` | Robust mixture prior |
| `sceptical_prior()` | Sceptical prior centred at null |
| `calibrate_power_prior()` | Calibrated power prior from historical data |
| `prior_report()` | Regulatory prior justification report |

---

# References

Box, G. E. P. (1980). Sampling and Bayes' inference in scientific modelling and
robustness. *Journal of the Royal Statistical Society A*, 143, 383–430.

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

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

O'Hagan, A., Buck, C. E., Daneshkhah, A., Eiser, J. R., Garthwaite, P. H.,
Jenkinson, D. J., Oakley, J. E., & Rakow, T. (2006). *Uncertain Judgements:
Eliciting Experts' Probabilities*. Wiley.

Oakley, J. E. & O'Hagan, A. (2010). SHELF: the Sheffield Elicitation
Framework. University of Sheffield.

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.