R/conflict_sensitivity.R

In bayprior: Bayesian Prior Elicitation and Diagnostics for Clinical Trials

#' Compute prior-data conflict diagnostics
#'
#' Evaluates conflict between a specified prior and observed data using
#' multiple complementary diagnostics: Box's (1980) predictive p-value,
#' the surprise index (standardised distance), Kullback-Leibler divergence,
#' and the Bhattacharyya overlap coefficient between the prior and the
#' (normalised) likelihood.
#'
#' @param prior       A \code{bayprior} object.
#' @param data_summary Named list describing the observed data:
#'   \describe{
#'     \item{\code{type}}{\code{"binary"}, \code{"continuous"},
#'       \code{"poisson"}, or \code{"survival"}.}
#'     \item{\code{x}}{Number of events (binary / poisson / survival)
#'       or observed mean (continuous).}
#'     \item{\code{n}}{Sample size (binary / continuous), total exposure
#'       (poisson: person-time), or total follow-up time (survival).}
#'     \item{\code{sd}}{Observed standard deviation (continuous only).}
#'   }
#' @param alpha Numeric. Significance level for the Box p-value flag.
#'   Default \code{0.05}.
#'
#' @return An object of class \code{bayprior_conflict} containing:
#'   \describe{
#'     \item{\code{box_pvalue}}{Box's prior predictive p-value.}
#'     \item{\code{surprise_index}}{Standardised distance between prior mean
#'       and observed data.}
#'     \item{\code{kl_prior_likelihood}}{KL divergence from prior to likelihood.}
#'     \item{\code{overlap}}{Bhattacharyya overlap coefficient in \[0, 1\].}
#'     \item{\code{conflict_severity}}{One of \code{"none"}, \code{"mild"},
#'       \code{"severe"}.}
#'     \item{\code{conflict_flag}}{Logical; \code{TRUE} if
#'       \code{box_pvalue < alpha}.}
#'     \item{\code{recommendation}}{Plain-language guidance string.}
#'     \item{\code{data_summary}}{The data summary passed in.}
#'     \item{\code{prior}}{The input prior.}
#'   }
#'
#' @references
#' Box, G. E. P. (1980). Sampling and Bayes' inference in scientific modelling
#' and robustness. \emph{Journal of the Royal Statistical Society A}, 143,
#' 383-430.
#'
#' @examples
#' prior <- elicit_beta(mean = 0.30, sd = 0.10, method = "moments",
#'                      label = "Response rate")
#' cd <- prior_conflict(prior, list(type = "binary", x = 18, n = 40))
#' print(cd)
#'
#' @importFrom rlang %||% abort
#' @export
prior_conflict <- function(prior, data_summary, alpha = 0.05) {

  if (!inherits(prior, "bayprior")) {
    rlang::abort("`prior` must be a bayprior object.")
  }

  type <- data_summary$type %||% "binary"
  n    <- data_summary$n
  x    <- data_summary$x

  # Approximate prior as Normal for analytic diagnostics
  prior_mean <- prior$fit_summary$mean
  prior_sd   <- prior$fit_summary$sd

  # Likelihood parameters -- approximate as Normal for analytic diagnostics
  if (type == "binary") {
    obs_mean <- x / n
    obs_se   <- sqrt(obs_mean * (1 - obs_mean) / n)
  } else if (type %in% c("poisson", "survival")) {
    # Poisson: x events over exposure n; rate = x/n, SE via delta method
    # Survival: x events over total follow-up n; hazard = x/n
    obs_mean <- x / n
    obs_se   <- sqrt(x) / n        # delta method: SE(x/n) = sqrt(x)/n
  } else {
    obs_mean <- x
    obs_se   <- data_summary$sd / sqrt(n)
  }

  # Guard against obs_se = 0 (e.g. all successes / all failures),
  # which would cause division-by-zero in pred_sd, z, kl, and overlap.
  obs_se <- max(obs_se, 1e-8)

  # Box's prior predictive p-value
  pred_sd <- sqrt(prior_sd^2 + obs_se^2)
  z       <- (obs_mean - prior_mean) / pred_sd
  box_p   <- 2 * stats::pnorm(-abs(z))

  # Surprise index
  surprise <- abs(z)

  # KL divergence (normal approximation); see .kl_normal() for convention
  kl <- .kl_normal(prior_mean, prior_sd, obs_mean, obs_se)

  # Bhattacharyya overlap
  overlap <- .bhattacharyya_normal(prior_mean, prior_sd, obs_mean, obs_se)

  # Severity classification
  severity <- dplyr::case_when(
    box_p >= alpha               ~ "none",
    box_p < alpha & surprise < 3 ~ "mild",
    TRUE                         ~ "severe"
  )

  recommendation <- switch(severity,
    none = glue::glue(
      "No evidence of prior-data conflict (Box p = {round(box_p, 3)}). ",
      "The prior appears consistent with the observed data."
    ),
    mild = glue::glue(
      "Mild prior-data conflict detected (Box p = {round(box_p, 3)}, ",
      "surprise = {round(surprise, 2)}). ",
      "Consider reporting a sensitivity analysis with a more diffuse prior."
    ),
    severe = glue::glue(
      "Severe prior-data conflict detected (Box p = {round(box_p, 4)}, ",
      "surprise = {round(surprise, 2)}). ",
      "Re-elicitation or use of a robust/sceptical prior is strongly recommended."
    )
  )

  structure(
    list(
      box_pvalue          = box_p,
      surprise_index      = surprise,
      kl_prior_likelihood = kl,
      overlap             = overlap,
      conflict_severity   = severity,
      conflict_flag       = box_p < alpha,
      recommendation      = as.character(recommendation),
      data_summary        = data_summary,
      prior               = prior,
      prior_mean          = prior_mean,
      prior_sd            = prior_sd,
      obs_mean            = obs_mean,
      obs_se              = obs_se,
      alpha               = alpha
    ),
    class = "bayprior_conflict"
  )
}


