---
title: "Prior-Data Conflict Diagnostics"
author: "Ndoh Penn"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Prior-Data Conflict Diagnostics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 7, fig.height = 4,
                      message = FALSE, warning = FALSE)
library(bayprior)
```

# Why Conflict Diagnostics Matter

Before updating a prior with trial data it is essential to check whether the
two are compatible. Substantial prior-data conflict indicates one of three
things: the prior was mis-specified, the data are anomalous, or the model is
wrong. Regulators increasingly require evidence that conflict was assessed and
addressed.

**bayprior** provides five complementary univariate diagnostics (one of them
an exploratory companion to the Box p-value) and one
multivariate diagnostic, following Box (1980) and related work. The conflict
module supports **four data types**: binary, continuous, Poisson/count, and
survival.

---

# Supported Data Types

| Data type | `type =` argument | Conjugate update | Typical endpoint |
|---|---|---|---|
| Binary | `"binary"` | Beta–Binomial | Response rate, ORR |
| Continuous | `"continuous"` | Normal–Normal | Mean difference, HbA1c |
| Poisson / count | `"poisson"` | Gamma–Poisson | Adverse event rate |
| Survival | `"survival"` | Gamma–Exponential | Hazard rate, OS, PFS |

---

# Univariate Conflict: `prior_conflict()`

## Binary Data

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

```{r binary-none}
# No-conflict: data consistent with prior
cd_none <- prior_conflict(
  prior        = prior,
  data_summary = list(type = "binary", x = 13, n = 40),
  alpha        = 0.05
)
print(cd_none)
```

```{r binary-mild}
# Mild conflict
cd_mild <- prior_conflict(
  prior        = prior,
  data_summary = list(type = "binary", x = 20, n = 40)
)
print(cd_mild)
```

```{r binary-severe}
# Severe conflict
cd_severe <- prior_conflict(
  prior        = prior,
  data_summary = list(type = "binary", x = 35, n = 40)
)
print(cd_severe)
```

---

# The Diagnostics Explained

## 1. Box's Prior Predictive P-value

$$p_{\text{Box}} = 2\Phi\!\left(-\left|\frac{\hat\theta - \mu_\pi}
{\sqrt{\sigma_\pi^2 + \text{SE}^2}}\right|\right)$$

**Interpretation:** $p < 0.05$ flags conflict at the 5% level.

### Exact prior predictive version

The formula above is the Normal approximation to the prior and to the
likelihood. Where the prior predictive distribution has a closed form,
`prior_conflict()` uses it instead by default: Beta-Binomial (a Beta prior
with `type = "binary"` data) and Gamma-Poisson / Negative Binomial (a Gamma
prior with `type = "poisson"` or `"survival"` data). The approximation is
usually adequate for moderate samples and weakly informative priors, but
for an informative prior and a small sample it can flag conflict much more
often than its nominal level when there is none (for example when zero
events are observed, which the approximation treats as having zero sampling
variability). The exact version does not have this problem: it is the
prior predictive check of Box (1980) in the form used by Evans and Moshonov
(2006), the total prior predictive probability of outcomes no more probable
than the observed one. Use
`exact = FALSE` to obtain the approximation, for instance to reproduce
results from bayprior 0.4.0 or earlier.

```{r exact-binary}
cd_approx <- prior_conflict(
  prior        = prior,
  data_summary = list(type = "binary", x = 35, n = 40),
  exact        = FALSE
)
cd_exact <- prior_conflict(
  prior        = prior,
  data_summary = list(type = "binary", x = 35, n = 40)   # default: exact
)
data.frame(
  method     = c(cd_approx$box_pvalue_method, cd_exact$box_pvalue_method),
  box_pvalue = signif(c(cd_approx$box_pvalue, cd_exact$box_pvalue), 3),
  s_value    = round(c(cd_approx$s_value, cd_exact$s_value), 1)
)
```

For prior and data combinations without a closed-form check, the Normal
approximation is used (with `exact = TRUE`, with a message), and
`$box_pvalue_method` records which version was used (`"exact"` or
`"normal_approx"`). The other diagnostics below are not affected by `exact`.

## 2. S-value (Surprisal)

$$S = -\log_2(p_{\text{Box}}) \quad \text{(bits)}$$

`prior_conflict()` also reports `s_value`, a reparametrisation of the Box
p-value above onto a continuous, additive scale of evidence. It is
deliberately **not** a thresholded diagnostic: prior-data
conflict checks are typically not powered for a binary accept/reject
decision, so `s_value` is reported as an exploratory companion to
`box_pvalue` rather than replacing it -- see the note under
[Threshold Summary](#threshold-summary) below for why it has no severity
band of its own.

**Interpretation:** As a rough anchor, $S = 4.3$ bits corresponds to $p =
0.05$ (as surprising as 4-5 consecutive coin flips landing the same way),
and $S = 6.6$ bits to $p = 0.01$. Larger is more surprising; there is no
recommended cutoff.

Note `s_value` and `surprise_index` (next) are easily confused by name but
measure different things -- see `?prior_conflict` for the distinction.

## 3. Surprise Index

$$S = \left|\frac{\hat\theta - \mu_\pi}
{\sqrt{\sigma_\pi^2 + \text{SE}^2}}\right|$$

**Interpretation:** $S > 2$ moderate; $S > 3$ high surprise. When the Box
p-value is exact, `surprise_index` is the Normal-equivalent distance of that
p-value, $\Phi^{-1}(1 - p/2)$, so it agrees with `box_pvalue` and the
severity; `surprise_index_normal` keeps the distance above.

## 4. KL Divergence

$$\text{KL}(P \| Q) = \log\frac{\sigma_Q}{\sigma_P}
+ \frac{\sigma_P^2 + (\mu_P - \mu_Q)^2}{2\sigma_Q^2} - \frac12$$

**Interpretation:** $> 1$ indicates large information distance.

## 5. Bhattacharyya Overlap Coefficient

$$\text{BC} = \exp\!\left(
-\frac{(\mu_P - \mu_Q)^2}{4(\sigma_P^2 + \sigma_Q^2)}\right)
\cdot \left(\frac{2\sigma_P\sigma_Q}
{\sigma_P^2 + \sigma_Q^2}\right)^{1/2}$$

**Interpretation:** 1 = identical distributions; 0 = no overlap.

## Threshold Summary {#threshold-summary}

```{r thresholds, echo=FALSE}
knitr::kable(data.frame(
  Diagnostic        = c("Box p-value", "Surprise index",
                        "KL divergence", "Bhattacharyya overlap"),
  `No conflict`     = c(">= 0.05",   "< 2",    "< 0.5", "> 0.6"),
  `Mild conflict`   = c("0.0027-0.05", "2-3",   "0.5-1", "0.3-0.6"),
  `Severe conflict` = c("< 0.0027",  "> 3",    "> 1",   "< 0.3"),
  check.names = FALSE
), align = "lrrr")
```

Only the Box p-value determines `conflict_severity` (at the default
`alpha = 0.05`; "severe" is $p < 2\Phi(-3) \approx 0.0027$, the two-sided
tail area of a 3-standard-deviation surprise). The other three columns are
rough guides to reading the standardised distance, the divergence and the
overlap; they are conventions of this package, not calibrated cut-offs, and
we know of no published calibration for them. A p-value above `alpha` is the
absence of a flag, not evidence that the prior and the data agree: with a
small sample the check has little power to detect conflict.

`s_value` is intentionally absent from this table: it is reported as a
continuous, exploratory companion to `box_pvalue`, not as a fifth
thresholded diagnostic -- see [S-value (Surprisal)](#s-value-surprisal)
above.

---

# Prior-Likelihood-Posterior Overlay

```{r overlay-none}
plot_prior_likelihood(
  prior,
  data_summary   = list(type = "binary", x = 13, n = 40),
  show_posterior = TRUE
)
```

*No conflict: prior and likelihood overlap substantially.*

```{r overlay-severe}
plot_prior_likelihood(
  prior,
  data_summary   = list(type = "binary", x = 35, n = 40),
  show_posterior = TRUE
)
```

*Severe conflict: prior (blue) and likelihood (orange) widely separated.*

---

# Continuous Endpoint Conflict

```{r continuous}
prior_norm <- elicit_normal(
  mean = 0.0, sd = 0.3, method = "moments",
  label = "Log odds ratio"
)

