Characterizing spatial correlation#
Stationarity and the semivariogram#
In geostatistics, a spatial process or random field is a collection of random variables [Cressie, 1993]:
where \(D \subset \mathbb{R}^d\) is the spatial domain, \(d \in \mathbb{N}\) is the dimension of the space, and \(Z(\mathbf{x})\) is a random variable that represents the value of some spatially varying quantity (e.g. temperature) at location \(\mathbf{x}\). Importantly, the term process refers to this entire stochastic system, not a single observation at any specific location [Cressie, 1993, Stein, 1999]. A spatial process is said to be second-order stationary if its expected value is constant and its covariance depends only on the separation between two locations:
where \(m\) is a constant expected process mean and \(\mathbb{C}(\cdot)\) the covariance function [Cressie, 1993]. Here \(\mathbf{h}\) denotes the lag vector (direction and magnitude) between two locations and \(h\) denotes the Euclidean length (scalar lag distance). Either \(\mathbf{h}\) or \(h\) can be used depending on whether directionality matters, which is discussed later (isotropy and anisotropy). For a specific pair of separate locations indexed by \(i\) and \(j\) (\(i, j \in \{1, \dots, n_s \}\), where \(n_s\) is the number of samples) define \(\mathbf{h}\) and \(h\) by:
where \(x_{i,k}\) denotes the \(k\)th coordinate of \(\mathbf{x}_i\) and \(k = 1 \dots d\) indexes spatial coordinates (e.g., \(d=1\) for a line, \(d=2\) for a plane and \(d=3\) for three‑dimensional space). The specific pairwise scalar distance used in empirical sums and binning is written \(h_{ij}=\|\mathbf{x}_j-\mathbf{x}_i\|\), whereas plain \(h\) denotes a generic lag argument in theoretical expressions such as \(\gamma(h)\) or \(\mathbb{C}(h)\).
Theoretical vs. empirical expressions
While \(Z(\mathbf{x})\) denotes the random field, the quantity \(z(\mathbf{x})\) can be used to express a realized observation at location \(\mathbf{x}\). This distinction makes clear whether an expression is a theoretical statement or an empirical computation.
This information helps characterize spatial dependence by using the semivariogram, which describes how the average squared difference between values changes with distance. The semivariance at lag \(h\) is defined as:
where \(\mathbb{C}(0)\) is the variance at zero lag (the process variance). In this context, the semivariogram and covariance function contain equivalent information about spatial correlation [Chilès, 2012, Webster and Oliver, 2007]. In this document, the scalar \(h\) denotes a generic lag distance.
Isotropy and anisotropy#
A spatial process is said to be isotropic if its spatial dependence structure is the same in all directions; that is, the covariance and semivariogram functions depend only on the distance between locations, not on the direction. Mathematically, this means that the semivariogram and covariance are functions of \(h\), hence Eq. (3.14). This assumption simplifies variogram analysis, as all pairs of points separated by the same distance are treated equivalently, regardless of their orientation.
In contrast, a process is anisotropic if its spatial dependence varies with direction. In this case, the covariance and semivariogram depend on both the magnitude and orientation of \(\mathbf{h}\):
Directional (azimuthal) variogram analysis is required to detect and model anisotropy. This involves grouping pairs of locations not only by lag distance but also by their relative orientation (azimuth angle), allowing the spatial correlation structure to be characterized as a function of both distance and direction.
The intrinsic model#
However, real-world datasets often do not fulfill the requirements for full second-order stationarity. The intrinsic hypothesis provides a weaker condition by requiring only that the increments of \(Z(\mathbf{x})\), \(Z(\mathbf{x}) - Z(\mathbf{x}+\mathbf{h})\), be second-order stationary. Under this weaker condition, the semivariogram is now defined as:
which means that, for a given \(\mathbf{h}\), all possible pairs of locations in the spatial domain separated by \(\mathbf{h}\) are considered [Chilès, 2012, Montero, 2015]. The expectation \(\mathbb{E}[\cdot]\) then represents the average of these squared differences over all such pairs. Thus, the intrinsic model generalizes the covariance-based approach, allowing spatial correlation to be characterized through the semivariance even when the covariance function is not well-defined. In practice, however, \(Z(\mathbf{x})\) is not known everywhere; only a finite set of observations at sampled locations is available.
To estimate the semivariogram from data, the theoretical expectation is replaced with an average over all observed pairs of points separated by lag \(h\). This leads to the empirical semivariogram estimator:
where \(z(\mathbf{x}_i)\) and \(z(\mathbf{x}_j)\) denote the observed (realized) values of the variable at sample locations \(\mathbf{x}_i\) and \(\mathbf{x}_j\), respectively. Each \(\mathbf{x}_i\) (or \(\mathbf{x}_j\)) is a location in the sampling set; \(z(\mathbf{x}_i)\) is the concrete measurement recorded at that location (scalar unless a multivariate observation is specified). This estimator provides a practical method for computing the semivariogram from data, enabling empirical characterization of spatial dependence.
Practical notes on empirical computation#
When computing the empirical semivariogram, the choice of lag bins is important. Lag bins are typically defined by dividing the range of observed pairwise distances into intervals of fixed width, or by using quantiles to ensure a similar number of pairs per bin. The bin width and maximum lag should be chosen to balance resolution and statistical reliability; bins with too few pairs can yield unstable estimates. In practice, a tolerance is used so that all pairs with \(h_{ij}\) within a specified interval around \(h\) are included in the bin for \(h\). Uneven sampling or missing data may result in some bins having very few pairs, which can be mitigated by adjusting bin widths or excluding sparse bins from model fitting. Formally, \(N(h)\) is the count of unique pairs falling in the bin for \(h\):
where \(\Delta h\) is the bin width and \(\mathbf{1}\{\cdot\}\) is the indicator function.
Standardizing empirical variograms#
The classic empirical semivariogram assumes homoscedasticity (homogeneous variance) across space. Many real datasets exhibit heteroscedasticity, where local variance changes with location due to patchiness, measurement error, or sampling design. This can bias or destabilize semivariogram estimates. To reduce sensitivity to local variance differences and outliers, the empirical semivariogram from Eq. (3.7) is normalized by the product of their local standard deviations:
where \(\sigma(\mathbf{x}_i)\) and \(\sigma(\mathbf{x}_j)\) are local standard deviations at \(\mathbf{x}_i\) and \(\mathbf{x}_j\). This standardization improves robustness to local variance differences and outliers, providing more reliable estimates of spatial correlation in heterogeneous datasets [Goovaerts, 1997, Webster and Oliver, 2007].
Semivariogram features and interpretation#
Theoretical semivariogram models are parametric functions designed to capture the essential features of spatial correlation observed in empirical semivariograms. These models are characterized by several key parameters:
Nugget (\(C_0\)): The value of the semivariogram at zero lag. The nugget represents measurement error or micro-scale variation that occurs at distances smaller than the sampling interval. A large nugget indicates substantial noise or unresolved spatial structure.
Sill (\(C\)): The plateau value that the semivariogram approaches at large lag distances. The sill reflects the total variance of the process. In bounded models, the semivariogram levels off at the sill; in unbounded models, it may increase indefinitely.
Range (\(a\)): The lag distance at which the semivariogram reaches the sill (or a specified fraction of the sill, such as 95%). The range indicates the scale of spatial correlation: pairs of points separated by less than the range are correlated, while those farther apart are effectively uncorrelated.
Model shape: The functional form of the model determines how quickly spatial correlation decays with distance. For example, some models rise linearly near the origin and flatten at the sill (spherical), others rise rapidly and approach the sill asymptotically (exponential), and some are very smooth near the origin (Gaussian). The choice of model shape reflects the expected smoothness and continuity of the underlying spatial process.
Anisotropy: In some cases, spatial correlation may vary with direction. Anisotropic models allow the range or sill to depend on the orientation of the lag vector, capturing directional effects in the data.
Interpreting these features provides insight into the spatial structure of the data, the scale of spatial dependence, and the presence of measurement error or micro-scale variability [Chilès, 2012, Cressie, 1993].
A comprehensive list of theoretical variogram models, including their mathematical equations and descriptive properties, is provided in the semivariogram equations reference document. Readers should consult that document for details on specific models, their parameterizations, and guidance on model selection for different types of spatial processes.
Diagnostics and interpretation#
After fitting a theoretical variogram model, it is important to assess the quality of the fit. Visual inspection of the empirical and fitted curves can reveal systematic deviations, such as underfitting at short or long lags, or failure to capture anisotropy. Quantitative diagnostics, such as cross-validation or analysis of kriging residuals, can help evaluate whether the chosen model adequately represents spatial dependence. Poor fits may indicate the need for alternative model shapes, inclusion of anisotropy, or further data cleaning. Interpreting the fitted parameters in the context of the data and sampling design is essential for drawing reliable conclusions about spatial structure.
Semivariogram model fitting#
After estimating the empirical semivariogram from spatial data, a theoretical variogram model is fit to the empirical values to obtain a smooth, interpretable representation of spatial correlation. Theoretical models are parametric functions of the form:
where \(C_0\) is the nugget (the value at zero lag), \(C_1\) is the partial sill (\(C - C_0\)), and \(\mathscr{f}(h; \theta)\) is a correlation function that depends on lag \(h\) and model-specific parameters \(\theta\) (e.g., range, smoothness, periodicity) [Matérn, 1986, Stein, 1999, Whittle, 1954]. The process of fitting a variogram model involves finding the set of parameters \(\theta\) that best match the empirical semivariogram. This is typically achieved by minimizing the weighted sum of squared differences between the empirical semivariogram values \(\gamma_\text{empirical}(h_i)\) and the theoretical model \(\gamma_\text{model}(h_i; \theta)\) at each lag \(h_i\):
where \(n\) is the number of lag bins and \(b\) indexes those lag bins (\(b=1, \dots, n\)) [Cressie, 1993, Webster and Oliver, 2007]. The weights \(w_b\) are associated with each lag bin. Values of \(w_b\) are often chosen to be proportional to the number of data pairs in the \(b\)‑th bin, \(N(h_b)\), so that lag bins with more data have greater influence on the fit. When no standardization is applied, these weights are defined as:
In contrast, when \(w_b\) are normalized such that they sum to 1:
Assumptions and limitations#
Variogram analysis assumes at least intrinsic stationarity of increments and, for second‑order stationarity, a constant mean [Chilès, 2012, Cressie, 1993]. If the mean exhibits spatial trend, detrend or model the mean prior to variogram estimation [Goovaerts, 1997]. Variogram estimation is sensitive to outliers, edge effects, and irregular sampling; handle these issues explicitly and document the choices made [Webster and Oliver, 2007]. Standardized estimators, robust estimators (e.g. Cressie–Hawkins), and careful binning can mitigate some practical problems, but they do not replace the need to evaluate assumptions and diagnostics [Cressie, 1993, Montero, 2015].