Bayesian Spatiotemporal Hierarchical Model with INLA

Fits a full Bayesian spatiotemporal hierarchical model using INLA (Integrated Nested Laplace Approximation). The model decomposes disease risk into structured spatial effects (BYM2, BYM, Besag, or IID), structured temporal effects (RW1, RW2, AR1, or IID), and optionally a space-time interaction following the four Knorr-Held (2000) interaction types. Climate and environmental covariates from a climasus_df can be included as fixed effects.

Usage

sus_mod_spacetime_bayes(
  df,
  outcome,
  W,
  time_col = "date",
  time_unit = c("year", "month", "week", "auto"),
  covariates = NULL,
  offset = NULL,
  family = c("poisson", "nbinomial", "binomial", "gaussian"),
  spatial_model = c("bym2", "bym", "besag", "iid"),
  temporal_model = c("rw1", "rw2", "ar1", "iid_time"),
  interaction_type = c("none", "I", "II", "III", "IV"),
  pc_prior_u = 0.5,
  pc_prior_alpha = 0.01,
  compute_waic = TRUE,
  compute_cpo = FALSE,
  exceedance_threshold = 1,
  n_samples = 1000L,
  seed = 42L,
  lang = "pt",
  verbose = TRUE
)

Arguments

df: A data.frame or climasus_df at the aggregate, spatial, or climate pipeline stage. Must contain code_muni (7-digit IBGE code), the outcome column, the time_col column, all columns named in covariates, and optionally an offset column. The data must form a panel (one row per area-time combination).
outcome: Character. Name of the response (count) column in df.
W: A climasus_weights object produced by sus_mod_spatial_weights(). Must contain either a $W dense matrix or a $nb neighbour list used to build the INLA graph.
time_col: Character. Name of the column containing time information. Can be a Date, integer, or character column. Default "date".
time_unit: Character. Temporal aggregation unit used to build the integer time index. One of "year", "month", "week", or "auto" (automatic detection based on date range). Default "auto".
covariates: Character vector of covariate column names in df to include as fixed effects. Default NULL (intercept-only model).
offset: Character. Name of a population or expected counts column in df (natural scale). Entered as log(offset) in the linear predictor. Default NULL.
family: Character. Likelihood family. One of "poisson" (default), "nbinomial" (negative binomial), "binomial", or "gaussian".
spatial_model: Character. Spatial random effect structure. One of "bym2" (default), "bym", "besag", or "iid". See Spatial model options.
temporal_model: Character. Temporal random effect structure. One of "rw1" (default), "rw2", "ar1", or "iid_time".
interaction_type: Character. Space-time interaction type following Knorr-Held (2000). One of "none" (default), "I", "II", "III", or "IV". See Knorr-Held interaction types.
pc_prior_u: Positive numeric. PC prior parameter $u$ for the BYM2 spatial precision: $P(\sigma > u) = \alpha$. Default 0.5.
pc_prior_alpha: Numeric in (0, 1). PC prior parameter $\alpha$ (tail probability). Default 0.01.
compute_waic: Logical. If TRUE (default), compute WAIC via control.compute = list(waic = TRUE).
compute_cpo: Logical. If TRUE, compute CPO/PIT values via control.compute = list(cpo = TRUE). Default FALSE (slower).
exceedance_threshold: Positive numeric. Compute $P(RR > threshold)$ for each space-time cell using the posterior marginals. Default 1.0.
n_samples: Positive integer. Number of posterior samples drawn via INLA::inla.posterior.sample() to compute RR credible intervals and exceedance probabilities. Default 1000L.
seed: Integer. Random seed for reproducibility. Default 42L.
lang: Character. Language for CLI messages: "pt" (default), "en", "es".
verbose: Logical. If TRUE (default), emit progress messages via cli.

Value

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

