Too many slides…

Tuberculosis in context

• Mortality: Tuberculosis (TB) remains a major global health threat, second in infectious disease mortality only to COVID-19 (as of 2023).

• Over 10.1 million peoplewere diagnosed and 1.23 million lost their life in 2024 alone due to TB ¹ adding to a historic legacy of over a billion deaths in the last two centuries.

• The WHO estimates that 22 billion dollars are required per year for prevention, diagnosis, and treatment in 2024 (but only 5.9 billion available).

• The burden is deeply unequal: while high-burden nations face crises, most countries in the Global North maintain incidence rates below 10 cases per 100,000.

TB in RS

• In Brazil: 40.4 cases per 100,000 inhabitants in 2024.

• Rio Grande do Sul (RS) historically reported high incidence

  • 38.7 to 45.6 (from 2011 to 2022)

• Dependence: Studies demonstrate strong spatial dependence of TB infections in Brazil, but temporal and spatiotemporal structures have been largely overlooked.

• Risk Factors: TB risk factors include densely populated areas, poverty, substance abuse, and incarceration (Cortez et al. 2021).

Disease mapping models

• Nonseparable SPT models are more complex as they consider that the spatial and temporal correlations might be intertwined (Cressie and Wikle 2015, pg. 309–321).

• Separable models encode the SPT covariance as the product of a temporal and a spatial terms.

  • Spatial models: CAR (Besag 1974), ICAR, BYM (Besag et al. 1991), Leroux (Leroux et al. 2000), DAGAR (Datta et al. 2019), RENeGe (Cruz-Reyes et al. 2023).

• Advantages of separable models: Computational efficiency & positive-definiteness of the covariance function.

Methodology & Objectives

• set-indexed Gaussian Process (GP): we propose using the newly developed set-indexed GPs for the spatial portion of the model (Godoy et al. 2026a; b) along with some computational improvements.

• Reliable incidence estimates:

  • Smaller municipalities benefit from borrowed strength from neighbors, improving estimate reliability.
  • Results enable the calculation of standardized incidence ratios to pinpoint high-risk areas.

• Forecasting: Predicted TB incidence rates one year ahead offer crucial insights for proactive public health planning.

Challenges

1. Quantifying isotropic spatial covariance

2. Validity of the proposed models

3. Computation: Given the high computational complexity of a “full” GP, we rely on a NNGP (Datta et al. 2016)

4. How to compare models in this context?

RF in spatial statistics

• A random field (RF): \(\{ \omega(\mathbf{s}) \; : \; \mathbf{s} \in D \}\), where \(D\) is an index set.

• RF are used extensively in spatial statistics (Cressie 1993), where sample units are conceptualized as elements of an index-set \(D\).

• RF assumptions change heavily according to the spatial structure/geometry of the observed spatial data.

A different look into random fields

• A more flexible index set: The class of class of non-empty, closed and bounded sets in \(D\) (denoted \(\mathcal{C}_D\)).

  • \(\{\omega(\mathbf{s}) \, : \, \mathbf{s} \in \mathcal{C}_D\}\) (Godoy et al. 2026b).
  • Isotropic GPs are defined with their covariance function depending on a distance between sets.

• Successful in practice: Competitive with areal models, often better than models for data fusion (Godoy et al. 2026b).

• Lacking theoretical foundation: no formal proof of validity of covariance functions.

Distances: from points to sets

• Hausdorff distance \[ h(A_1, A_2) = \inf \{ r \geq 0 \, : \, A_1 \subseteq {\rm B}_r(A_2), A_2 \subseteq {\rm B}_r(A_1) \}, \] where \(A_1 \subset D\) and \(A_2 \subset D\) are two non-empty sets.

• Ball-Hausdorff distance: Let \((D, d)\) be a length-space. The distance is defined as: \[ bh(A_1, A_2) = d(c(A_1), c(A_2)) + \lvert R(A_1) - R(A_2) \rvert, \] where \(A_1, A_2 \in \mathcal{C}_D\).

