Spatial Data

• Taking into account spatial dependence possibly present in data is a foremost aspect of spatial statistics.

• Statistical inference depends heavily on the spatial structure/geometry of the observed spatial data.

General Spatial Model

The model below is flexible enough to analyze most spatial data structures (Cressie 1993)

\[\{ Z(\mathbf{s}) \; : \; \mathbf{s} \in D \},\]

where \(D\) is an index set.

• Areal models: index set is countable.

• Geostatistics: index set is a continuum.

• Gaussian Process (GP): a stochastic process such that all its finite-dimensional marginal distributions are multivariate normal distributions.

Research Questions

• The main research questions we are interested in are:

  1. Can we propose a model for spatial data that accomodates areal, point-referrenced, and fused data?

  2. If so, is this model competitive when compared to specialized models?

• Proposal: A GP defined on a flexible index set.

• Main challenge: Defining a valid correlation function.

Distances between Sets

• Metric space: \((D, d)\), where \(D\) is a spatial region of interest.

• Distance between a point and a set: \(d(x, A) = \inf_{a \in A} d(x, a)\), where \(d(x, y)\) is the distance between any two elements \(x, y \in D\).

• Directed Hausdorff distance: \[{\vec h}(A, B) = \sup_{a \in A} d(a, B)\]

• Hausdorff distance: \[h(A, B) = \max \left \{ \vec{h}(A, B), \vec{h}(B, A) \right \}\]

Hausdorff Distance for Spatial Data Analysis

• Study region: In spatial statistics, \(D\) is tipically a closed and bounded subset of \(\mathbb{R}^2\).

• In this context, the Hausdorff distance is a metric (Sendov 2004).

The HGP

• Isotropic stochastic process: \(\{ Z(\mathbf{s}) \, : \, \mathbf{s} \in \mathcal{B}(D) \},\) where \(D \subset \mathbb{R}^2\).

• Flexible index set: \(\mathcal{B}(D)\) contains all closed and bounded subsets of \(D\).

• The HGP: The stochastic process \(Z(\cdot)\) is an HGP if:

  • It is a GP

  • Its correlation function is a function of the Hausdorff distance between sets (in our case, spatial units).

• This process is suitable to analyze point-referrenced, areal, and mixed spatial data.

Ensuring Validity: Correlation Functions for the HGP

For a valid process, its correlation function must satisfy the following properties:

• Diminish with increasing distance: \(\lim_{h \to \infty}r(h) = 0\).

• Bounded from above by 1: \(r(0) = 1\).

• Positive-definiteness: yields positive-definite correlation matrices for all its finite-dimensional marginal distributions.

• Unfortunately, functions that are guaranteed to be positive definite on \((\lVert \cdot \rVert_2, \mathbb{R}^2)\) are not necessarily positive definite on other metric spaces (Li et al. 2023).

The Powered Exponential Correlation (PEC) Function

• PEC function: \[r(h; \phi, \nu) = \exp\left \{ - \frac{h^{\nu}}{\phi^{\nu}}\right \},\] where \(\nu\) is a smoothness parameter and \(\phi\) governs the range of the spatial dependence.

• Parametrization: We reparametrize this function with \(\rho = {\log(10)}^{1 / \nu} \phi\).

• Interpretation: \(\rho\) is the distance at which the spatial correlation reduces to \(0.10\).

Visualizing the PEC function

HGP Recap

• Flexibility: A process that handles point-referrenced, areal, and mixed spatial data by construction.

• Hausdorff distance: Enables HGP’s correlation function to account for the shape, size, and orientation of spatial objects.

• Validity: Using a PEC function ensures the HGP is a valid process.

flowchart LR
    A(Flexible index set) ==> B(GP)
    B ==> E(((HGP)))
    C(Hausdorff distance) ==> D(PEC)
    D ==> E

Notation

• Set of spatial locations: \(\mathcal{S} = \{ \mathbf{s}_1, \ldots, \mathbf{s}_n \}\), where \(\mathbf{s}_i \in \mathcal{B}(D)\) and \(D \subset \mathbb{R}^2\).

• Response variable: \(\mathbf{Y} = {(Y(\mathbf{s}_1), \ldots, Y(\mathbf{s}_n))}^\top\).

• Covariates: \(\mathbf{X} = (\mathbf{X}_1, \ldots, \mathbf{X}_p)\), where \(\mathbf{X}_{j} = (x_{j1}, \ldots, x_{jn})^\top\) and \(\mathbf{x}_i = (x_{1i}, \ldots, x_{pi})\).

