Spatial Statistics

"Everything is related to everything else, but near things are more related than distant things." — Tobler's First Law of Geography

Motivations for Our Spatial Statistics Research:

To study real phenomena, one might ideally take measurements at infinitely many spatial coordinates. This indexed collection of measurements is what we call a spatial field. In practice, however, dense measurements are often infeasible or impossible to obtain (for example, subsurface topography). Instead, we model the field as a random process. An ideal random field captures the essential structure of a complicated phenomenon using only a few physically meaningful parameters. We assume the field is Gaussian because Gaussian fields are fully characterized by their first- and second-order statistics. Of particular importance are the second-order statistics, since they describe correlation structure in space.

Spatial correlations can be interpreted in terms of oscillatory structure. This motivates Fourier analysis and the Wiener–Khinchin relationship between covariance and spectrum. We additionally assume isotropy as the null hypothesis. We use a zero mean and a Matérn covariance model. We choose Matérn covariance because its parameters —variance, smoothness, and range— are physically meaningful.

Although the spatial and spectral domains contain equivalent information, the spectral domain is often computationally advantageous because it approximately diagonalizes the problem. For an infinite stationary field, different wavenumbers are uncorrelated in the spectral domain. The major complication is that real observations are finite and sampled. Finite spatial windows introduce spectral leakage/blurring (“boundary effects”) and induce correlations between nearby wavenumbers. Increasing the spatial domain helps reduce these effects because larger windows better capture long-range correlations and reduce spectral leakage.

To make likelihood-based inference computationally feasible, we often use the Whittle likelihood, which approximates the full covariance likelihood in the spectral domain. To assess parameter resolution and uncertainty, we study quantities such as the score function, Hessian, Fisher information matrix, and residual diagnostics (which will follow a chi-squared with 2 degrees of freedom due to the quadratic form of the likelihood). The Fisher information matrix measures how sensitive the likelihood is to changes in model parameters. Larger values generally indicate parameters that are better resolved by the data.



Key Spatial Questions:

  • It gives us information about spatial relationships
  • It allows us to model random fields (stochastic processes in multiple dimensions)
  • It can be applied to a number of problems (interpolation, extrapolation, simulation, modeling, growing-domain, etc.)

Spatial statistics is the study of an indexed family of random variables, $Z(\textbf{s})$. A random variable can be thought of as a mapping from outcomes to numbers. The Z notation for a random variable is just a nod to the fact that we have the ability to work with complex values. The "index" refers to the spatial vector s that can live in any dimension. For instance, a cartesian map with 5 points on it will be represented as $(\textbf{s}_1,\textbf{s}_2,\textbf{s}_3,\textbf{s}_4,\textbf{s}_5)$ where each $\textbf{s}_i=(\phi_i,\lambda_i)$ for latitude and longitude. Lastly, something to be aware of is that functions of random variables are random variables themselves.

A Gaussian process is a generalization of a random variable. More specifically, it is a distribution over functions. Instead of each "draw" or "realization" being a number, each realization is an entire function called a sample path. Four 1D sample paths in the time/spatial domain are pictured below, each generated with the same mean and covariance structure.

Gaussian Process Sample Paths
Covariance Structure

The covariance kernel/matrix/function provides the underlying structure of the field (tunable through parameters such as correlation length or variance). Choosing a covariance structure can be a bit arbitrary at times. In general, you want to model it off of the properties you observe in your time series/data (is it periodic? stationary? smooth?...). I refer the interested reader to the Kernel Cookbook . The only rules for choosing a kernel is that it has to be symmetric and positive semi-definite. To incorporate randomness, you are essentially transforming white noise with the covariance kernel:

$z \sim N(0,1)$

