Fit the BYM2 expected-count model

Estimates expected counts for each district-month using a negative binomial or Poisson model with BYM2 spatial effects, fitted with INLA.

Usage

spi_expected(
  cases,
  population,
  adjacency,
  covariates = NULL,
  season = c("harmonic", "rw2", "monthly", "none"),
  year_effect = c("iid", "none", "rw1"),
  overdispersion = c("nb", "iid", "none", "auto"),
  prior_phi = list(U = 0.5, alpha = 0.5),
  prior_precision = list(U = 1, alpha = 0.01),
  prior_precision_year = list(U = 1, alpha = 0.01),
  n_draws = 1000L,
  log_transform = NULL,
  id_col = "district_id",
  pop_col = "pop_u15",
  keep_draws = TRUE,
  verbose = TRUE,
  check = TRUE,
  debug = FALSE,
  seed = 42L,
  num_threads = if (is.null(seed)) NULL else "1:1"
)

Arguments

cases: Tibble with columns: the district identifier (see id_col, character), month (Date), count (integer). One row per district-month. Zero-count rows must be present explicitly.
population: Tibble with the district identifier and pop (double). The time unit is inferred: if month (Date) is present, population is joined month-by-month and the offset is log(pop); otherwise year (integer) is required and the offset is log(pop / 12).
adjacency: nb object OR sf object. If sf, adjacency is computed internally via spi_adjacency() using id_col. Values of the id column must match those in cases.
covariates: Tibble with columns: the district identifier (see id_col), month (Date) OR year (int), plus one or more numeric covariate columns. If NULL, the model uses population, spatial, seasonal, and year terms according to the selected settings. Covariates are standardised internally (mean = 0, sd = 1). Covariate adjustment changes the expected counts and therefore the comparison made by SPI. For example, including travel time makes the expected count depend partly on access to care. Coefficients describe associations after accounting for the other model terms; they do not establish cause and effect. Default: NULL.
season: Character. Seasonal specification: "harmonic" (1st + 2nd order sin/cos, 4 terms), "rw2" (cyclic 2nd-order random walk, 12 knots), "monthly" (12 monthly fixed effects, January omitted), "none" (no seasonal component). Default: "harmonic".
year_effect: Character. Optional year effect to account for differences in reported counts between years (e.g. changes during insecurity, programme transitions, or a pandemic). "iid" (default) estimates each year's effect independently; this is the setting used in the accompanying study. It allows reporting to rise and fall between years without assuming a smooth trend. "none" omits the year effect; "rw1" smooths effects across adjacent years.
overdispersion: Character. How to allow for count variation beyond a Poisson distribution: "nb" (negative binomial likelihood, default, the paper specification), "iid" (Poisson-lognormal, iid N(0, sigma^2) on log scale per district-month), "none" (plain Poisson, not recommended for sparse data), or "auto" (fit all three and select a model using spi_compare_overdispersion(), then refit it at the requested n_draws). Default: "nb".
prior_phi: Named list with elements U and alpha giving the BYM2 mixing parameter penalised-complexity (PC) prior P(phi < U) = alpha. Default: list(U = 0.5, alpha = 0.5) (50% prior probability that phi is below 0.5).
prior_precision: Named list with elements U and alpha giving the marginal SD PC prior P(1/sqrt(tau) > U) = alpha. Default: list(U = 1, alpha = 0.01) (1% chance SD exceeds 1 on log scale).
prior_precision_year: Named list with elements U and alpha giving the PC prior for the year random effect precision. Default list(U = 1, alpha = 0.01) matches the BYM2 spatial-precision prior: a 1% prior probability that the log-scale SD exceeds 1. Ignored when year_effect = "none".
n_draws: Integer. Samples from the fitted model, used to calculate credible intervals for SPI. Each draw samples the model terms together. Use 1000 for analysis, 100 for quick checks. Default: 1000L.
log_transform: Character vector. Covariate column names to log(1 + x) transform. Typical: c("facility_index", "conflict_events"). Default: NULL.
id_col: Character. Name of the district identifier column in cases, population, and covariates. Must match the id_col used when building adjacency via spi_adjacency(). The function renames this column to district_id internally and renames it back on output. Default: "district_id".
pop_col: Character. Name of the denominator column in population. For AFP surveillance the at-risk population is under-15, so the default is "pop_u15". Pass "pop" if your tibble carries total population, "pop_u16" for under-16, or any other column name to switch population columns without renaming them first. The function uses pop internally and reports an error if the requested column is missing. Default: "pop_u15".
keep_draws: Logical. If TRUE (default), the full n_draws x n_district_months draws matrix is returned in result$draws. Set FALSE for large jobs to keep only the summary quantiles and save memory.
verbose: Logical. Print progress messages. Set FALSE to suppress them. Default: TRUE.
check: Logical. Run spi_check_inputs() on cases, population, and adjacency before fitting and abort on error-level issues. Default TRUE. Set FALSE only if the inputs have already been checked. Automatic model selection and spi_compare_overdispersion() use this to avoid repeating the same check for each fit.
debug: Logical. If TRUE, runs INLA in verbose mode (prints its raw stdout/stderr, including VB-correction notes), prints model diagnostics (CPO failure rate, hyperparameter posterior summary, sanity check that median expected count is comparable to median observed count), and skips muting INLA's compiled-binary messages. Use to investigate unexpected results or convergence problems. Default: FALSE.
seed: Integer. Random seed forwarded to inla.posterior.sample(), which controls the posterior sampling. Must be a single non-negative whole number; NULL draws a fresh seed each call. Seeding alone is not sufficient for run-to-run reproducibility – see the Reproducibility section. Default: 42L.
num_threads: Threads for the INLA fit, passed to INLA::inla() and inla.posterior.sample() as num.threads (INLA's "A:B" form, e.g. "1:1" or "12:1"). Defaults to "1:1" whenever seed is set, because a multithreaded fit is not bit-reproducible even at a fixed seed; pass NULL to inherit INLA::inla.getOption("num.threads") and use multiple threads, with some variation between runs. Ignored when seed = NULL unless given explicitly.

Value

Object of class spi_expected. A list containing:

draws: Matrix [n_draws x n_district_months] of posterior draws of exp(eta), with colnames of the form "<id>|YYYY-MM". NULL if keep_draws = FALSE.
summary: Tibble with the id column, month, count, pop, expected_median, expected_q05, expected_q10, expected_q90, expected_q95, expected_mean.
formula: INLA formula used.
model: Fitted INLA object.
adjacency: spdep nb object.
hyperparameters: Tibble with posterior summaries of tau, phi, sigma.
priors: List of prior settings used.
cov_params: Named list of (mean, sd) used to standardise each covariate, used by predict().
data: Input data tibble (for downstream functions).
id_col: The id column name, echoed for downstream use.
call: Matched call.

Details

Model specification. The expected count model uses negative binomial (default) or Poisson regression with:

Log person-time offset
BYM2 spatial random effects (Riebler et al. 2016)
Seasonal and year effects, controlled by season and year_effect
Optional covariates

Input checks. Before fitting, the function checks:

Required columns present and correct types
No negative counts or NAs
All district ids in cases appear in adjacency
No duplicate district x month combinations
Population coverage check
Adjacency graph connectivity check
Covariate validation

Computational notes.

Posterior sampling can be time-consuming for large datasets
Use fewer draws for initial exploration
INLA uses sparse matrix methods for efficiency

INLA diagnostic messages (debug = TRUE).

vb.correction aborted / iterative process seems to diverge – INLA's variational Bayes correction is an optional refinement on top of the Laplace approximation. When it fails to converge, INLA may return the Laplace approximation. Review convergence and other diagnostics before relying on the fit; this message alone does not establish whether the approximation is adequate.
Matrix is not positive definite during fitting – usually means the BYM2 precision matrix is singular; check that the adjacency graph is symmetric and that disconnected components are handled (see spi_adjacency()).

CPO / PIT and overdispersion = "iid". When the model has an iid effect per observation, conditional predictive ordinate (CPO) and the probability integral transform (PIT) can be unreliable. INLA may flag many observations as failure = 1. This can follow from the model specification and does not by itself establish poor fit. For CPO-based diagnostics, use the default overdispersion = "nb" or overdispersion = "none".

Paper specification

The accompanying study uses season = "harmonic", year_effect = "iid", overdispersion = "nb", with no covariates – the defaults below. The published index, spi_index(), is built on fits with these defaults.

Choosing a seasonal specification

Surveillance counts often show within-year cycles driven by transmission biology, climate, school terms, or reporting patterns. A seasonal term that is too restrictive can leave seasonal variation unexplained; one that is too flexible can fit random variation.

Quick diagnostic. Before fitting, plot the monthly average count to see whether a cycle exists and what shape it has:


cases |>
  dplyr::mutate(m = lubridate::month(month)) |>
  dplyr::group_by(m) |>
  dplyr::summarise(mean = mean(count, na.rm = TRUE)) |>
  plot(type = "b")

How to choose:

"none" – no visible cycle, or counts look flat across months. Use for aseasonal outcomes (e.g. neonatal tetanus, DHIS2 reporting completeness when administratively driven).
"harmonic" (default, 4 terms: sin/cos at 12 and 6 month periods) – one or two smooth peaks per year, using four coefficients.
"rw2" – cyclic 2nd-order random walk over 12 months. Smooth but arbitrary shape; lets the prior do the smoothing rather than imposing a sinusoid. Use when seasonality is real but asymmetric or multi-modal (e.g. cholera in some settings).
"monthly" – 11 free monthly fixed effects (January as reference). Estimates each month separately, without smoothing. It needs enough data (about three or more years across most districts) to distinguish monthly patterns from random variation.

Compare plausible choices. Examine the fitted seasonal pattern, diagnostics, and DIC or WAIC together. Small differences in a single criterion may not distinguish the models clearly.

Disease-specific starting points (acute case-based surveillance):

AFP / polio	harmonic
Measles	harmonic or rw2
Meningitis (African belt)	harmonic
Yellow fever	harmonic
Cholera	rw2 or monthly
Neonatal tetanus	none
DHIS2 reporting completeness	none or rw2

These are starting points, not prescriptions – always verify against the monthly-average plot for your data and country, since seasonality varies by climate zone and surveillance system.

Reproducibility

Reproducible results depend on three settings:

Which posterior configuration each draw comes from. inla.posterior.sample() picks this with R's own RNG, so set.seed() governs it.
The latent field drawn within that configuration. This comes from an RNG inside INLA's compiled code that R's set.seed() cannot reach, so the seed has to be handed to INLA directly.
The fit itself. INLA::inla() is not bit-reproducible when it runs multithreaded: the hyperparameter mode and integration points shift with thread scheduling, so the draws are taken from a slightly different posterior. num_threads governs this.

spi_expected() handles (1) and (2) whenever seed is non-NULL, restoring the caller's RNG state afterwards. (3) is why num_threads defaults to "1:1" alongside a seed.

Changes in these settings can affect SPI values and field-guide labels. A district's label can also change when results for other districts change, because extent depends on their SPI values.

The defaults use one thread for reproducibility. To use INLA's global thread setting, pass num_threads = NULL. Set seed = NULL to use a new seed on each call.

Small numerical differences remain possible with one thread and a fixed seed. Repeated fits have differed by up to 1.6e-6 in relative terms, much less than the differences observed with an unset seed.

References

Riebler A, et al. (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research, 25(4), 1145-1165.

Examples

if (FALSE) { # \dontrun{
data(synth_surveillance, package = "spi")

fit <- spi_expected(
  cases = synth_surveillance$cases,
  population = synth_surveillance$population,
  adjacency = synth_surveillance$boundaries,
  covariates = synth_surveillance$covariates,
  season = "harmonic"
)

print(fit)
} # }