• Spatial random effects: \(\mathbf{Z} = {(Z(\mathbf{s}_1), \ldots, Z(\mathbf{s}_n))}^\top\).

• Density functions (for response and latent variables): \(p(\cdot \mid \cdot)\)

• Prior and posterior distributions: \(\pi(\cdot)\) and \(\pi(\cdot \mid \cdot)\)

Spatial GLMM

A generalized linear mixed model (GLMM) can be written as \[\begin{aligned} & Y(\mathbf{s}_i) \mid \mathbf{x}_i, z(\mathbf{s}_i) \overset{{\rm ind}}{\sim} f(\cdot \mid \mu_i, \omega) \\ & g(\mu_i) = \mathbf{x}_i \boldsymbol{\beta} + z(\mathbf{s}_i). \end{aligned}\]

• Probability distribution: \(f(\cdot)\)

• Link function: \(g(\cdot)\)

• Conditional mean: \(\mu_i = \mathbb{E}[Y(\mathbf{s}_i) \mid \mathbf{x}_i, z(\mathbf{s}_i)]\)

• Model parameters: \(\theta = {\{\boldsymbol{\beta}^\top, \boldsymbol{\sigma}^\top, \delta^\top, \omega^\top \}}^\top\)

• Joint density: \(p(\mathbf{y} \mid \mathbf{z}, \theta) = \prod_{i = 1}^n f(y(\mathbf{s}_i) \mid \mu_i, \omega)\)

Bayesian Estimation

• Priors: \[\begin{align*} & \mathbf{Z} \sim \mathrm{HGP}\{0, v(\cdot), r(\cdot) \} \\ & \boldsymbol{\beta} \sim \mathcal{N}(\mathbf{0}, 10^2 \mathbf{I}) \quad \sigma \sim t_{+}(3) \\ & \rho \sim \mathrm{Exp}(a_\rho) \quad \nu \sim \mathrm{Beta}(a_\nu, b_\nu) \end{align*}\]

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

• MCMC sampler: No-U-Turn (Homan and Gelman 2014).

• Convergence assessment: traceplots and split-\({\hat{R}}\) (Vehtari et al. 2021).

Goodness-of-fit and Predictions

• Goodness-of-fit criteria: DIC, WAIC, and LOOIC (lower values indicate better fit)

• Unobserved outcome at \(m\) locations: \(\mathbf{Y}^\ast = (y(\mathbf{s}^\ast_1), \ldots, y(\mathbf{s}^\ast_m))^\top\).

• Posterior predictive distributions: \(p(\mathbf{y}^{\ast} \mid \mathbf{y}) = \int p(\mathbf{y}^{\ast} \mid \mathbf{z}^{\ast}, \theta) p(\mathbf{z}^{\ast} \mid \mathbf{z}, \theta) \pi(\theta \mid \mathbf{y}) {\rm d} \theta\)

• Monte Carlo samples from \(p(\mathbf{y}^{\ast} \mid \mathbf{y})\) are obtained using the draws from the posterior \(\pi(\theta \mid \mathbf{y})\)

• Predictions assessment: Interval Score (IS) and RMSP (lower values indicate better fit)

Tuberculosis in Brazil

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

• Brazilian TB Trends: Brazil’s national TB incidence declined only slightly from 2006–2015. Rio Grande do Sul (RS) reported significantly higher incidence than the national average in 2021, with the eastern region even more affected.

• Spatial Dependence: Studies demonstrate strong spatial dependence of TB infections (and co-infections) in Brazil, both nationwide and in specific cities. Time dependence has been largely overlooked in the literature.

• Risk Factors: TB risk factors include densely populated areas, poverty, substance abuse, and incarceration.

Objectives

• Reliable incidence estimates:
  • TB incidence across the study period, incorporating location-specific factors.
  • Smaller municipalities benefit from “borrowed strength” from neighbors (in the Hausdorff distance sense).
  • Results enable the calculation of SIRs to pinpoint high-risk areas.
• Forecasts:
  • TB incidence rates one year ahead.
  • These predictions offer crucial insights for proactive public health planning.

Model

• Sample units: 54 municipalities, across 11 years (2011 to 2021). We use 2022 to assess the quality of predictions.

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

• Population: \(P_t(\mathbf{s}_i)\).

• Five covariates and two way interactions between presence of prison and continuous variables (except IDESE).

\[\begin{aligned} & (Y_t(\mathbf{s}_i) \mid \mathbf{X}_{t}(\mathbf{s}_i), z_t(\mathbf{s}_i)) \overset{{\rm ind}}{\sim} \text{Poisson}(P_t(\mathbf{s}_i) \mu_{it}) \\ & \log(\mu_{it}) = \beta_0 + \mathbf{X}^\top_{t}(\mathbf{s}_i) \beta + z_t(\mathbf{s}_i) \end{aligned}\]