#' Derive a single-family working prior from a mixture (internal)
#'
#' Sensitivity grids are defined over a single distribution family's
#' hyperparameters, so mixture priors (from \code{\link{aggregate_experts}}
#' or \code{\link{robust_prior}}) need a single-family stand-in. This
#' moment-matches the mixture's actual pooled mean and SD (from
#' \code{prior$fit_summary}) to the dominant component's family, so the
#' working prior reflects the full pooled information rather than
#' discarding all but one component. Falls back to the dominant component
#' alone, with a warning, if the family cannot be moment-matched from
#' mean/SD (e.g. \code{"exponential"}, \code{"weibull"}) or if
#' \code{fit_summary$sd} is unavailable (e.g. logarithmic pooling).
#'
#' Shared by \code{\link{sensitivity_grid}} and the Shiny sensitivity
#' module (\code{mod_sensitivity.R}) so that the parameter ranges shown in
#' the UI are always centred on the same working prior actually analysed.
#'
#' See \code{\link{aggregate_experts}} for the pooled mean/variance
#' formula, and \code{\link{elicit_beta}}, \code{\link{elicit_normal}},
#' \code{\link{elicit_gamma}}, \code{\link{elicit_lognormal}} for each
#' family's moment-matching identities.
#'
#' @param prior A \code{bayprior} object, mixture or single-family.
#' @param quiet Logical. If \code{TRUE}, suppress the informational message
#'   emitted when moment-matching succeeds (warnings for failures/fallback
#'   are still shown). Useful for reactive UI contexts. Default \code{FALSE}.
#'
#' @return A single-family \code{bayprior} object.
#' @keywords internal
.mixture_working_prior <- function(prior, quiet = FALSE) {

  if (!prior$dist %in% c("mixture", "log_pool")) return(prior)

  dominant    <- which.max(prior$weights)
  fam         <- prior$components[[dominant]]$dist
  fams        <- unique(vapply(prior$components, function(x) x$dist, character(1)))
  pooled_mean <- prior$fit_summary$mean
  pooled_sd   <- prior$fit_summary$sd

  if (length(fams) > 1) {
    rlang::warn(paste0(
      "[bayprior] mixture components have different distribution ",
      "families (", paste(fams, collapse = ", "), "). Moment-matching ",
      "the pooled mean/SD against the dominant component's family ('",
      fam, "')."
    ))
  }

  moment_matcher <- switch(fam,
    beta      = elicit_beta,
    normal    = elicit_normal,
    gamma     = elicit_gamma,
    lognormal = elicit_lognormal,
    NULL
  )

  matched <- if (is.null(moment_matcher) || is.null(pooled_sd) || is.na(pooled_sd)) {
    NULL
  } else {
    tryCatch(
      moment_matcher(mean = pooled_mean, sd = pooled_sd, method = "moments",
                     label = prior$label,
                     expert_id = "Pooled (moment-matched)"),
      error = function(e) NULL
    )
  }

  if (is.null(matched)) {
    rlang::warn(paste0(
      "[bayprior] could not moment-match the pooled '", fam, "' mixture ",
      "(mean = ", round(pooled_mean, 4), ", SD = ",
      round(pooled_sd %||% NA_real_, 4), "). Falling back to the dominant ",
      "component only (weight = ", round(prior$weights[dominant], 3),
      ", expert_id = '", prior$components[[dominant]]$expert_id, "'). ",
      "Results reflect this single component, not the full pooled prior."
    ))
    return(prior$components[[dominant]])
  }

  if (!quiet) {
    message(paste0(
      "[bayprior] using a moment-matched working prior (", toupper(fam),
      ", mean = ", round(pooled_mean, 4), ", SD = ", round(pooled_sd, 4),
      ") that represents the full pooled mixture, not a single component."
    ))
  }
  matched
}