# No conflict: observed log OR consistent with prior
prior_conflict(
  prior_norm,
  data_summary = list(type = "continuous", x = 0.15, sd = 0.20, n = 80)
)
```

---

# Poisson / Count Data Conflict {#poisson}

Poisson/count data arises for endpoints such as adverse event rates, where the
observed data consists of a count of events over an exposure period
(person-time). The conjugate update uses a **Gamma-Poisson** model:

$$\text{Prior: Gamma}(\alpha, \beta) \quad\Rightarrow\quad
\text{Posterior: Gamma}(\alpha + x,\; \beta + n)$$

where $x$ = observed events and $n$ = total person-time.

```{r poisson-setup}
# Gamma prior on adverse event rate
# Expert believes the AE rate is ~0.15 events per person-year (SD 0.06)
prior_ae <- elicit_gamma(
  mean      = 0.15,
  sd        = 0.06,
  method    = "moments",
  label     = "Adverse event rate (per person-year)",
  expert_id = "Safety Expert"
)
print(prior_ae)
```

```{r poisson-noconflict}
# Observed: 12 AEs over 100 person-years => rate 0.12 — consistent with prior
cd_pois_none <- prior_conflict(
  prior        = prior_ae,
  data_summary = list(type = "poisson", x = 12, n = 100),
  alpha        = 0.05
)
print(cd_pois_none)
```

```{r poisson-conflict}
# Observed: 40 AEs over 100 person-years => rate 0.40 — much higher than prior
cd_pois_conf <- prior_conflict(
  prior        = prior_ae,
  data_summary = list(type = "poisson", x = 40, n = 100)
)
print(cd_pois_conf)
```

```{r poisson-overlay}
plot_prior_likelihood(
  prior_ae,
  data_summary   = list(type = "poisson", x = 40, n = 100),
  show_posterior = TRUE
)
```

---

# Survival / Time-to-Event Data Conflict {#survival}

Survival data conflict arises for endpoints like OS and PFS hazard rates. The
observed data consists of event counts over a total follow-up time
(person-months). The conjugate update uses a **Gamma-Exponential** model:

$$\text{Prior: Gamma}(\alpha, \beta) \quad\Rightarrow\quad
\text{Posterior: Gamma}(\alpha + d,\; \beta + T)$$

where $d$ = observed events and $T$ = total follow-up time.

An Exponential prior (constant hazard) is equivalent to $\text{Gamma}(1, \lambda)$.

```{r survival-setup}
# Exponential prior on hazard rate
# Expert believes median OS ≈ 20 months => hazard ≈ 0.05
prior_hz <- elicit_exponential(
  mean      = 0.05,
  method    = "moments",
  label     = "OS hazard rate",
  expert_id = "Clinical Expert"
)
print(prior_hz)
```

```{r survival-noconflict}
# Observed: 20 deaths over 400 person-months => hazard 0.05 — consistent
cd_surv_none <- prior_conflict(
  prior        = prior_hz,
  data_summary = list(type = "survival", x = 20, n = 400),
  alpha        = 0.05
)
print(cd_surv_none)
```

```{r survival-conflict}
# Observed: 60 deaths over 400 person-months => hazard 0.15 — much higher
cd_surv_conf <- prior_conflict(
  prior        = prior_hz,
  data_summary = list(type = "survival", x = 60, n = 400)
)
print(cd_surv_conf)
```

```{r survival-overlay}
plot_prior_likelihood(
  prior_hz,
  data_summary   = list(type = "survival", x = 60, n = 400),
  show_posterior = TRUE
)
```

Also works with a Gamma prior on the hazard rate:

```{r survival-gamma}
prior_hz_gamma <- elicit_gamma(
  mean = 0.05, sd = 0.02, method = "moments",
  label = "OS hazard rate (Gamma prior)"
)
prior_conflict(
  prior_hz_gamma,
  data_summary = list(type = "survival", x = 20, n = 400)
)
```

---

# Compatibility Validation

bayprior refuses a prior family that cannot describe the data type, rather
than silently falling back to a Normal approximation of the prior's moments.
A Beta prior describes a probability, so it needs binary data; a Gamma,
lognormal, exponential or Weibull prior describes a positive rate or mean and
cannot be combined with binary data; a Normal prior works with every data
type:

```{r compat-error}
prior_beta <- elicit_beta(mean = 0.15, sd = 0.06, method = "moments",
                           label = "Rate")
