Posterior Exceedance Probabilities from a Spatiotemporal Bayesian Model

Given a fitted climasus_spacetime_bayes object, computes for each municipality $\times$ time cell the posterior probability that the relative risk (RR) exceeds one or more critical thresholds: $P(RR > t \mid \text{data})$.

Exceedance probabilities are a cornerstone of Bayesian disease surveillance because they translate uncertain posterior estimates into decision-relevant quantities. For example, $P(RR > 1.5) > 0.80$ can trigger an alert for a municipality in a given week, whereas the raw posterior mean alone may understate uncertainty.

Usage

sus_mod_spacetime_exceedance(
  fit,
  thresholds = c(1, 1.5, 2),
  aggregate_time = NULL,
  municipalities = NULL,
  lang = "pt",
  verbose = TRUE
)

Arguments

fit: A climasus_spacetime_bayes object produced by sus_mod_spacetime_bayes(). Must contain a $rr data frame with columns code_muni, time_idx, rr_mean, rr_lower95, rr_upper95.
thresholds: Numeric vector of positive RR thresholds for which $P(RR > t)$ is computed. Default c(1.0, 1.5, 2.0).
aggregate_time: Character or NULL. If "year" or "month", aggregate the temporal dimension before reporting exceedances. Requires that time_idx be coercible to Date (formats "YYYY-MM-DD" or "YYYY-MM"). Default NULL (no aggregation).
municipalities: An sf object with at least a code_muni column and geometry. If provided, geometry is joined to the exceedance table so the result can be mapped directly. Requires the sf package. Default NULL.
lang: Character. Language for CLI messages: "pt" (default), "en", or "es".
verbose: Logical. If TRUE (default), print progress messages via cli.

Value

A named list of class c("climasus_spacetime_exceedance", "list") with the following slots:

$exceedance: data.frame with columns code_muni, time_idx (or time_period if aggregated), rr_mean, and one p_gt_X column per threshold (e.g. p_gt_1_0, p_gt_1_5, p_gt_2_0).
$thresholds: Numeric vector of thresholds used.
$n_exceed: data.frame with columns threshold and n_cells_exceed counting how many cells have $P(RR > t) > 0.80$.
$call: The matched call.

Approximation method

Exact computation requires access to the full MCMC sample or INLA posterior marginals. When only the posterior mean and 95\ stored (as in the $rr slot), exceedance probabilities are approximated via a log-normal approximation:

$$ \sigma_{\log} = \frac{\log(\text{RR}_{\text{upper95}}) - \log(\text{RR}_{\text{lower95}})}{2 \times 1.96} $$

$$ P(RR > t) \approx \Phi\!\left( \frac{\log(\text{RR}_{\text{mean}}) - \log(t)}{\sigma_{\log}} \right) $$

where $\Phi$ is the standard normal CDF. This is equivalent to assuming that $\log(RR)$ is approximately normally distributed in the posterior, which is a standard approximation for count-data hierarchical models (Richardson et al., 2004).

For exact computation when the full posterior is available, use INLA::inla.pmarginal() on the fitted INLA marginals directly.

Temporal aggregation

When aggregate_time is "year" or "month", the function averages $\log(RR_{\text{mean}})$ and $\sigma_{\log}$ across cells within each municipality $\times$ time period, then recomputes exceedance probabilities for the aggregated values. This is appropriate for surveillance summaries but not for causal inference.

References

Richardson, S., Thomson, A., Best, N., & Elliott, P. (2004). Interpreting posterior relative risk estimates in disease-mapping studies. Environmental Health Perspectives, 112(9), 1016-1025.

Lawson, A. B. (2018). Bayesian Disease Mapping: Hierarchical Modeling in Spatial Epidemiology (3rd ed.). CRC Press.

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). CRC Press.

Examples

if (FALSE) { # \dontrun{
library(climasus4r)

# Fit a spatiotemporal Bayesian model (requires climasus_spacetime_bayes)
# fit <- sus_mod_spacetime_bayes(df = agg, outcome = "deaths", ...)

# Compute exceedance probabilities for RR thresholds 1.0, 1.5, and 2.0
exc <- sus_mod_spacetime_exceedance(
  fit        = fit,
  thresholds = c(1.0, 1.5, 2.0),
  lang       = "pt"
)

print(exc)
head(exc$exceedance)
exc$n_exceed

# Aggregate to annual exceedances and join municipality geometry
library(geobr)
shp <- geobr::read_municipality(code_muni = "all", year = 2022)
exc_annual <- sus_mod_spacetime_exceedance(
  fit            = fit,
  thresholds     = c(1.0, 1.5, 2.0),
  aggregate_time = "year",
  municipalities = shp,
  lang           = "pt"
)
} # }