#' Sensitivity grid over prior hyperparameters
#'
#' Evaluates how posterior inferences change as prior hyperparameters vary
#' over a specified grid. This is the core function for demonstrating
#' robustness of trial conclusions to prior choice.
#'
#' @param prior       A \code{bayprior} object (the reference prior).
#' @param data_summary Named list as for \code{\link{prior_conflict}}.
#' @param param_grid  Named list of numeric vectors, one per hyperparameter
#'   to vary. Names must match hyperparameter names in \code{prior$params}.
#'   Example: \code{list(alpha = seq(1, 8, 0.5), beta = seq(2, 20, 1))}.
#' @param target Character vector. Which posterior quantities to compute.
#'   Any of \code{"posterior_mean"}, \code{"posterior_sd"},
#'   \code{"prob_efficacy"}.
#' @param threshold Numeric. Efficacy threshold used in
#'   \code{Pr(theta > threshold)}. Default \code{0.30}.
#'
#' @details
#' The hyperparameter grid is defined over a single distribution family, so
#' when \code{prior} is a mixture (e.g. from \code{\link{aggregate_experts}}
#' or \code{\link{robust_prior}}), a single-family \emph{working prior} is
#' derived first, via the internal \code{.mixture_working_prior()} helper.
#' As of this version, that working prior is obtained in two steps: first,
#' the mixture's pooled mean and SD are computed exactly (see the mean/var
#' formula in \code{\link{aggregate_experts}}); second, those pooled moments
#' are matched to the dominant component's distribution family using that
#' family's own moment-matching identities (see \code{\link{elicit_beta}},
#' \code{\link{elicit_normal}}, \code{\link{elicit_gamma}}, or
#' \code{\link{elicit_lognormal}} for the specific formula used). The grid
#' therefore reflects the full pooled information, not just one component.
#' A message reports the working prior used. If the dominant component's
#' family cannot be moment-matched from mean/SD alone (\code{"exponential"}
#' or \code{"weibull"}), the function falls back to the dominant component
#' by weight and issues an explicit warning identifying which component
#' was used and why.
#'
#' @return An object of class \code{bayprior_sensitivity}.
#'
#' @examples
#' prior <- elicit_beta(mean = 0.30, sd = 0.10, method = "moments",
#'                      label = "Response rate")
#' sa <- sensitivity_grid(
#'   prior,
#'   data_summary = list(type = "binary", x = 14, n = 40),
#'   param_grid   = list(alpha = seq(1, 8, 0.5), beta = seq(2, 20, 1))
#' )
#' plot_tornado(sa)
#' plot_sensitivity(sa, target = "posterior_mean")
#'
#' @importFrom rlang %||% abort warn
#' @export
sensitivity_grid <- function(prior,
                              data_summary,
                              param_grid,
                              target    = c("posterior_mean", "posterior_sd",
                                            "prob_efficacy"),
                              threshold = 0.30) {

  target <- match.arg(target, several.ok = TRUE)

  type <- data_summary$type %||% "binary"
  n    <- data_summary$n
  x    <- data_summary$x

  # -- For mixture priors, derive a single-family working prior ----------------
  # See .mixture_working_prior() for the moment-matching logic; shared with
  # the Shiny sensitivity module so UI defaults stay consistent with what
  # is actually analysed here.
  working_prior <- .mixture_working_prior(prior)

  # -- Auto-remap param_grid names to prior hyperparameter names ----------------
  #
  # The Shiny UI generates generic names ("param1", "param2", ?) that will
  # never literally match prior hyperparameter names ("alpha"/"beta", etc.).
  # We resolve this with a three-step strategy:
  #
  #   1. Exact match  -- use as-is (names already correct).
  #   2. Positional remap -- if no names match but the *count* equals the number
  #      of prior hyperparameters, rename param_grid entries to the prior's
  #      param names in order and warn the caller.
  #   3. Abort -- counts also differ; the mapping is genuinely ambiguous.
  #
  prior_param_names <- names(working_prior$params)
  grid_param_names  <- names(param_grid)
  valid_names       <- intersect(grid_param_names, prior_param_names)

  if (length(valid_names) == 0) {
    # No exact matches -- attempt positional remap
    if (length(grid_param_names) == length(prior_param_names)) {
      # Downgraded to a message (not a warning) because this is expected
      # behaviour when the Shiny UI passes generic names like param1/param2.
      # rlang::warn() causes a Shiny warning banner; message() does not.
      message(paste0(
        "[bayprior] sensitivity_grid: remapping param_grid names positionally: ",
        paste(grid_param_names, "->", prior_param_names, collapse = ", "), "."
      ))
      names(param_grid) <- prior_param_names
      valid_names       <- prior_param_names
    } else {
      rlang::abort(paste0(
        "Cannot map param_grid to prior hyperparameters: names don't match ",
        "and counts differ.\n",
        "  param_grid names (", length(grid_param_names), "): ",
        paste(grid_param_names,  collapse = ", "), "\n",
        "  Prior param names (", length(prior_param_names), "): ",
        paste(prior_param_names, collapse = ", "), "\n",
        "Either rename param_grid entries to match the prior's hyperparameter ",
        "names, or supply exactly ", length(prior_param_names),
        " grid vector(s) in the same order."
      ))
    }
  }

  # Build full Cartesian grid (after any name remapping above)
  grid_df <- do.call(expand.grid, param_grid)

  # Identify reference row closest to working prior's parameters.
  # mapply() returns a plain vector when param_grid has only one entry,
  # which would make rowSums() produce wrong scalar results. matrix() with
  # explicit nrow handles both single- and multi-parameter cases correctly.
  dist_matrix <- matrix(
    mapply(
      function(col, ref) (grid_df[[col]] - ref)^2,
      valid_names,
      working_prior$params[valid_names]
    ),
    nrow = nrow(grid_df)
  )
  dists   <- rowSums(dist_matrix)
  ref_row <- which.min(dists)

  # Evaluate posterior summaries at each grid point
  results <- purrr::map_dfr(seq_len(nrow(grid_df)), function(i) {

    row    <- grid_df[i, , drop = FALSE]
    params <- as.list(row)

    # Build temporary prior at this grid point
    tmp_prior <- tryCatch(
      .make_bayprior(working_prior$dist, params, working_prior$method,
                     working_prior$expert_id, working_prior$label,
                     working_prior$input),
      error = function(e) NULL
    )
    if (is.null(tmp_prior)) {
      out <- as.list(row)
      for (t in target) out[[t]] <- NA_real_
      return(as.data.frame(out))
    }

    # Compute posterior via conjugate update
    post <- tryCatch(
      .conjugate_update(tmp_prior, data_summary),
      error = function(e) NULL
    )
    if (is.null(post)) {
      out <- as.list(row)
      for (t in target) out[[t]] <- NA_real_
      return(as.data.frame(out))
    }

    post_s <- post$fit_summary
    out    <- as.list(row)

    if ("posterior_mean" %in% target) out$posterior_mean <- post_s$mean
    if ("posterior_sd"   %in% target) out$posterior_sd   <- post_s$sd
    if ("prob_efficacy"  %in% target) {
      out$prob_efficacy <- tryCatch({
        if (post$dist == "beta") {
          stats::pbeta(threshold, post$params$alpha, post$params$beta,
                       lower.tail = FALSE)
        } else if (post$dist == "normal") {
          stats::pnorm(threshold, post$params$mu, post$params$sigma,
                       lower.tail = FALSE)
        } else if (post$dist == "gamma") {
          stats::pgamma(threshold, post$params$shape, post$params$rate,
                        lower.tail = FALSE)
        } else {
          # Mixture or unknown: normal approximation from fit_summary
          stats::pnorm(threshold, post_s$mean, post_s$sd,
                       lower.tail = FALSE)
        }
      }, error = function(e) NA_real_)
    }

    as.data.frame(out)
  })

  # Detect the case where every grid point produced NA (e.g. because
  # .make_bayprior() or .conjugate_update() failed for all rows) and abort
  # with a clear message before calling range() or building a colorscale --
  # otherwise this cascades into "no non-missing arguments to min/max"
  # warnings and a broken Plotly colorscale downstream.
  all_na_targets <- Filter(function(t) {
    v <- results[[t]]
    is.null(v) || all(is.na(v))
  }, target)

  if (length(all_na_targets) > 0) {
    rlang::abort(paste0(
      "All grid evaluations returned NA for: ",
      paste(all_na_targets, collapse = ", "), ".\n",
      "This usually means .make_bayprior() or .conjugate_update() failed at ",
      "every grid point. Check that:\n",
      "  1. param_grid values produce valid hyperparameters (e.g. alpha > 0).\n",
      "  2. The prior distribution supports conjugate updating with this data type.\n",
      "  3. data_summary$type is set correctly ('binary', 'continuous', 'poisson', or 'survival')."
    ))
  }

  # Warn about target columns missing from results entirely (distinct from
  # a target column that exists but is all-NA, handled above).
  missing_targets <- setdiff(target, names(results))
  if (length(missing_targets) > 0) {
    rlang::warn(paste0(
      "The following targets are missing from results entirely and will have ",
      "influence score 0: ", paste(missing_targets, collapse = ", ")
    ))
  }

  # Influence scores: range of each target across the grid.
  # Returns 0 when all values are NA or non-finite (after the all-NA guard above,
  # this only triggers for partial-NA columns, which is legitimate).
  influence <- vapply(target, function(t) {
    v <- results[[t]]
    if (is.null(v) || all(is.na(v))) return(0)
    fin <- v[is.finite(v)]
    if (length(fin) == 0) return(0)
    diff(range(fin))
  }, numeric(1))

  structure(
    list(
      grid             = results,
      param_grid       = param_grid,
      target           = target,
      reference_row    = ref_row,
      influence_scores = influence,
      threshold        = threshold,
      prior            = prior
    ),
    class = "bayprior_sensitivity"
  )
}