$fixed: data.frame with columns term, mean, sd, lower95, upper95, mode (posterior fixed-effect summaries).
$rr: data.frame with columns code_muni, time_idx, rr_mean, rr_lower95, rr_upper95, p_exceed (relative risk per space-time cell with exceedance probability).
$spatial_re: data.frame with columns code_muni, phi_mean, phi_sd (posterior summaries of structured spatial random effect).
$temporal_re: data.frame with columns time_idx, psi_mean, psi_sd (posterior summaries of structured temporal random effect).
$interaction_re: data.frame with columns code_muni, time_idx, gamma_mean, gamma_sd, or NULL if interaction_type = "none".
$fitted: Numeric vector of posterior mean fitted values (response scale).
$waic: Named list with $waic (scalar) and $p_eff (effective number of parameters), or NULL if compute_waic = FALSE.
$dic: Named list with $dic (scalar) and $p_d (effective parameters), or NULL.
$model_spec: Named list with $spatial_model, $temporal_model, $interaction_type, $family.
$call: The matched call.
$n_areas: Integer. Number of unique municipalities.
$n_times: Integer. Number of unique time periods.

Model structure

The linear predictor for area $i$ at time $t$ is:

$$\eta_{it} = \alpha + \mathbf{x}_{it}^{\top}\boldsymbol{\beta} + \phi_i + \psi_t + \delta_{it} + \log(E_{it})$$

where $\phi_i$ is the spatial random effect, $\psi_t$ the temporal random effect, $\delta_{it}$ the optional space-time interaction, and $E_{it}$ the expected counts offset.

Spatial model options

"bym2": BYM2 model (Riebler et al. 2016): scaled ICAR + IID with a single precision and a mixing parameter. Recommended for disease mapping. PC priors applied.
"bym": Original BYM model (Besag, York & Mollie 1991): ICAR + independent Gaussian components with separate precision parameters.
"besag": Intrinsic CAR (Besag 1974): purely structured spatial random effect.
"iid": Independent normal random effects (no spatial structure).

Knorr-Held interaction types

"none": No space-time interaction.
"I": IID interaction: $\delta_{it} \sim N(0, \sigma^2_\delta)$.
"II": Structured in time, unstructured in space: separate RW1 per area.
"III": Structured in space, unstructured in time: separate ICAR per time period.
"IV": Structured in both space and time (inseparable interaction).

PC priors

For the BYM2 spatial model, penalised-complexity (PC) priors are used for the precision hyperparameter: $P(\sigma > u) = \alpha$ where pc_prior_u and pc_prior_alpha control the tail probability. Default values u = 0.5, alpha = 0.01 are weakly informative and suitable for standardised disease mapping.

References

Knorr-Held, L. (2000). Bayesian modelling of inseparable space-time variation in disease risk. Statistics in Medicine, 19(17-18), 2555-2567.

Riebler, A., Sorbye, S. H., Simpson, D., & Rue, H. (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research, 25(4), 1145-1165.

Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models using integrated nested Laplace approximations. Journal of the Royal Statistical Society B, 71(2), 319-392.

Blangiardo, M., & Cameletti, M. (2015). Spatial and Spatio-temporal Bayesian Models with R-INLA. Wiley.

Besag, J., York, J., & Mollie, A. (1991). Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics, 43(1), 1-20.

Examples

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

# Build spatial weights for Rio Grande do Norte
shp <- geobr::read_municipality(code_muni = "RN", year = 2022)
W   <- sus_mod_spatial_weights(shp, return_matrix = TRUE)

# Prepare aggregated health-climate panel data (areas x time periods)
# agg <- sus_data_aggregate(cleaned_df) |>
#          sus_climate_inmet() |>
#          sus_climate_aggregate(temporal_strategy = "exact")

# BYM2 spatiotemporal model with Poisson family
result <- sus_mod_spacetime_bayes(
  df               = agg,
  outcome          = "deaths",
  W                = W,
  time_col         = "date",
  time_unit        = "month",
  covariates       = c("temp_mean", "prec_total"),
  offset           = "population",
  family           = "poisson",
  spatial_model    = "bym2",
  temporal_model   = "rw1",
  interaction_type = "I",
  pc_prior_u       = 0.5,
  pc_prior_alpha   = 0.01,
  compute_waic     = TRUE,
  n_samples        = 1000L,
  seed             = 42L,
  lang             = "pt"
)

print(result)
result$fixed
result$rr
result$spatial_re
result$temporal_re
result$waic
} # }