7 Continuous Spatial Variation 7 Continuous Spatial Variation 7.2 Model Fitting Via Maximum Likelihood

7.1 Geostatistical models

Definition 7.1.

Let $S(x)$ , $x\in\mathbb{R}^{2}$ be a real-valued random variable. We say that $S$ is a spatial Gaussian process if for any finite collection of points on the plane, $x_{1},\ldots,x_{n}\in\mathbb{R}^{2}$ , the joint distribution of the random vector $[S(x_{1}),\ldots,S(x_{n})]$ is multivariate Gaussian.

Definition 7.2.

Let $S$ be a Gaussian process. If in addition,

1.

$\mathbb{E}[S(x)]=m$ for all $x\in\mathbb{R}^{2}$ and
2.

$\mathrm{corr}[S(x),S(x+u)]=\rho(u)$

then $S$ is called a weakly stationary or a second-order stationary Gaussian process. In this chapter, we will refer to $S$ as simply a stationary Gaussian process and we will write,

S(x)\sim{\rm SGP}\{m,\sigma^{2},\rho(u)\},

where $\sigma^{2}=\mathbb{V}[S(x)]$ . $\rho(u)$ is called the correlation function.

Example 7.1.

An example correlation function is the exponential correlation function

\rho(u;\phi)=\exp(-u/\phi).

The parameter $\phi$ controls how quickly the spatial dependence between points at a distance $d$ apart drops off as a function of the distance apart: in a large value of $\phi$ represents long-range correlations and a very small $\phi$ means that points close together are nearly independent.

A typical model for measurements $\{Y_{1},\ldots,Y_{n}\}$ at locations $\{x_{1},\ldots,x_{n}\}$ , given covariate information, $\{z_{1}(x_{i}),\ldots,z_{p}(x_{i})\}$ at each $x_{i}$ , is

$\displaystyle Y_{i}$	$\displaystyle=$	$\displaystyle\mu(x_{i})+S(x_{i})+\varepsilon_{i},$	(7.1)
$\displaystyle\mu(x_{i})$	$\displaystyle=$	$\displaystyle\sum_{k=1}^{p}z_{k}(x_{i})\beta_{k},$	(7.2)
$\displaystyle\varepsilon_{i}$	$\displaystyle\sim$	$\displaystyle\mathrm{N}(0,\tau^{2}),$

where $m=\mathbb{E}[S(x_{i})]=0$ . Models of this kind are often called geostatistical models (Cressie, 1991). This is an indirect reference to their historical development in connection with spatial prediction problems in the mining industry.

Example 7.2.

Typical geostatistical problem: use data $Y_{i}$ from locations $x_{i}$ to predict either $S(x)$ itself, or a functional such as the average value of $S$ over region $A$ ,

T=\frac{1}{|A|}\int_{A}S(x)\mathrm{d}x,

where $|A|$ is the area of region $A$ . Note that prediction locations, $x$ , are often on a regular grid over the observation window of interest: this allows us to produce raster images of the prediction surface.

Example 7.3.

The dataset camg in the geoR package (having loaded the package, simply type data(camg) at the console) is an example of a geostatistical dataset containing, among other things, the concentration of Magnesium in a set of soil samples taken over a geographical region. Figure 7.1 shows the magnesium soil content (0-20cm below surface) measured at various locations in a field.

Figure 7.1: Link, Caption: Magnesium soil content 0-20cm below surface. The polygon denotes the study area and the circles are locations where soil samples were taken; the size of each circle is proportional to the concentration of magnesium observed in each soil sample. In this example, the locations of the samples (incidentally) form a fairly regular grid over the observation window

Figure 7.1 is an example of a typical geostatistical dataset: we observe a quantity of interest (and possibly some covariates) at a set of locations. A key property of geostatistical data is that the locations of the data are fixed by design and are not informative about the underlying spatial process; although it is possible to perform inference for the case where the locations themselves are modelled by a stochastic process (Diggle et al., 2010; Taylor et al., 2015, 2018, 2019, 2020).

How do we go about fitting a geostatistical model? There are three main methods:

Method 1: Use maximum likelihood to obtain estimates of $\beta$ , $\sigma$ , $\phi$ and $\tau^{2}$ .
Method 2: Use the variogram (see below) to obtain estimates of $\beta$ , $\sigma$ , $\phi$ and $\tau^{2}$ .
Method 3: Use Bayesian methods to obtain estimates of $\beta$ , $\sigma$ , $\phi$ and $\tau^{2}$ : write down priors for these parameters and produce inferential statements from the posterior, $\pi(\beta,\sigma,\phi,\tau|\text{data})$ . We do not cover Bayesian estimation for Geostatistical models here.

Having obtained estimates of the parameters, we can predict the process $Y$ over the spatial region of interest (i.e. in places where we do not have data). Note that in the case that our model includes covariates $\{z_{k}(x_{i})\}_{k=1}^{p}$ at each spatial location $x_{i}$ , we would need to know the values of the covariates at each of the prediction locations in order to be able to predict $Y$ (though predicting $S$ is still possible without covariates). The process of forming predictions of $Y$ (or $S$ ) is known as kriging.