#' Credible interval sensitivity over prior hyperparameters
#'
#' Evaluates how the posterior credible interval width and bounds change as
#' prior hyperparameters vary over a specified grid. This is the preferred
#' function for demonstrating that key regulatory conclusions (e.g., whether
#' the CrI excludes a null value) are robust to prior choice.
#'
#' @param prior        A \code{bayprior} object (the reference prior).
#' @param data_summary Named list as for \code{\link{prior_conflict}}.
#' @param param_grid   Named list of numeric vectors, one per hyperparameter.
#' @param cri_level    Numeric in (0, 1). Credible interval level. Default
#'   \code{0.95}.
#' @param threshold    Optional numeric. Computes \code{Pr(theta > threshold)}
#'   at each grid point if supplied.
#'
#' @details
#' As with \code{\link{sensitivity_grid}}, when \code{prior} is a mixture,
#' a single-family working prior is derived via \code{.mixture_working_prior()}
#' -- moment-matching the mixture's pooled mean/SD to the dominant
#' component's family -- rather than analysing the dominant component alone.
#'
#' @return An object of class \code{bayprior_sensitivity} whose grid contains
#'   columns \code{cri_lower}, \code{cri_upper}, and \code{cri_width}, plus
#'   optionally \code{posterior_mean}, \code{posterior_sd}, and
#'   \code{prob_efficacy}.
#'
#' @examples
#' prior <- elicit_beta(mean = 0.30, sd = 0.10, method = "moments",
#'                      label = "Response rate")
#' cri_sa <- sensitivity_cri(
#'   prior,
#'   data_summary = list(type = "binary", x = 14, n = 40),
#'   param_grid   = list(alpha = seq(1, 8, 0.5), beta = seq(2, 20, 1)),
#'   cri_level    = 0.95
#' )
#' plot_sensitivity(cri_sa, target = "cri_width")
#'
#' @importFrom rlang %||% abort
#' @export
sensitivity_cri <- function(prior,
                             data_summary,
                             param_grid,
                             cri_level = 0.95,
                             threshold = NULL) {

  if (!inherits(prior, "bayprior"))
    rlang::abort("`prior` must be a bayprior object.")
  if (!is.numeric(cri_level) || cri_level <= 0 || cri_level >= 1)
    rlang::abort("`cri_level` must be a number strictly between 0 and 1.")

  alpha_lo <- (1 - cri_level) / 2
  alpha_hi <- 1 - alpha_lo

  # Moment-match a mixture prior's pooled mean/SD to a working prior in the
  # dominant component's family (see .mixture_working_prior() for details
  # and the same fix applied to sensitivity_grid()), rather than discarding
  # all but the dominant component.
  working_prior <- .mixture_working_prior(prior)

  # Remap generic param names (e.g. "param1", "param2") to prior param names
  prior_param_names <- names(working_prior$params)
  grid_param_names  <- names(param_grid)
  valid_names       <- intersect(grid_param_names, prior_param_names)

  if (length(valid_names) == 0) {
    if (length(grid_param_names) == length(prior_param_names)) {
      message(paste0(
        "[bayprior] sensitivity_cri: remapping param_grid names positionally: ",
        paste(grid_param_names, "->", prior_param_names, collapse = ", "), "."
      ))
      names(param_grid) <- prior_param_names
      valid_names       <- prior_param_names
    } else {
      rlang::abort(paste0(
        "Cannot map param_grid to prior hyperparameters. ",
        "Names don't match and counts differ.\n",
        "  param_grid: ", paste(grid_param_names, collapse = ", "), "\n",
        "  Prior params: ", paste(prior_param_names, collapse = ", ")
      ))
    }
  }

  grid_df <- do.call(expand.grid, param_grid)

  # Reference row: closest grid point to working prior's current parameters
  ref_diffs <- vapply(valid_names, function(nm) {
    (grid_df[[nm]] - (working_prior$params[[nm]] %||% 0))^2
  }, numeric(nrow(grid_df)))

  # Always treat as matrix to avoid apply() dimension errors on single-param grids
  ref_diffs <- matrix(ref_diffs, nrow = nrow(grid_df))
  ref_row   <- which.min(rowSums(ref_diffs))

  # Evaluate CrI at each grid point using purrr::map_dfr (never apply())
  target_cols <- c("cri_lower", "cri_upper", "cri_width", "posterior_mean",
                   "posterior_sd")
  if (!is.null(threshold)) target_cols <- c(target_cols, "prob_efficacy")

  results <- purrr::map_dfr(seq_len(nrow(grid_df)), function(i) {
    row    <- grid_df[i, , drop = FALSE]
    params <- as.list(row)

    tmp_prior <- tryCatch(
      .make_bayprior(working_prior$dist, params, working_prior$method,
                     working_prior$expert_id, working_prior$label,
                     working_prior$input),
      error = function(e) NULL
    )

    na_row <- as.data.frame(c(as.list(row),
                              list(cri_lower     = NA_real_,
                                   cri_upper     = NA_real_,
                                   cri_width     = NA_real_,
                                   posterior_mean = NA_real_,
                                   posterior_sd   = NA_real_)))
    if (!is.null(threshold)) na_row$prob_efficacy <- NA_real_
    if (is.null(tmp_prior)) return(na_row)

    post <- tryCatch(
      .conjugate_update(tmp_prior, data_summary),
      error = function(e) NULL
    )
    if (is.null(post)) return(na_row)

    s <- post$fit_summary
    out <- as.data.frame(as.list(row))

    # Compute posterior quantiles for CrI bounds
    bounds <- tryCatch({
      if (post$dist == "beta") {
        c(stats::qbeta(alpha_lo, post$params$alpha, post$params$beta),
          stats::qbeta(alpha_hi, post$params$alpha, post$params$beta))
      } else if (post$dist == "normal") {
        c(stats::qnorm(alpha_lo, post$params$mu, post$params$sigma),
          stats::qnorm(alpha_hi, post$params$mu, post$params$sigma))
      } else if (post$dist == "gamma") {
        c(stats::qgamma(alpha_lo, post$params$shape, post$params$rate),
          stats::qgamma(alpha_hi, post$params$shape, post$params$rate))
      } else {
        # Normal approximation for mixture / other
        c(stats::qnorm(alpha_lo, s$mean, s$sd),
          stats::qnorm(alpha_hi, s$mean, s$sd))
      }
    }, error = function(e) c(NA_real_, NA_real_))

    out$cri_lower      <- bounds[1]
    out$cri_upper      <- bounds[2]
    out$cri_width      <- bounds[2] - bounds[1]
    out$posterior_mean <- s$mean
    out$posterior_sd   <- s$sd

    if (!is.null(threshold)) {
      out$prob_efficacy <- tryCatch({
        if (post$dist == "beta")
          stats::pbeta(threshold, post$params$alpha, post$params$beta,
                       lower.tail = FALSE)
        else if (post$dist == "normal")
          stats::pnorm(threshold, post$params$mu, post$params$sigma,
                       lower.tail = FALSE)
        else if (post$dist == "gamma")
          stats::pgamma(threshold, post$params$shape, post$params$rate,
                        lower.tail = FALSE)
        else
          stats::pnorm(threshold, s$mean, s$sd, lower.tail = FALSE)
      }, error = function(e) NA_real_)
    }

    out
  })

  report_targets <- intersect(target_cols, names(results))

  influence <- vapply(report_targets, function(t) {
    v <- results[[t]]
    if (is.null(v) || all(is.na(v))) return(0)
    fin <- v[is.finite(v)]
    if (length(fin) == 0) return(0)
    diff(range(fin))
  }, numeric(1))

  structure(
    list(
      grid             = results,
      param_grid       = param_grid,
      target           = report_targets,
      reference_row    = ref_row,
      influence_scores = influence,
      cri_level        = cri_level,
      threshold        = threshold,
      prior            = prior
    ),
    class = "bayprior_sensitivity"
  )
}