What is a ball?

Balls on spheres?

Theoretical Validity

• Distance Validity: We proved that the ball-Hausdorff distance (\(bh\)) is Conditionally Negative Definite (CND) on length spaces.

• Hilbert Embedding: The semi-metric space \((\mathcal{C}_D, \sqrt bh)\) can be isometrically embedded in a Hilbert space.

• Positive Definiteness: By Schoenberg’s Theorem, this ensures that Matérn and Powered Exponential covariance families are valid (positive-definite) for sets (under a \(\sqrt bh\) distance).

Valid covariance functions for sets

• The Powered Exponential (PEXP) covariance function \[ K(d; \, \theta) = \sigma^2 \exp \left\{ - \left( \frac{d}{\phi} \right)^{\nu} \right\} \] is a valid family on \((\mathcal{C}_D, \sqrt bh)\) for \(\nu \in (0, 2]\).

• The Matérn covariance function \[ K(d; \, \theta) = \sigma^2 \frac{1}{2^{\nu - 1}\Gamma(\nu)} {\left(\frac{d}{\phi}\right)}^{\nu} K_{\nu} \left( \frac{d}{\phi} \right) \] is a valid family on \((\mathcal{C}_D, \sqrt bh)\).

Computational bottlenecks

• The ball-Hausdorff distance is computationally efficient and needs to computed only once.

• Bottleneck: For \(N\) municipalities, a spatial covariance matrix \(\mathbf{K} = \sigma^2 \mathbf{R}_s\) requires \(\mathcal{O}(N^3)\) operations.

• NNGP Acceleration: Using nearest-neighbor Gaussian processes to approximate \(\mathbf{R}_s\), we are able to compute Cholesky of \(\mathbf{R}^{-1}_s\) efficiently reducing the complexity to \(\mathcal{O}(M^3 N)\), where \(M\) is the number of nearest neighbors.

Data

