Spatial data misalignment
Spatial data
In health and environmental applications, we often need to analyze data available at different spatial and spatio-temporal resolutions and that come from different sources
Spatial data misalignment
The analysis of spatial data at different spatial resolutions entails a number of statistical challenges. These may occur in several inference problems:
Data fusion
Better predict a variable by combining data available at several spatial resolutions
Estimate air pollution by combining point- and area-level data
Interpolation
Predict a variable at locations or areal units different from those of its original collection
Downscale health outcomes from state to municipality level
Regression
Relationship between response variable and explanatory variables at different spatial scales
Relationship between dengue at county level and temperature given at point locations
Gotway and Young, JAMA, 2002
Data fusion
Ground measurements (points) + Satellite derived measurements (grid)
European Environment Agency (EEAA) https://www.eea.europe.eu. NASA Socioeconomic Data and Applications Center (SEDAC) https://sedac.ciesin.columbia.edu
Fast and flexible spatial modeling by assuming a spatially continuous variable underlying all observations modeled using a Gaussian random field
Moraga et al., Spatial Statistics, 2017
Fast and flexible spatial modeling
Assume there is a spatially continuous variable underlying all observations that can be modeled using a zero-mean Gaussian random field \(S=\{S(\boldsymbol{x}): \boldsymbol{x} \in D \subset \mathbb{R}^2 \}\)
\[Y(\boldsymbol{x})|S(\boldsymbol{x}) \sim N(\mu(\boldsymbol{x})+S(\boldsymbol{x}),\tau^2)\]
Point data observed at \(\boldsymbol{x}_i \in D\) \[E[Y(\boldsymbol{x}_i)] = \mu(\boldsymbol{x}_i)+S(\boldsymbol{x}_i)\]
Areal data arise as region averages in regions \(B_j \subset D\) \[E[Y(B_j)]=|B_j|^{-1}\int_{B_j}(\mu(\boldsymbol{x})+S(\boldsymbol{x}))d\boldsymbol{x},\ |B_j|>0\] Moraga et al., Spatial Statistics, 2017
Fast and flexible spatial modeling
Models for point and areal obs. do not represent the same process:
- Include calibration parameter for the common spatial random effect:
\[E[Y(\boldsymbol{x}_i)] = \mu(\boldsymbol{x}_i) + S(\boldsymbol{x}_i)\]
\[E[Y(B_j)] = |B_j|^{-1} \int_{B_j} (\mu(\boldsymbol{x}) + \alpha(\boldsymbol{x}) S(\boldsymbol{x}))dx\]
- Different covariates and spatial random effects for each type of data:
\[E[Y(\boldsymbol{x}_i)] = \mu_1(\boldsymbol{x}_i) + S(\boldsymbol{x}_i) + S_1(\boldsymbol{x}_i)\]
\[E[Y(B_j)] = |B_j|^{-1} \int_{B_j} (\mu_2(\boldsymbol{x}) + \alpha(\boldsymbol{x}) S(\boldsymbol{x}) + S_2(x))dx\]
Inference using INLA + custom matrix SPDE
Integrated nested Laplace approximations (INLA) is a computational approach to perform approximate Bayesian inference in latent Gaussian models
In the SPDE approach, the continuously indexed Gaussian random field \(S\) is represented as a discretely indexed Gaussian Markov random field (GMRF) by means of a finite basis function defined on a triangulation of the study region
\[S(\boldsymbol{x}) = \sum_{g=1}^G \psi_g(\boldsymbol{x}) S_g\]
\(\psi_g(\cdot)\) piecewise polynomial basis functions on each triangle
\(\{S_g \}\) zero-mean Gaussian distributed
\(G\) number of vertices in triangulation
Point observations
\(S(\boldsymbol{x})\) weighted average of the GMRF values at the vertices of the triangle containing the point. Weights = barycentric coordinates
\[S(\boldsymbol{x}) \approx \frac{T_{1}}{T}S_1 + \frac{T_{2}}{T}S_2 + \frac{T_{3}}{T}S_3\]
\(T_1, T_2, T_3\) areas subtriangles formed by \(\boldsymbol{x}\) and vertices. \(T\) area whole triangle
Areal observations
\(S(B)=|B|^{-1} \int_{B} S(\boldsymbol{x})d\boldsymbol{x}\) weighted average of the GMRF values at the vertices of the triangles within the area. Weights = \(\mbox{ (number vertices)}^{-1}\)
\[S(B) \approx \frac{1}{m} \sum_{g \in B} S_g\]
Projection matrix
\[S(\boldsymbol{x}_i) \approx \sum_{g=1}^G A_{ig} S_g\ \ \ \ \ \ \ \ \ \ \ S(B_j) \approx \sum_{g=1}^G A_{jg} S_g\]
\(A\) projection matrix that maps GMRF from observations to triangulation nodes
Row \(i\) in \(A\) of point observation: three non-zero values at columns that represent vertices of triangle containing point (= barycentric coordinates)
Row \(j\) in \(A\) of areal obs: non-zero values in m vertices inside area (= \(1/m\))
\[A = \begin{bmatrix}
A_{11} & A_{12} & A_{13} & \dots & A_{1G} \\
A_{21} & A_{22} & A_{23} & \dots & A_{2G} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
A_{n1} & A_{n2} & A_{n3} & \dots & A_{nG}
\end{bmatrix} = \begin{bmatrix}
0 & 0 & 1 & \dots & 0 \\
.2 & .2 & 0 & \dots & .6 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1/m & 1/m & 1/m & \dots & 0
\end{bmatrix}\]
Multivariate downscaling of air pollutants
\[\begin{split}
\text{PM}_{2.5}(\mathbf{s}) & = \alpha_1 + z_{1}(\mathbf{s}) + e_1(\mathbf{s}) \\
\text{PM}_{10}(\mathbf{s}) & = \alpha_2 + \lambda_1 z_{1}(\mathbf{s})+z_2(\mathbf{s})+ e_2(\mathbf{s})\\
\text{Ozone}(\mathbf{s}) & = \alpha_3 + \lambda_2 z_{1}(\mathbf{s})+\lambda_3 z_2(\mathbf{s})+z_3(\mathbf{s}) +e_3(\mathbf{s})
\end{split}\]
Improve resolution from 10 to 2 km, reveal dependencies among pollutants
Rodriguez, Chacon and Moraga, The American Statistician, 2025
Spatio-temporal downscaling model
Spatio-temporal Gaussian field \(S = \{S(\boldsymbol{x}, t): \boldsymbol{x} \in D \subset \mathbb{R}^2, t \in T \subset \mathbb{R}^{+} \}\)
\[Y(\boldsymbol{x}, t) | S(\boldsymbol{x}, t) \sim N(\mu(\boldsymbol{x}, t) + S(\boldsymbol{x}, t), \tau^2)\]
- Observations at locations \(\boldsymbol{x_i}\) and times \(t_k\)
\[E[Y(\boldsymbol{x_i}, t_k)] = \mu(\boldsymbol{x_i}, t_k) + S(\boldsymbol{x_i}, t_k)\]
- Observations at areas \(B_j \subset D\) and periods of time \(\tau_l \in T\), averages of the process in space and time
\[E[Y(B_j, \tau_l)] = |B_j|^{-1} |\tau_l|^{-1} \int_{B_j} \int_{\tau_l}(\mu(\boldsymbol{x}, t) + S(\boldsymbol{x}, t)) d \boldsymbol{x} dt,\] where \(|B_j|>0\) and \(|\tau_l|>0\)
Precision disease mapping
Disease mapping is important to understand geographic and temporal patterns of diseases and allocate resources where most needed
Often, maps given at an areal resolution which difficulties decision-making
Map shows malaria prevalence in Mozambique. However, disease risk varies continuously in space & areal data unable to show how risk varies within areas
Areal estimates make difficult targeting health interventions and directing resources where most needed
Disaggregate area-level data
High-resolution estimates permit to find differences in disease risk within study regions, and identify areas and groups of people at higher risk
Bayesian spatial disaggregation model
Model assumes there is a spatially continuous variable underlying all observations that can be modeled using a zero-mean Gaussian random field
\[\begin{equation*}
\begin{aligned}
Y(\mathbf{x}) & \sim \pi \left( \theta(\mathbf{x}), \tau \right), \quad \mathbf{x} \in A \subset \mathbb{R}^2, \\
\theta(\mathbf{x}_i) & = g^{-1}\left(\mu(\mathbf{x}_i)+S\left(\mathbf{x}_i\right) \right), \quad i=1, \ldots, n, \\
\theta(B_i) & =\left|B_i\right|^{-1} \int_{B_i} g^{-1}(\mu (\mathbf{x}) + S(\mathbf{x})) d \mathbf{x}, \quad i=n+1, \ldots, n+m.
\end{aligned}
\end{equation*}\]
Zhong and Moraga, JABES, 2023
Zhong et al., JRSSA, 2024
Malaria prevalence in Madagascar
Predicting malaria prevalence by combining point and areal data
\[\begin{aligned}
Y(\mathbf{x}) & \sim \operatorname{Binomial}\left(N(\mathbf{x}), P(\mathbf{x})\right), \quad \mathbf{x} \in A \subset \mathbb{R}^2, \\
P(\mathbf{x_i}) & = \text{logit}^{-1}\left(\mu(\mathbf{x}_i)+S(\mathbf{x}_i)\right), \quad i=1, \ldots, n, \\
P(B_i) & = \left|B_i\right|^{-1} \int_{B_i} \text{logit}^{-1} \left(\mu (\mathbf{x}) + S(\mathbf{x}) \right) d \mathbf{x}, \quad i=n+1, \ldots, n+m.
\end{aligned}\]
Alahmadi and Moraga, SERRA, 2025
Bayesian spatial disaggregation model for the detection of disease clusters
A disease cluster is an unusual aggregation of cases occurring together in space. When cases are aggregated in areas, traditional methods for the detection of clusters such as scan statistics identify clusters formed by multiple areas, despite disease risk varying continuously in space
Exceedance probabilities from Bayesian spatial disaggregation models
We propose a method to detect clusters of any shape indep. of boundaries
First obtain risk surfaces from a Bayesian spatial disaggregation model
\[\begin{equation*}
\begin{aligned}
Y(\mathbf{x}) & \sim \pi \left( \theta(\mathbf{x}), \tau \right), \quad \mathbf{x} \in A \subset \mathbb{R}^2, \\
\theta(\mathbf{x}_i) & = g^{-1}\left(\mu(\mathbf{x}_i)+S\left(\mathbf{x}_i\right) \right), \quad i=1, \ldots, n, \\
\theta(B_i) & =\left|B_i\right|^{-1} \int_{B_i} g^{-1}(\mu (\mathbf{x}) + S(\mathbf{x})) d \mathbf{x}, \quad i=n+1, \ldots, n+m.
\end{aligned}
\end{equation*}\]
Then use exceedance probabilities to identify high-risk locations
\[P(\theta(\mathbf{x}) > \mbox{threshold})\]
Moraga and Alahmadi, Spatial Statistics, 2026
Detection of clusters
Through simulation, the disaggregation model showed high sensitivity and competitive specificity in a range of scenarios when compared to the circular and flexible scan statistic, and exceedance prob. from a Bayesian areal model
Real application detecting clusters of lung cancer in Pennsylvania. Clusters detected have shapes that are different from a set of areas
