---
title: "Sensitivity Analysis"
author: "Ndoh Penn"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Sensitivity Analysis}
  %\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 Sensitivity Analysis Is Required

The FDA's 2026 draft guidance, *Use of Bayesian Methodology in Clinical
Trials of Drug and Biological Products*, calls for sensitivity analyses
that demonstrate trial conclusions are robust to plausible variations in
prior specification.

**bayprior** provides two dedicated functions:

- `sensitivity_grid()` — evaluates posterior mean, SD, and efficacy probability
- `sensitivity_cri()` — focuses on credible interval width and bounds

> **Sensitivity analysis is fully independent of conflict diagnostics.** You can
> run it without first running `prior_conflict()` by supplying `data_summary`
> directly.

---

# Supported Data Types

| `type =` | Conjugate update | `param_grid` names |
|---|---|---|
| `"binary"` | Beta-Binomial | `alpha`, `beta` |
| `"continuous"` | Normal-Normal | `mu`, `sigma` |
| `"poisson"` | Gamma-Poisson | `shape`, `rate` |
| `"survival"` | Gamma-Exponential | `shape`, `rate` |

---

# Binary Data

```{r binary-setup}
prior    <- elicit_beta(mean = 0.30, sd = 0.10, method = "moments",
                        label = "Response rate")
data_obs <- list(type = "binary", x = 14, n = 40)
```

```{r binary-grid}
sa <- sensitivity_grid(
  prior        = prior,
  data_summary = data_obs,
  param_grid   = list(alpha = seq(1, 8, 1), beta = seq(2, 20, 2)),
  target       = c("posterior_mean", "prob_efficacy"),
  threshold    = 0.30
)
sa$influence_scores
```

```{r binary-tornado}
plot_tornado(sa)
```

```{r binary-heatmap}
plot_sensitivity(sa, target = "posterior_mean")
```

## Credible Interval Sensitivity

```{r binary-cri}
cri_sa <- sensitivity_cri(
  prior        = prior,
  data_summary = data_obs,
  param_grid   = list(alpha = seq(1, 8, 1), beta = seq(2, 20, 2)),
  cri_level    = 0.95
)
cri_sa$influence_scores
plot_sensitivity(cri_sa, target = "cri_width")
```

---

# Poisson / Count Data {#poisson}

Sensitivity with Poisson data uses Gamma-Poisson conjugate updating.

```{r poisson-grid}
prior_ae  <- elicit_gamma(mean = 0.15, sd = 0.06, method = "moments",
                           label = "AE rate (per person-year)")
data_pois <- list(type = "poisson", x = 18, n = 120)

sa_pois <- sensitivity_grid(
  prior        = prior_ae,
  data_summary = data_pois,
  param_grid   = list(shape = seq(2, 10, 1), rate = seq(5, 40, 5)),
  target       = c("posterior_mean", "prob_efficacy"),
  threshold    = 0.20
)
sa_pois$influence_scores
plot_tornado(sa_pois)
```

```{r poisson-heatmap}
plot_sensitivity(sa_pois, target = "posterior_mean")
```

---

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

Sensitivity with survival data uses Gamma-Exponential conjugate updating.

```{r survival-grid}
prior_hz  <- elicit_exponential(mean = 0.05, method = "moments",
                                 label = "OS hazard rate")
data_surv <- list(type = "survival", x = 30, n = 600)

sa_surv <- sensitivity_grid(
  prior        = prior_hz,
  data_summary = data_surv,
  param_grid   = list(shape = seq(1, 5, 0.5), rate = seq(5, 30, 5)),
  target       = c("posterior_mean", "prob_efficacy"),
  threshold    = 0.10
)
sa_surv$influence_scores
plot_tornado(sa_surv)
```

```{r survival-cri}
cri_surv <- sensitivity_cri(
  prior        = prior_hz,
  data_summary = data_surv,
  param_grid   = list(shape = seq(1, 5, 0.5), rate = seq(5, 30, 5)),
  cri_level    = 0.95
)
plot_sensitivity(cri_surv, target = "cri_width")
```

---

# Continuous Endpoints

