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.

• The model of choice 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.

Preliminaries: Isotropic GPs

• Gaussian process: A stochastic process such that all its finite-dimensional marginal distributions are Gaussian.

• Completely characterized (probabilistically) by two functions:

  • The mean function: \(m(\mathbf{s}_1) = \mathbb{E}[Z(\mathbf{s}_1)]\)
  • The covariance function: \[\mathrm{C}(Z(\mathbf{s}_1), Z(\mathbf{s}_2)) = \mathrm{Cov}(Z(\mathbf{s}_1), Z(\mathbf{s}_2)).\]

• A GP is isotropic if \(\mathrm{Cov}(Z(\mathbf{s}_1), Z(\mathbf{s}_2)) = \mathrm{C(d(\mathbf{s}_1, \mathbf{s}_2))}\), where \(d\) is a distance function.

Preliminaries: Sets and distances

• Point-referenced and areal spatial units are (closed and bounded) sets.

• We need to generalize distance between points to distance between sets.

• Ideally, this distance should:

  1. Take into account the shape, size, and orientation of spatial sample units.
  2. Be “spatially interpretable”.

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 non-empty, closed and bounded subset of \(\mathbb{R}^2\).

• In this context, the Hausdorff distance is a metric (Appendix C, Molchanov 2005).

The Hausdorff–Gaussian Process

• Index set: \(\mathcal{B}(D)\) represents the non-empty, 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})]\)

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

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

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

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: \(\rho = {\log(10)}^{1 / \nu} \phi\).

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

Visualizing the PEC function

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

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.

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

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)

Bayesian modeling recap

• HGP as a prior for the random effects disribution in a GLMM.

• Bayesian inference through MCMC.

• Uncertainty quantification of predictions through the posterior predictive distributions.

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

Disease mapping: Adjusted SIR

Air pollution in Ventura and Los Angeles

Air pollution: Estimation and prediction assessment

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 of the current work

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.

Work in progress (Theory)

• Radially positive-definite functions of the Hausdorff distance on non-Euclidean spaces(with Marcos Prates, Fernando Quintana and Jun Yan)

• HGP, big data and INLA (with Marcos Prates and Jun Yan)

• Nearest-neighbor mixture processes and the Hausdorff distance (with Bruno Sansó)

• Neural-networks and spatiotemporal data fusion (with Lucas Michelin, Marcos Prates, and Abhi Datta)

• Multiplier bootstrap and fast inference for auto-logistic models (with Yuan Huang and Jun Yan)

Work in progress (Application)

• Review paper on sports analytics (with Greg J. Matthews and Jun Yan)

• Combining Bayesian hierarchical models & population dynamics for species distribution models (with Malin Pinsky)

• R package to estimate risks in consumer automobile loans (with Jackson Lautier)

• Non-stationary geostatistical models for stock assessment in fisheries (with James Thorson)

• Non-linear mixed effect models to analyze plants heat-tolerance changes in face of climate change (with Georgia Hernández)

Other contributions

• Organization of weekly meetings and extracurricular projects for students.

• Co-advised a Master’s student in Brazil (Michelin et al. 2025)

• Many papers as a result of my work as statistical consultant. Highlights: Andreyeva et al. (2024) and Sieger et al. (2025)

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