res <- tryCatch(
  prior_conflict(prior_beta,
                 data_summary = list(type = "poisson", x = 15, n = 100)),
  error = function(e) conditionMessage(e))
res
```

---

# Mixture Prior Conflict

Mixture priors (robust priors, pooled consensus) are supported throughout.
The Box p-value uses the mixture's prior predictive distribution, the
weighted mixture of the components' predictives (exact for Beta or Gamma
components; computed by numerical integration if a component has no conjugate form,
with that component truncated to the parameter space). The other diagnostics use a
Normal approximation from the mixture's `fit_summary`:

```{r mixture-conflict}
e1 <- elicit_beta(mean = 0.25, sd = 0.08, method = "moments",
                  expert_id = "E1", label = "ORR")
e2 <- elicit_beta(mean = 0.40, sd = 0.12, method = "moments",
                  expert_id = "E2", label = "ORR")
mix <- aggregate_experts(list(E1 = e1, E2 = e2), weights = c(0.5, 0.5))

prior_conflict(mix, list(type = "binary", x = 18, n = 40))
```

---

# Multivariate Conflict: `conflict_mahalanobis()`

For trials with co-primary endpoints, the Mahalanobis distance provides an
omnibus multivariate conflict test:

```{r mahalanobis}
conflict_mahalanobis(
  prior_means = c(0.35, 0.60),
  prior_cov   = matrix(c(0.010, 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"),
  alpha       = 0.05
)
```

---

# Conflict Severity Classification

| Severity | Recommended action |
|---|---|
| **None** | Proceed — prior is consistent with data |
| **Mild** | Consider reporting both prior-weighted and likelihood-only estimates |
| **Severe** | Revise prior or use a robust/power prior; report sensitivity |

---

# References

- Evans, M. and Moshonov, H. (2006). Checking for prior-data conflict. *Bayesian Analysis*, 1(4).
- Box, G. E. P. (1980). Sampling and Bayes' inference in scientific modelling and robustness. *JRSS-A*, 143, 383–430.
- Bhattacharyya, A. (1943). On a measure of divergence between two statistical populations. *Bulletin of the Calcutta Mathematical Society*, 35, 99–109.
- Greenland, S. (2023). Divergence versus decision P-values: A distinction worth making in theory and keeping in practice. *Scandinavian Journal of Statistics*, 50(1), 54–88.