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.

Areal Data

Figure 1: Example of areal data.

Point-refereced Data

Figure 2: Example of point-referrenced data.

Fused Data

Figure 3: Example of fused 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

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

• GP components: A mean and a covariance functions defined as, respectively, \(\mathbb{E}[Z(\mathbf{s})] = m(\mathbf{s})\) and \({\rm Cov}(Z(\mathbf{s}_1), Z(\mathbf{s}_2)) = C(\mathbf{s}_1, \mathbf{s}_2)\).

• Common assumptions are stationarity and isotropy.

• Importance in spatial statistics: predominant foundation of geostatistical modeling.

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.

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.

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.

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).

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).

Beyond Traditional Modeling: 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.

Structure of an HGP

The HGP, denoted \(\mathrm{HGP}\{m(\mathbf{s}), v(\mathbf{s}), r(h)\}\), formulation depends on three functions:

• The mean and SD functions, defined as \[m(\mathbf{s}) = \mathbb{E}[Z(\mathbf{s})], \quad \text{and} \quad v(\mathbf{s}) = \sqrt{\mathrm{Var}(Z(\mathbf{s}))},\] respectively.

• Correlation function: \(r(h) = \mathrm{Cor}(Z(\mathbf{s}_1), Z(\mathbf{s}_2)),\) where \(h\) denotes the Hausdorff distance between \(\mathbf{s}_1, \mathbf{s}_2 \in \mathcal{B}(D)\).

• Induced covariance function: \(\mathrm{Cov}(Z(\mathbf{s}_1), Z(\mathbf{s}_2)) = v(\mathbf{s}_1) v(\mathbf{s}_2) r(h(\mathbf{s}_1, \mathbf{s}_2))\)

The \(v(\mathbf{s})\) function

A useful way to define this function is as follows: \[v(\mathbf{s}) = \exp \{ \alpha_0 + \alpha_1 w(\mathbf{s}) \},\] where \(w(\mathbf{s})\) is a covariate available for any \(\mathbf{s} \in \mathcal{B}(D)\).

• Special cases useful in practice:

  • Homoscedastic: \(w(\mathbf{s}) = 0\)

  • Data Fusion: \(w(\mathbf{s}) = \mathbb{1}(\mathcal{A}(\mathbf{s}) > 0)\).

  • Area dependent: \(w(\mathbf{s}) = \mathcal{A}(\mathbf{s})\)

• Although flexible, one has to be careful when choosing this function to ensure the process validity (Palacios and Steel 2006).

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

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

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.

• Next Up: Simulations and real-world applications will demonstrate the HGP’s robustness.

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

Outline

• Spatial GLMM and Bayesian Inference

• Simulation study for:

  • Areal data
  • Fused data

• Applications

  • Areal data and disease mapping
  • Fused data and air pollution

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\)

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

Bayesian Estimation

• Priors: \(\mathbf{Z} \sim \mathrm{HGP}\{0, v(\cdot), r(\cdot) \}\), \(\boldsymbol{\beta} \sim \mathcal{N}(\mathbf{0}, 10^2 \mathbf{I})\), \(\rho \sim \mathrm{Exp}(a_\rho)\), \(\nu = 0.7\).

  • Homoscedastic model (\(v(\mathbf{s}) = \sigma\)): \(\sigma \sim t_{+}(3)\)
  • Heteroscedastic model (\(\log(v(\mathbf{s})) = \alpha_0 + \alpha_1 w(\mathbf{s})\)): \((\alpha_0, \alpha_1)^\top \sim \mathcal{N}(\mathbf{0}, 10^2 \mathbf{I})\)

• 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)

  • \(\Delta_G = G(\text{M}) - G(\text{HGP})\), where \(G\) is a model comparison criterion and \(M\) is a competing method.

• 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})\)

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.

Respiratory Disease Hospitalization

• Sample units: 134 intermediate zones (IZ), where the \(i\)-th IZ is denoted \(\mathbf{s}_i\).

• Number of hospitalizations: \(\mathbf{Y} = {(Y(\mathbf{s}_1), \ldots, Y(\mathbf{s}_n))}^\top\).

• Expected number of hospitalizations based on the national age- and sex-standardized rates: \(E_i\).

• Percentage of people classified as income deprived: \(\mathbf{X}\).

\[\begin{aligned} & (Y(\mathbf{s}_i) \mid x_i, z(\mathbf{s}_i)) \sim \text{Poisson}(E_i \mu_i) \\ & \log(\mu_i) = \beta_0 + x_i \beta_1 + z(\mathbf{s}_i) \end{aligned}\]

Standardized Incidence Ratio (SIR)

Disease Mapping: GOF

Disease Mapping: Spatial Dependence

Disease Mapping: Insights

• HGP Advantages: The HGP provides a simpler and more intuitive understanding of spatial dependence compared to classical methods.

• Assumptions: The HGP’s core assumptions about spatial relationships differ from those used in specialized models for areal data.

• Spatial Correlation: In the HGP, spatial correlation decreases as the Hausdorff distance between areas increases.

• Specialized Models Complexity: Areal data models often rely on more complex spatial relationships based on detailed factors like shared borders and the number of neighboring units.

Data Analysis: Air Pollution in Ventura and Los Angeles

Air Pollution: GOF and 10-fold CV

Air Pollution: GOF and 10-fold CV

Air Pollution: Change-of-Support

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.

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” based on Hausdorff distance from neighbors, improving estimate reliability.
  • 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).

• The priors are the same as in the spatial setting. Here, we also placed a prior on \(\nu\).

\[\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

Forecast

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