Spatiotemporal Trend

Explanatory Variables

Random Effects

• No random effects

• Unstructured: Random intercept per municipality (iid Normal)

• Time dependent: \(\mathrm{AR}(1)\)

• Spatial: Homoscedastic HGP

• Separable spatiotemporal: \(\mathrm{AR}(1)\) + Homoscedastic HGP, Multivariate \(\mathrm{AR}(1)\) with HGP errors

• Separable spatiotemporal and unstructured: Multivariate \(\mathrm{AR}(1)\) with HGP errors + random intercept per municipality

Results: GOF and Predictive Performance

Results: Relative Risks

Spatiotemporal Dependence

Small Municipalities

Forecast

Population at Risk in 2022

Highlights

• Tailored an HGP extension for spatiotemporal disease mapping.

• Competitive with specialized models

• It helps to gain insights into spatiotemporal disease mapping through spatiotemporal correlation functions.

• More reliable estimates of risk factors

• Out-of-sample predictions to inform public policies

Conclusion

The HGP offers a powerful and versatile model for spatial analysis. Advantages include:

• Common framework for diverse spatial data types.

• Competitive performance when compared to specialized models.

• Directly interpretable spatial parameters.

• Reliable and accurate predictions.

Future work and Limitations

• Extensions:
  • non-Euclidean spaces
  • Big data
• Limitations:
  • Difficulties in obtaining spectral densities for correlation functions
  • Big “n” problem inherited from geostatistics
  • Unclear how to incorporate anisotropy

Areal Models: Limitations

• Out-of-sample predictions: Non-trivial, it may be necessary to refit the model.

• Interpretability: Relationship between the spatial dependence parameter and correlation between neighbors is counterintuitive.

• Heteroscedasticity: Marginal SD depend on the number of neighbors.

• Extensions for change of support and data fusion problems are cumbersome and there exists no consensus in the literature in how to approach such problems.

Data Fusion Models: Limitations

• Computing: Although satisfactory approximations exist, it becomes computationally prohibitive for moderately large datasets.

• Approximations: Inferences may drastically change depending on the precision of the approximations.

• Biases: It is hard to quantify the impact of the discretizations on inferences (Gonçalves and Gamerman 2018).

Areal Data and the CAR Model

• \(Z_k \mid Z_{-k} \sim \mathcal{N}\left (\frac{\psi} {m_k} \sum_{i \sim k} Z_i, \frac{\tau^2}{m_k} \right)\)

• Special case: ICAR.

• Extensions of the ICAR model, providing a measure of spatial dependence have been developed (Dean et al. 2001; Leroux et al. 2000).

• The directed acyclic graph auto-regressive (DAGAR) model (Datta et al. 2019) proposes a new way to construct precision matrices using a DAG derived from the original adjacency matrix.

Fused Data

• GP at the point-referrenced level and a stochastic integral at the areal level (AGP) (Moraga et al. 2017): \[Z(\mathbf{s}_i) = \begin{cases} {\mathcal{A}(\mathbf{s}_i)}^{-1} \int_{\mathbf{s}_i} Z(\mathbf{s}) \mathrm{d}\mathbf{s}, & \text{if } \mathcal{A}(\mathbf{s}_i) > 0, \\ Z(\mathbf{s}_i), & \; \text{otherwise,} \end{cases}\]

• Covariance (Johnson et al. 2020): \[\mathrm{Cov}(Z(\mathbf{s}_i), Z(\mathbf{s}_j)) = {(\mathcal{A}(\mathbf{s}_i) \mathcal{A}(\mathbf{s}_j))}^{-1} \int_{\mathbf{s}_j} \int_{\mathbf{s}_i} C(\mathbf{x}, \mathbf{y}) \mathrm{d} \mathbf{x} \mathrm{d} \mathbf{y}\]

• Arbitrary precision: The stochastic integrals are evaluated numerically and depend on a user-defined grid.

Theoretical Foundation: Positive Definiteness of the PEC

• The proposition above guarantees the validity of the HGP equipped with a PEC function with \(\nu \in (1/2, 1)\).

• The proof is based on embedding the Hausdorff distance into a high-dimensional \(L_1\) normed Euclidean space, and using the fact that the exponential correlation function is positive definite on this space.

Practical Validation: Empirical Assessment of Positive Definiteness

Simulation Studies: Setup

• Goal: Evaluate HGP’s ability to adapt to different scenarios and compare its performance to specialized models.