$C(t,t') = LL^T$

${cov}(Lz)=L^Tcov(z)L=LIL^T=LL^T=C(t,t')$

One can imagine that if you sample more and more, a distribution over the sample paths will emerge. That is called your latent function distribution, which I like to think about as the background function. We choose a Gaussian distribution because it is uniquely determined by only two measurements (mean and covariance), is closed under conditioning, and is closed under marginalization.

Five Key Facts
  1. A Gaussian process is completely determined by its first and second moments
  2. All finite subsets of Gaussian processes must be jointly Gaussian
  3. Covariance controls structure
  4. More data reduces uncertainty
  5. Conditioning updates beliefs

Convergence analysis allows us to make statements about long-term asymptotics of a sequence of random variables $X_n$. This is amazing because it can tell us if random events settle into patterns (whether about the mean, probability, or distribution).

In probability theory, there exist three types of convergence of sequences of random variables (ordered here from weakest to strongest):

  • Convergence in distribution (law): $F_{X_n}(x)\rightarrow F_X(x)$ (think: the histograms match)
    • If you roll a fair die many, many times, you'd expect the histogram to converge to the uniform distribution.
  • Convergence in probability: $P(|X_n-X|>\epsilon)\rightarrow0$ (think: they are close)
    • If you were to plot the height of more and more people, the true mean human height would emerge.
  • Almost sure convergence: $P(X_n\rightarrow X)=1$ (think: we know how the sample paths will converge)
    • If you studied how much a specific animal eats per day, that will fluctuate. But one day, you will be almost sure it will be 0.

The CLT and WLLN both study the same object, namely the sample mean $\bar X_n=\frac1N\sum X_i$, but ask two different questions about it:

  • CLT: What is the long-term distribution of the sample mean? $ \frac{\sqrt n}{\sigma} (\bar X_n-\mu) \rightarrow N(0,1) $
  • WLLN: What does the sample mean get closer and closer to? $ \bar X_n \rightarrow \mu $

First, it is important to motivate this question. Many estimation methods, such as the Whittle likelihood, rely on the assumption of stationarity to simplify the covariance structure in the spectral domain. Stationarity implies that the statistical properties of the process do not change with shifts in space. As a result, the covariance function depends only on the distance between two points, not on their absolute locations. This leads to a spectral density that is a function of wavenumber alone, without any cross-terms that would indicate correlation between different wavenumbers. In other words, under stationarity, the spectral density is diagonal in the wavenumber domain, which means that different wavenumbers are uncorrelated.

$\langle dH(k) dH(k`) \rangle = S(k) \delta(k-k`) dk dk`$

Breaking down the Matérn equation, we will be able to see that it is always non-negative. The Matérn covariance function is given by:

$C_{\mathbf{\theta}}(\mathbf r) &= \sigma^2 \frac{2^{1-\nu}}{\Gamma(\nu)} \left( \frac{2\sqrt{\nu}} {\pi\rho} \mathbf{r} \right)^\nu K_\nu \left( \frac{2\sqrt{\nu}} {\pi\rho} \mathbf{r} \right)$

We first note that all the parameters are non-negative ($\rho >0$, $\nu >0$, $\sigma^2 >0$). A potential source of negativity is the modified Bessel function of the second kind, $K_\nu$. However, $K_\nu(z)$ is positive for all $z > 0$ and $\nu > -1$. Therefore, the entire expression for the Matérn covariance function is non-negative for all valid parameter values and spatial lags $\mathbf{r}$ (which are always positive). Finally, the Gamma function $\Gamma(\nu)$ is given by:

$\Gamma(\nu) = \int_0^\infty t^{\nu-1} e^{-t} dt$

Which you can verify is always positive. Interestingly, negative processes do exist, but they are not commonly used in spatial statistics.

To figure out if your data is isotropic or an-isotropic, you can ask "is there a preferred direction in the spatial correlation structure or not?" A way you can test for this is by analyzing the residuals in wavenumber space. Ideally, the residuals $\chi_{\theta}(\mathbf{k}) = \hat{S}^{-1}(\mathbf{k}) |H(\mathbf{k})|^2$ should follow a $\chi^2$ statistic with 2 degrees of freedom ($\sim \frac{\chi_2^2}{2}$). You can this visually by plotting $\chi^2(\mathbf{k})$ and looking for structure.

Another way is to look at a quantile-quantile plot. If the points deviate from the line that correspons to ranks of a \chi_{\theta}(\mathbf{k}), that is a sign of an-isotropy. Below are example of isotropic and an-isotropic datasets illustrating these ideas. (These are taken from: Simons, Frederik J., Olivia L. Walbert, Arthur P. Guillaumin, Gabriel L. Eggers, Kevin W. Lewis, and Sofia C. Olhede. 2026. “Maximum-Likelihood Estimation of the Matérn Covariance Structure of Isotropic Spatial Random Fields on Finite, Sampled Grids.” Geophysical Journal International 245(2):ggag044. doi:10.1093/gji/ggag044)
Isotropic Data Anisotropic Data

A question that is good to always have in the back of one's mind pertaining to our research is "is isotropy a good enough model of the lithosphere?"

The two ways we focus on are:

  • Cholesky decomposition: This method was used in the Simons & Olhede 2013 paper (https://doi.org/10.1093/gji/ggt056) and uses the sgp.m MATLAB code (https://github.com/csdms-contrib/slepian_juliet/blob/master/sgp.m). It operates on the blurred spectrum but does not account for wavevector correaltions due to boundary effects.

  • Circulant embedding: This method involves embedding the covariance matrix in a larger circulant matrix in the spatial domain (mimicking periodicity) and then using the Fast Fourier Transform to generate realizations. It does account for exact blurring. It was used in the Simons et. al. 2026 paper (https://doi.org/10.1093/gji/ggag044) and uses the MATLAB code (https://github.com/csdms-contrib/slepian_juliet/blob/master/randgpn.m).

    • The Whittle likelihood is an approximation to the full covariance likelihood that is computationally more tractable, especially for large datasets. It is based on the idea of transforming the data into the spectral domain, where the covariance structure becomes simpler under the assumption of stationarity. The Whittle likelihood for a parameterized Matérn covariance model can be expressed as:

    $L(\theta) = -\sum_{\mathbf{k}} \left[ \log S_{\theta}(\mathbf{k}) + \frac{|H(\mathbf{k})|^2}{S_{\theta}(\mathbf{k})} \right]$

    Where $S_{\theta}(\mathbf{k})$ is the spectral density corresponding to the Matérn covariance function parameterized by $\theta$, and $H(\mathbf{k})$ is the Fourier transform of the observed data. The first term in the summation accounts for the model complexity, while the second term measures the fit of the model to the data. By maximizing this likelihood with respect to the parameters $\theta$, we can estimate the parameters of the Matérn covariance function that best explain the observed spatial data.
    We choose to use the debiased Whittle likelihood (a quasi-likelihood), which corrects for bias introduced by finite sample sizes and spectral blurring, because it delivers asympototically unbiased, gaussian estimates of Matérn parameters under various sampling senarios (Simons et. al. 2026). The following venn diagram compares the Whittle likelihood to debiased Whittle likelihood: Likelihood Comparison The debiased Whittle likelihood is given by:
    $\bar{L}(\theta) = -\sum_{\mathbf{k}} \left[ \log \bar{S}_{\theta}(\mathbf{k}) + \frac{|H(\mathbf{k})|^2 }{\bar{S}_{\theta}(\mathbf{k})} \right]$

    There are three approaches to asymptotic analysis in spatial covariance estimation that we focus on:

  • Growing Domain: This approach considers the behavior of estimators as the spatial domain increases indefinitely (holding $\Delta x$ and $\Delta y$ fixed). It is particularly relevant for understanding how well our methods perform when we have access to more and more spatial data.
  • Infill (fixed-domain): This approach studies the behavior of estimators as the spatial sampling becomes denser within a fixed domain (holding the domain size ($(N_x \Delta x, N_y \Delta y)$) fixed while $\Delta x$ and $\Delta y$ approach zero). It is useful for understanding the effects of increasing spatial resolution on parameter estimation.
  • Reduced missingness: Essentially, this eperiment replaces missing samples with data.


    • The main implications for the Matern parameters specifically are:
    1) For the growing domain, as the grid grows, relative estimator bias, variance, covariance, and RMSE all decrease.
    2) For the infill, as the grid becomes denser, a reduction of parameter variance for $\nu$, bias in $\nu$ and $\rho$ that didn't subside even with more data, reduced varaince of $\sigma^2$, and mixed covariances are observed.
    3) For the reduced missingness, essentially, the more observations of $\rho$ that fit in the grid experiment, the better.
    For more information, see Walbert, O. L., Simons, F. J., Guillaumin, A. P., & Olhede, S. C. (2026). Irregularly and incompletely sampled random fields in the Earth sciences: Analysis and synthesis of parameterized covariance models (arXiv:2604.02179). arXiv. https://doi.org/10.48550/arXiv.2604.02179

    The spectral window refers to the fact that all observed data is finite and sampled, which leads to spectral blurring/leakage. The role of the spectral window in spectral covariance estimation is to account for the effects of this blurring on the observed periodogram. When we compute the periodogram from finite data, it does not directly reflect the true spectral density of the underlying process due to the convolution with the spectral window. Therefore, when estimating parameters of a covariance model, we need to consider how the spectral window modifies the observed periodogram and adjust our estimation methods accordingly (e.g., using a debiased Whittle likelihood that incorporates the expected value of the windowed periodogram).

    The primary role of the expectation of the windowed periodogram in parametrized covariance estimation is to mitigate the effects of spectral leakage, which occurs when the finite spatial window causes energy from one wavenumber to leak into others. The expectation of the windowed periogram approaches the true spectral density. The more averages of periodograms you have, the closer the expectation of the windowed periodogram will be to the true spectral density. Our research "proved" this by averaging 100 periograms and also by showing that the plot of the expectation of the windowed periodogram divided by the blurred spectrum was static. Windowed Periodogram
    The variance of the windowed periodogram, on the other hand, estimates the blurred power spectral density (the diagonals of plots b-d are rewrapped below). Windowed Periodogram Variance

    1) Acknowledge that the expectation of the spectral desnity $E[\hat{S}(k)]$ is the blurred spectrum and so you should also blur the model spectrum $S_{\theta}(k)$ when you are doing parameter estimation.
    2) Actually capture the off-diagonals thorugh circulant embedding.

    The Vecchia approximation is a method for approximating the likelihood of a Gaussian process by decomposing the joint distribution into a product of conditional distributions. It is particularly useful for large datasets where computing the full covariance matrix is computationally infeasible. The approximation works by ordering the data points and then conditioning each point on a subset of the previous points in the order. This results in a sparse precision matrix, which allows for efficient computation of the likelihood and parameter estimation. The Vecchia approximation can be seen as a generalization of the nearest neighbor Gaussian process and is widely used in spatial statistics for modeling large spatial datasets.



    back