```{r continuous}
prior_cont <- elicit_normal(mean = 0.0, sd = 0.3, method = "moments",
                             label = "Log odds ratio")
sa_cont <- sensitivity_grid(
  prior        = prior_cont,
  data_summary = list(type = "continuous", x = 0.20, sd = 0.25, n = 60),
  param_grid   = list(mu = seq(-0.5, 0.5, 0.1), sigma = seq(0.1, 0.8, 0.1)),
  target       = c("posterior_mean", "posterior_sd")
)
plot_tornado(sa_cont)
```

---

# Interpreting Influence Scores

| Score | Sensitivity | Implication |
|---|---|---|
| < 0.05 | Not sensitive | Prior has negligible influence |
| 0.05 – 0.15 | Moderate | Report alongside primary estimate |
| > 0.15 | Sensitive | Consider robust prior; emphasise data |

The influence score is the range of a posterior quantity across the
hyperparameter grid, so it depends on the grid you chose and on the scale of
the quantity. The bands above are a rule of thumb for quantities on a 0 to 1
scale (a posterior mean of a proportion, a posterior probability); they are
not a calibrated or published cut-off, and they do not apply as they stand
to a log-odds, a rate or a hazard ratio. Report the grid alongside the scores,
and judge them against the size of a difference that would change a decision.
For comparison, the `priorsense` package flags prior sensitivity at 0.05 on a
different, normalised measure (a cumulative Jensen-Shannon distance); its
authors describe that cut-off as a guideline to adapt to the application
(Kallioinen et al., *Statistics and Computing* 34, article 57), and it is not comparable with the scores here.

---

# Choosing the grid

A grid over the raw hyperparameters (`alpha`, `beta`) mixes two different
questions: how *informative* the prior is, and where it is *centred*. Scaling
both by 0.5 to 2 for a Beta(6.7, 11.9) prior, for example, covers prior means
from 0.12 to 0.69, which is a much wider range of priors than is usually
plausible, and the influence scores then mostly reflect that width.

For Beta, Gamma and Normal priors the app therefore builds its default grid
in terms of the prior itself: prior strength (the effective sample size, for
a Beta prior `alpha + beta`) from 0.5 to 2 times the elicited value, and
optionally the prior mean shifted by up to a number of prior standard
deviations. In code, pass a data frame with one row per grid point; rows are
used as given, so hyperparameters can move together:

```{r strength-grid}
prior <- elicit_beta(mean = 0.36, sd = 0.10, method = "moments")
S <- prior$params$alpha + prior$params$beta   # effective sample size
m <- prior$fit_summary$mean
sd_prior <- prior$fit_summary$sd

grid <- expand.grid(mean = m + c(-1, 0, 1) * sd_prior,
                    strength = S * c(0.5, 1, 2))
param_grid <- data.frame(alpha = grid$mean * grid$strength,
                         beta  = (1 - grid$mean) * grid$strength,
                         prior_mean = grid$mean,
                         prior_strength = grid$strength)

sa_design <- sensitivity_grid(
  prior,
  data_summary = list(type = "binary", x = 18, n = 40),
  param_grid   = param_grid
)
sa_design$influence_scores
```

State what the grid covers when you report influence scores. Neither choice
is "correct": a wider or narrower range of plausible priors is a judgement
about the application, and the scores should be read against it.

---

# Mixture Prior Sensitivity

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

sa_mix <- sensitivity_grid(
  prior        = mix,
  data_summary = list(type = "binary", x = 14, n = 40),
  param_grid   = list(alpha = seq(1, 8, 1), beta = seq(2, 16, 2)),
  target       = "posterior_mean"
)
plot_tornado(sa_mix)
```

Note: for mixture priors, the grid varies the dominant component's parameters.
A compatibility warning is shown in the Shiny app.

---

# Regulatory Reporting

Include in the clinical study report:

1. Influence score table — one row per posterior quantity
2. Tornado plot — ordered from most to least sensitive
3. Heatmap — for the primary efficacy quantity
4. Classification statement — sensitive / not sensitive per quantity

All three are generated automatically by `prior_report()` when a
`bayprior_sensitivity` object is supplied.

---

# References

- FDA (2026). *Use of Bayesian Methodology in Clinical Trials of Drug and Biological Products* (Draft Guidance for Industry).
- ICH E9(R1) (2019). *Addendum on Estimands and Sensitivity Analysis*.