• Sample units: 497 municipalities, across 11 years (2011 to

• Number of TB cases: \(Y_t(\mathbf{s}_i)\) at location \(\mathbf{s}_i\) and time \(t\)

• Population: \(P_t(\mathbf{s}_i)\) (ranging from 1155 to 1.5M)

• Three covariates: penitenciary binary variable, homicide rate and population density (\(\mathbf{X}_t(\mathbf{s}_i)\))

Rates over time

Model

\[\begin{align*} Y_t(\mathbf{s}_i) & \mid \mathbf{X}_{t}(\mathbf{s}_i), \psi_t(\mathbf{s}_i) \overset{{\rm ind}}{\sim} f\big(\cdot \mid P_t(\mathbf{s}_i) \mu_t(\mathbf{s}_i), \boldsymbol{\tau}\big) \\ \mu_t(\mathbf{s}_i) & = \mathbb{E}[Y_t(\mathbf{s}_i) \mid \mathbf{X}_{t}(\mathbf{s}_i), \psi_t(\mathbf{s}_i)] \\ & = g^{-1}\big(\alpha + \mathbf{X}_{t}(\mathbf{s}_i)^\top \boldsymbol{\beta} + \psi_t(\mathbf{s}_i)\big) \end{align*}\]

• Expected incidence rate: \(\mu_t(\mathbf{s}_i)\)

• Latent random effect: \(\psi_t(\mathbf{s}_i)\)

• \(f(\cdot \mid \cdot)\) a pmf, possibly with additional parameters \(\mathbf{\tau}\).

• Intercept & regression coefficients: \(\pi(\alpha) \propto 1\) & \(\boldsymbol{\beta} \sim \mathcal{N}(\mathbf{0}, 100^2 \mathbf{I})\)

• Additional parameters: default priors from INLA.

Likelihood

• Poisson

• Negative Binomial

• Zero-inflated Poisson

• Zero-inflated negative binomial

Random effects

• Temporal dependence: \(\mathbf{R}_t(\gamma)\) is the correlation matrix of an AR(1) process, using a PC prior (Sørbye and Rue 2017) where \(\Pr(\lvert \gamma \rvert > 0.8) = 0.1\).

• Marginal variance: PC prior (Simpson et al. 2017) such that \(\Pr(\sigma > .05) = .05\).

Spatial component of the model

• We consider the following classic models as our baseline: ICAR, BYM (Besag et al. 1991), Leroux (Leroux et al. 2000)

  • using the (INLA) default hyperpriors for their parameters

• In addition, we employ the proposed set-indexed GP with a power exponential covariance using the NNGP for faster inference.

• Given the difficulty in estimating \(\nu\) (the power parameter), we fit models with \(\nu \in \{0.5, 1, 1.5, 2\}\).

• To assess the sensitivity to the choice of the number of nearest neighbors, we fit our model with \(M \in \{5, 10, 15, 20, 25\}\).

• In all cases, the prior on the range (\(\rho\)) was set so that: \(\Pr(\rho > \rho_0) = 0.01\).

Computational considerations

• Separable Bottleneck: For \(N\) municipalities and \(T\) years, the full SPT covariance \(\sigma^2 \mathbf{R}_s \otimes \mathbf{R}_t\) requires \(\mathcal{O}(N^3 T^3)\) operations.

• Kronecker Exploitation: Using properties of the Kronecker product, complexity drops to \(\approx \mathcal{O}(N^3 + T^3)\).

• Results relying on NNGP: High-dimensional spatiotemporal inference becomes feasible in few minutes under a complexity of \(\approx \mathcal{O}(M^3 N + T^3)\).

Bayesian Inference & Model Assessment

• Posterior: \(\pi(\boldsymbol{\theta} \mid \mathbf{y}, \boldsymbol{\psi}) \propto p(\mathbf{y} \mid \boldsymbol{\psi}, \boldsymbol{\theta}) p(\boldsymbol{\psi} \mid \boldsymbol{\theta}) \pi(\boldsymbol{\theta})\)

• INLA: Van Niekerk et al. (2023), Van Niekerk and Rue (2024)

• Parameter estimates: Quantiles of the posterior distribution

• Posterior predictive distributions: \(p(\mathbf{y}^{\ast} \mid \mathbf{y})\)

• Model assessment: Classic model assessment metrics (LPML) & approximate leave group out CV (Adin et al. 2024; Liu et al. 2025)

Modeling Strategy

1. Select a random effect structure

2. Select a likelihood

3. Assess sensitivity to the number of neighbors

4. Comparison across spatial models

Summary of model assessment

• The separable spatiotemporal structure was the best across the board. Regarding the likelihood, Poisson was the best.

• For the number of neighbors \(M\):

  • Smallest \(M\) yields best models according to LPML, WAIC, and DIC;
  • Largest \(M\) yields best models according to the leave group out CV (different group structures considered)

• Spatial structure:

  • set-GP with low smoothness (\(\nu = 0.5\)) was the best according to LPML, WAIC, and DIC
  • set-GP with high smoothness (\(\nu = 2.0\)) was the best according to leave group out CV

Spatial dependence

Relative Risks

Incidence: Small municipalities

Incidence: Large municipalities

Raw vs adjusted

Closing remarks

• Theoretical Foundation: Mathematically sound spatial covariance structures.

• Computational Scalability: High-dimensional spatiotemporal inference feasible in minutes.

• Flexibility: Models allow for tunable smoothness and is robust to changing spatial boundaries over time.

• Public Health Impact: More reliable risk factor estimates and out-of-sample forecasting to inform proactive TB control policies.

• Performance: Our approach matches or outperforms classic areal models (CAR, BYM, Leroux) while offering a more general framework.

Is there a package?

Model structure assessment

\(M\) Assessment

Comparison of spatial models

Sensitivity analysis I

Sensitivity analysis II

Sensitivity analysis III