# Define %||% explicitly so it is available regardless of whether rlang is
# attached (it may only be in Imports, not Depends). Functions in this file
# that carry @importFrom rlang %||% will also satisfy R CMD CHECK.
`%||%` <- rlang::`%||%`


# KL(P||Q): KL divergence from P = Normal(m1, s1) to Q = Normal(m2, s2).
# Formula: log(s2/s1) + (s1^2 + (m1-m2)^2) / (2*s2^2) - 1/2
# Argument order is (from_mean, from_sd, to_mean, to_sd) -- kept explicit
# here to prevent future sign/argument confusion.
.kl_normal <- function(m1, s1, m2, s2) {
  log(s2 / s1) + (s1^2 + (m1 - m2)^2) / (2 * s2^2) - 0.5
}

.bhattacharyya_normal <- function(m1, s1, m2, s2) {
  t1 <- 0.25 * log(0.25 * (s1^2 / s2^2 + s2^2 / s1^2 + 2))
  t2 <- 0.25 * (m1 - m2)^2 / (s1^2 + s2^2)
  exp(-(t1 + t2))
}

.conjugate_update <- function(prior, data_summary) {
  type <- data_summary$type %||% "binary"
  n    <- data_summary$n
  x    <- data_summary$x

  # -- Beta / binary ----------------------------------------------------------
  if (prior$dist == "beta" && type == "binary") {
    a_post <- prior$params$alpha + x
    b_post <- prior$params$beta  + (n - x)
    return(.make_bayprior("beta", list(alpha = a_post, beta = b_post),
                          "posterior", prior$expert_id, prior$label, list()))
  }

  # -- Normal -----------------------------------------------------------------
  if (prior$dist == "normal") {
    obs_mean  <- x
    obs_se    <- (data_summary$sd %||% prior$fit_summary$sd) / sqrt(n)
    prior_var <- prior$params$sigma^2
    lik_var   <- obs_se^2
    post_var  <- 1 / (1 / prior_var + 1 / lik_var)
    post_mean <- post_var * (prior$params$mu / prior_var + obs_mean / lik_var)
    return(.make_bayprior("normal", list(mu = post_mean, sigma = sqrt(post_var)),
                          "posterior", prior$expert_id, prior$label, list()))
  }

  # -- Exponential / Poisson or survival ---------------------------------------
  # Exponential prior on rate is equivalent to Gamma(1, rate) parameterisation.
  # Conjugate update with Poisson/survival data: posterior is still Gamma.
  # Map to Gamma update: shape = 1 + x, rate = rate + n.
  if (prior$dist == "exponential" && type %in% c("poisson", "survival", "binary")) {
    # Represent as Gamma(1, lambda) for conjugate update
    lambda     <- prior$params$rate
    shape_post <- 1 + x
    rate_post  <- lambda + n
    return(.make_bayprior("gamma", list(shape = shape_post, rate = rate_post),
                          "posterior", prior$expert_id, prior$label, list()))
  }

  # -- Weibull -- Normal approximation ------------------------------------------
  # No closed-form conjugate update for Weibull. Approximate via a Normal
  # distribution matched to the posterior mean and SD using the prior's
  # fit_summary as the prior parameters.
  if (prior$dist == "weibull") {
    prior_mean <- prior$fit_summary$mean
    prior_sd   <- prior$fit_summary$sd
    if (type %in% c("poisson", "survival")) {
      obs_mean_w <- x / n
      obs_se_w   <- sqrt(x) / n
    } else if (type == "binary") {
      obs_mean_w <- x / n
      obs_se_w   <- sqrt(obs_mean_w * (1 - obs_mean_w) / n)
    } else {
      obs_mean_w <- x
      obs_se_w   <- (data_summary$sd %||% prior_sd) / sqrt(n)
    }
    obs_se_w   <- max(obs_se_w, 1e-8)
    prior_var  <- prior_sd^2
    lik_var    <- obs_se_w^2
    post_var   <- 1 / (1 / prior_var + 1 / lik_var)
    post_mean  <- post_var * (prior_mean / prior_var + obs_mean_w / lik_var)
    return(.make_bayprior("normal",
                          list(mu = post_mean, sigma = sqrt(post_var)),
                          "posterior", prior$expert_id, prior$label, list()))
  }

  # -- Gamma / Poisson or survival ---------------------------------------------
  # Gamma(shape, rate) prior on rate lambda; Poisson(lambda * n) likelihood.
  # Posterior: Gamma(shape + x, rate + n).
  # For survival data: Gamma prior on hazard; Exponential(lambda) lifetimes.
  # Posterior: Gamma(shape + events, rate + total_follow_up).
  if (prior$dist == "gamma" && type %in% c("poisson", "survival")) {
    shape_post <- prior$params$shape + x
    rate_post  <- prior$params$rate  + n
    return(.make_bayprior("gamma", list(shape = shape_post, rate = rate_post),
                          "posterior", prior$expert_id, prior$label, list()))
  }

  # -- Gamma ------------------------------------------------------------------
  # For continuous data, x is the observed mean (not a total count), so the
  # total must be recovered before updating the Gamma rate. Prefer an
  # explicit x_sum field when available (e.g. Poisson count data) and fall
  # back to n * x otherwise.
  if (prior$dist == "gamma" && type == "continuous") {
    x_sum      <- data_summary$x_sum %||% (x * n)
    shape_post <- prior$params$shape + x_sum
    rate_post  <- prior$params$rate  + n
    return(.make_bayprior("gamma", list(shape = shape_post, rate = rate_post),
                          "posterior", prior$expert_id, prior$label, list()))
  }

  # -- Mixture -- update each component and re-weight by marginal likelihood ---
  if (prior$dist == "mixture") {
    components <- prior$components
    weights    <- prior$weights

    # Update each component individually
    post_components <- lapply(components, function(comp) {
      tryCatch(.conjugate_update(comp, data_summary), error = function(e) NULL)
    })

    # Drop components that failed to update
    keep <- !vapply(post_components, is.null, logical(1))
    if (!any(keep)) {
      rlang::abort("Could not update any mixture component with the supplied data.")
    }
    post_components <- post_components[keep]
    weights         <- weights[keep]

    # Re-weight components by marginal likelihood
    log_marg <- vapply(seq_along(post_components), function(i) {
      comp     <- components[keep][[i]]
      obs_mean <- if (type == "binary") x / n else x
      obs_se   <- if (type == "binary") {
        sqrt(obs_mean * (1 - obs_mean) / n)
      } else {
        (data_summary$sd %||% comp$fit_summary$sd) / sqrt(n)
      }
      stats::dnorm(obs_mean,
                   mean = comp$fit_summary$mean,
                   sd   = sqrt(comp$fit_summary$sd^2 + obs_se^2),
                   log  = TRUE)
    }, numeric(1))

    log_post_wts <- log(weights) + log_marg
    post_weights <- exp(log_post_wts - max(log_post_wts))
    post_weights <- post_weights / sum(post_weights)

    # Mixture posterior summary
    post_means <- vapply(post_components, function(p) p$fit_summary$mean, numeric(1))
    post_sds   <- vapply(post_components, function(p) p$fit_summary$sd,   numeric(1))
    mix_mean   <- sum(post_weights * post_means)
    mix_sd     <- sqrt(sum(post_weights * (post_sds^2 + (post_means - mix_mean)^2)))

    return(structure(
      list(
        dist        = "mixture",
        params      = list(weights = post_weights),
        components  = post_components,
        weights     = post_weights,
        method      = "posterior",
        expert_id   = prior$expert_id,
        label       = prior$label,
        input       = list(),
        fit_summary = list(
          mean = mix_mean,
          sd   = mix_sd,
          q025 = mix_mean - 1.96 * mix_sd,
          q500 = mix_mean,
          q975 = mix_mean + 1.96 * mix_sd
        )
      ),
      class = "bayprior"
    ))
  }

  # -- Generic fallback: Normal approximation ----------------------------------
  # No exact conjugate formula exists for this prior/data-type pairing (e.g.
  # a Beta prior with continuous data -- Beta is conjugate to Binomial, not
  # to a continuous likelihood). Previously this aborted, which contradicted
  # the compatibility warning shown to users elsewhere in the package, which
  # explicitly states the analysis will "proceed using a Normal
  # approximation." Generalises the same approach already used for Weibull
  # above (and the same approximation prior_conflict() uses for its
  # diagnostics): approximate both prior and likelihood as Normal via the
  # prior's fit_summary, and return a Normal posterior. This keeps the
  # promised approximation behaviour consistent between the numeric
  # diagnostics and any posterior-dependent plot or downstream computation,
  # including each component of a mixture prior.
  prior_mean <- prior$fit_summary$mean
  prior_sd   <- prior$fit_summary$sd

  if (is.null(prior_mean) || is.null(prior_sd) ||
      is.na(prior_mean)  || is.na(prior_sd)) {
    rlang::abort(glue::glue(
      "Conjugate update not implemented for dist='{prior$dist}' with ",
      "type='{type}', and no fit_summary mean/SD is available for a ",
      "Normal-approximation fallback."
    ))
  }

  if (type %in% c("poisson", "survival")) {
    obs_mean <- x / n
    obs_se   <- sqrt(x) / n
  } else if (type == "binary") {
    obs_mean <- x / n
    obs_se   <- sqrt(obs_mean * (1 - obs_mean) / n)
  } else {
    obs_mean <- x
    obs_se   <- (data_summary$sd %||% prior_sd) / sqrt(n)
  }
  obs_se <- max(obs_se, 1e-8)

  prior_var <- prior_sd^2
  lik_var   <- obs_se^2
  post_var  <- 1 / (1 / prior_var + 1 / lik_var)
  post_mean <- post_var * (prior_mean / prior_var + obs_mean / lik_var)

  rlang::inform(glue::glue(
    "No exact conjugate update exists for dist='{prior$dist}' with ",
    "type='{type}'; approximated the posterior as Normal(mean = ",
    "{round(post_mean, 4)}, sd = {round(sqrt(post_var), 4)})."
  ))

  .make_bayprior("normal", list(mu = post_mean, sigma = sqrt(post_var)),
                 "posterior (normal approximation)", prior$expert_id,
                 prior$label, list())
}