• 200 simulated datasets for each combination of data generation process and scenario.

• Estimation assessment: MAPE, RMSE, frequentist coverage of the credible interval (CP).

• Predictions assessment: bias, RMSP, interval score (Gneiting and Raftery 2007, IS), frequentist coverate of the prediction intervals (CPP).

Simulation Studies: Scenarios

Simulation Study: Areal Data

• Data model: Poisson GLMM, single covariate, log link function, and a spatial random effect.

• We fit Poisson GLMMs with a homoscedastic HGP prior on the random effect.

• Assessment:

  • Parameter estimation when generating data from HGP random effects.
  • How good are the estimates of the regression coefficient under random effect misspecification?
  • Goodness-of-fit comparison with DAGAR priors on the random effect.

Areal Random Effects

• Data Generation:

  • Besag-York-Mollié (BYM, Besag et al. 1991): \(\mathbf{Z} \sim \mathcal{N}(\mathbf{0}, \sigma^{2}((1 - \zeta) \mathbf{I} + \zeta \mathbf{Q}^{-})\), where \(\mathbf{Q}\) is the precision matrix from the ICAR model and \(\zeta\) is a mixing parameter.
  • Aggregated GP (AGP): \(Z(\mathbf{s}_i) = {\mathcal{A}(\mathbf{s}_i)}^{-1} \int_{\mathbf{s}_i} K(\mathbf{s}) \mathrm{d} \mathbf{s}\), where \(\{ K(\mathbf{s}) \, : \, \mathbf{s} \in D \}\) is a (point-level) zero-mean, second order stationary and isotropic GP.

• Competing Method:

  • DAGAR: \(\mathbf{Z} \sim \mathcal{N}(\mathbf{0}, \Lambda(\tau, \psi))\), where \(\Lambda\) is a precision matrix, \(\tau\) is the conditional precision, and \(\psi\) measures spatial dependence.

Areal Simulation: Estimation

• Generating data from HGP:

  • CP of the hardest parameter to estimate (\(\rho\)) was \(92.5\%\) in the worst case scenario.
  • For all the other parameters, the coverage was extremely close to the nominal level.

• Regression parameter under misspecification:

  • Largest (worst) MAPE was \(6.6\%\).
  • CP of the credible interval suggests reliable estimates.

Areal Simulation: GOF

HGP outperforms DAGAR when simulating data from a AGP random effect.

Simulation Study: Data Fusion

• Data model: Gaussian GLMM, no covariates, identity link function, and a spatial random effect.

• We fit the models with a heteroscedastic HGP as prior on the random effect.

• Competing method: AGP₁ and AGP₂ priors on the random effect. The latter employes a finer grid/mesh resolution to approximate integrals numerically.

• Parameters estimation when generating data from HGP random effects.

• Model comparison: Predictions assessment.

Competing Method: Further Details

• The model employs a Matérn covariance function (at point level) with smoothness parameter fixed at 1. The other two parameters are the SD \(\sigma\) and the practical range \(\rho\). Both are analogous to the parameters with the same name in the HGP model.

• Inferences are conducted using the INLA (Rue et al. 2009) and the stochastic partial differential equation (Lindgren et al. 2011, SPDE) method.

Data Fusion Simulation: Estimation

• The HGP inherits same problems from geostatistical models.

• There are no weakly consistent estimators for the variance and spatial dependence parameters (Zhang 2004)

• Predictions are still optimal (Zhang 2004)

• The small scale variation parameter had a low coverage in the credible intervals (\(79.5\%\) for scenario 2).

• The MAPE for the \(\sigma\), \(\sigma_a\), and \(\rho\) parameters fluctuated around \(10\%\) in both scenarios.

Data Fusion Simulation: Prediction

• Uncertainty of the predictions: The HGP handles uncertainty in its predictions better than the AGP₁ model, even when the latter reflects the true data generation process.

• HGP as an AGP Approximation: The HGP delivers satisfactory approximations to an AGP model.

• HGP Advantages: The HGP offers potential computational advantages for large datasets while maintaining clear parameter interpretation.

Highlights

• Flexibility & Interpretability: Models diverse spatial data types and provides meaningful spatial parameters, regardless of unit size or shape.

• Predictions: The HGP delivers accurate predictions

• Areal Data: The model delivers valuable insights on the spatial dependence and accounts for sample units’ size, shape, and orientation.

• Data Fusion: The HGP simplifies the modeling process eliminating the need of a user-defined grid or mesh. Consequently, avoids computationally intensive numerical integrations with arbitrary precision, making inference seamless.