Spatial Statistics

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

• General spatial model (Cressie 1993): \(\{ Z(\mathbf{s}) \; : \; \mathbf{s} \in D \}\), where \(D\) is an index set.

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

Classic Models for Spatial Data

• Areal data

  • CAR (Besag 1974), ICAR, BYM (Besag et al. 1991), DAGAR (Datta et al. 2019).

• Point-referenced data

  • Gaussian processes (GP), Nearest-neighbor mixture processes (Zheng et al. 2023)

• Mixed/fused data

  • “Aggregated” GP (AGP) (Moraga et al. 2017): \[ Z(\mathbf{s}_i) = \begin{cases} Z(\mathbf{s}_i), & \text{if } \mathcal{A}(\mathbf{s}_i) = 0 \\ {\mathcal{A}(\mathbf{s}_i)}^{-1} \int_{\mathcal{A}(\mathbf{s}_i)} Z(\mathbf{s}) \mathrm{d}\mathbf{s}, & \text{if } \mathcal{A}(\mathbf{s}_i) > 0 \end{cases} \]

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: An isotropic 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 (Appendix C, Molchanov 2005).

The HGP

• Index set: \(\mathcal{B}(D)\) represents the closed and bounded subsets of \(D \subset \mathbb{R}^2\).

• \(Z(\mathbf{s}) \sim \mathrm{HGP}\{m(\mathbf{s}), v(\mathbf{s}), r(h)\}\), where \(\mathbf{s} \in \mathcal{B}(D)\).

• Mean function: \(m(\mathbf{s}) = \mathbb{E}[Z(\mathbf{s})]\)

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

• SD function: \(v(\mathbf{s}) = \sqrt{\mathrm{Var}(Z(\mathbf{s}))}\)

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

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

We may defined: \[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)\).

• Useful special cases:

  • Homoscedastic: \(w(\mathbf{s}) = 0\) (consequence, \(\sigma = \exp \{ \alpha_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 the process’ validity

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 (PD): yields positive definite correlation matrices for all its finite-dimensional marginal distributions.

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

PD and Empirical Assessment

• Definition: A real function \(k\) on a metric space \((D, d)\) is PD if it is continuous and satisfies:

  \[ \sum_{i = 1}^{m} \sum_{j = 1}^{m} a_i a_j k(\mathbf{s}_i, \mathbf{s}_j) \geq 0 \]

  for real numbers \(a_1, \ldots, a_m\) and \(\mathbf{s}_1, \ldots, \mathbf{s}_m \in D\).

• Necessary condition: All the eigen values of the resulting correlation matrix must be non-negative.

• Empirical assessment: We can check for PD by examining the smallest eigenvalue of the correlation matrix across the parameter space.

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.

Visualizing the PEC function

Empirical Assessment of PEC’s 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.

• PEC Function Validity: A valid HGP depends on selecting an appropriate PEC correlation function.

• Empirical assessments: When analyzing a new dataset or trying a new correlation function, it is imperative to empirically assess the validity of the function.

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

Spatial GLMM

A generalized linear mixed effects 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, \boldsymbol{\gamma}) \\ & g(\mu_i) = \mathbf{x}_i \boldsymbol{\beta} + Z(\mathbf{s}_i). \end{aligned}\]

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

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

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

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

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

Priors

• Independent normal priors for the regression coefficients: \(\boldsymbol{\beta} \sim \mathcal{N}(\mathbf{0}, 10 \mathbf{I})\)

• HGP prior for the latent random effects: \(\mathbf{Z} \sim \mathrm{HGP}\{0, v(\cdot), r(\cdot) \}\)

• Exponential prior for the spatial dependence parameter: \(\rho \sim \mathrm{Exp}(a_\rho)\), where \(a_{\rho} = - \log(p_{\rho}) / \rho_0\).

  • \(a_\rho\) is chosen such that \(\mathbb{P}(\rho \geq \rho_0) = p_\rho\).

• Homoscedastic variance: \(v(\mathbf{s}) = \sqrt{\mathrm{Var}(Z(\mathbf{s}))} = \sigma \sim t_{+}(3)\)

• Heteroscedastic HGP: \(\alpha \sim \mathcal{N}(\mathbf{0}, \mathbf{I})\), where \(v(\mathbf{s}) = \exp \{ \alpha_0 + \sum_i \alpha_i w_i(\mathbf{s}) \}\).

Bayesian Inference & Model Assessment

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

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

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

• Goodness-of-fit criteria: LOOIC (lower values indicate better fit)

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

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

Respiratory Disease Hospitalization in Glasgow

Respiratory Disease Hospitalization in Glasgow

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

• 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: \(x_i\).

Disease Mapping: Estimation & GOF

Disease Mapping: Spatial Dependence

Air Pollution in Ventura and Los Angeles

Air Pollution: Change-of-Support

Applications: key findings

• The proposed method consistently demonstrated performance comparable to specialized models tailored for areal and fused data.

• Unlike traditional areal models, the HGP’s marginal variances are independent of the number of neighbors.

• The HGP model simplifies data fusion by bypassing the need to define arbitrary grids for numerical integral evaluation, eliminating this step each time the joint probability distribution of the data and parameters is calculated.

• Across both applications, the HGP provides an interpretable spatial dependence parameter and a spatial correlation function.

Highlights

The HGP has proven to be a powerful model that offers:

• Versatility: accomodates diverse spatial data types.

• Performance: competitive against models designed for specific spatial data types.

• Reliable predictions: prediction intervals with near nominal frequentist coverage.

• Our conclusions are further supported by a comprehensive simulation study, detailed in our available preprint.

Future work and Limitations

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