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

7.3 Model Fitting Via The Variogram

Using (7.1) as the model of interest, our ultimate goal is to predict the process $Y$ at locations where we do not have data. In the mining industry (the origin of geostatistics), this concept has obvious utility: based on core samples extracted from a set of locations, where should we set up our mine (a costly and time-intensive operation) so it will be most productive? In environmental epidemiology, we are often concerned with estimating the burden or risk of disease: similar concepts can be applied.

Spatial prediction (or kriging) via variogram estimation (as opposed to maximum likelihood estimation) usually involves at least four stages: (i) produce initial estimates of $\beta$ using ordinary least squares (OLS) regression (ii) produce estimates of $\tau$ , $\sigma$ and $\phi$ (iii) re-estimate $\beta$ using EGLS (iv) produce kriged estimates of $Y$ . Step (iii) is optional, but if we are interested in making inference about the parameters, $\beta$ , then we do really need to take account of the fact that the observations are not independent (which is assumed by OLS) before presenting the results of our analysis. Steps (ii) and (iii) can also be iterated.

One of the main ingredients for producing spatial predictions of $Y$ is an estimate of the surface $S(x)$ and its variability at locations $x$ where we do not necessarily have data. In order to obtain estimates of $S(x)$ and $\mathbb{V}[S(x)]$ we first estimate the second order properties of $S+\varepsilon$ . The variogram is an exploratory tool for doing this; it is also used in longitudinal data analysis.

The variogram is defined as:

V(x,x^{\prime})=\frac{1}{2}\mathbb{V}\{S(x)-S(x^{\prime})\}

For a stationary process, this quantity (known as the semivariance) can be estimated as

V(u_{ij})=\frac{1}{2}\mathbb{E}[(Y_{i}-Z_{i}\hat{\beta})-(Y_{j}-Z_{j}\hat{% \beta})]^{2},

for each distance $u_{ij}=\|x_{i}-x_{j}\|$ , where $Z_{i}$ is the covariate data for the $i$ th individual and $\hat{\beta}$ has been computed using OLS (for example). The resulting point estimates i.e. $\frac{1}{2}[(Y_{i}-Z_{i}\hat{\beta})-(Y_{j}-Z_{j}\hat{\beta})]^{2}$ , computed using all possible pairs $(i,j)$ , can be plotted against $u_{ij}$ as a ‘variogram cloud’ (Figure 7.2, left plot), or averaged for similar values of $u_{ij}$ to produce a ‘binned variogram’ (Figure 7.2, right plot). In this example, we have not used covariates, so Figure 7.2 shows the variogram for $\frac{1}{2}\mathbb{E}[Y_{i}-Y_{j}]^{2}.$

Figure 7.2: Link, Caption: Estimating the variogram from the magnesium soil content data, see Example 7.3. The left plot shows a variogram cloud, which is computed from the data: for all pairs of points

(i,j)

, their physical distance

u_{ij}

(on the

x

-axis) is plotted against the empirical semi-variance,

\frac{1}{2}[(Y_{i}-Z_{i}\hat{\beta})-(Y_{j}-Z_{j}\hat{\beta})]^{2}

(on the

y

-axis) - this forms a cloud of points. In order to compute an estimate of the variogram, we can bin the distances on the

x

-axis and plot these against the mean of the empirical semivariances within each bin – this leads to a binned variogram, shown on the right. The binned variogram reveals a broadly increasing relationship, as might be expected: we are working with a hypothesis that things closer together are more similar than things further apart. The next stage in such an analysis might be to postulate a model for the relationship in the variogram, and thus estimate how quickly dependency decays with respect to distance.

The variogram helps to provide plausible initial estimates of the parameters of the process, as it can be shown that

V(u)=\tau^{2}+\sigma^{2}\{1-\rho(u)\}.

Fitting the stationary Gaussian model with exponential correlation function using ordinary least squares regression, we obtain:

\hat{\sigma}^{2}=35.2,\;\hat{\phi}=126.267,\;\hat{\tau}^{2}=8.303

$\sigma^{2}$ is known as the variance parameter, $\phi$ as the range parameter, $\tau^{2}$ as the measurement error parameter or nugget effect, and $\sigma^{2}+\tau^{2}$ as the sill.

Figure LABEL:mgvarianceplot.pdf below shows the shape of the fitted covariance function, i.e. $\sigma^{2}\rho(u)$ with parameters set at the ordinary least squares estimates.

Figure 7.3: Link, Caption: Showing the fitted covariance function (covariance on the

y

-axis plotted against distance,

u

on the

x

-axis) from the analysis of the magnesium soil content data, Example 7.3. As distance increases, the covariance decreases. Note that for interpretability, it is usually better to visualise the correlation function, rather than the covariance function, as it is scale-independent.

Example 7.6.

Exam Question 2014

Suppose that it is desired to fit the following model to a set of geostatistical data

	$\displaystyle Y_{i}$	$\displaystyle=$	$\displaystyle X_{i}\beta+S(x_{i})+Z_{i},$
	$\displaystyle Z_{i}$	$\displaystyle\sim$	$\displaystyle\mathrm{N}(0,\tau^{2}),$

where $\{Y_{1},\ldots,Y_{n}\}$ are measurements at locations $\{x_{1},\ldots,x_{n}\}$ ; $X_{i}$ is a vector of covariate information at each $x_{i}$ ; and $S(x_{i})$ is the value of a zero-mean second order stationary spatial Gaussian process at $x_{i}$ . Suppose $S$ has an exponential correlation function with variance parameter $\sigma^{2}$ and spatial decay parameter $\phi$ . Figure 7.4 shows the variogram used to produce the estimates $\sigma^{2}$ , $\phi$ and $\tau^{2}$ used to create $\Sigma$ . Use this plot to suggest estimates for $\sigma^{2}$ , $\phi$ and $\tau^{2}$ stating any results you use.

Figure 7.4: Link, Caption: Plot of fitted variogram for Example 7.6. This plot illustrates an exponential model for the variogram (blue line) with points from the raw binned variogram as blue dots; the numbers are not relevant for this example.

Solution:

The variogram helps to provide plausible initial estimates of the parameters of the process $S$ as it can be shown that $V(u)=\tau^{2}+\sigma^{2}\{1-\rho(u)\}$ , where $\rho$ is the correlation function. Since $\rho(u)$ tends to zero as $u\rightarrow\infty$ , the sill is an estimate of $\tau^{2}+\sigma^{2}$ and since $\rho(0)=1$ , the intercept is an estimate of $\tau^{2}$ , hence from the plot $\sigma^{2}=68118-11915=56203$ , $\tau^{2}=11915$ (the intercept). For $\phi$ choose a point on the line and note the value of $u$ and $V(u)$ , rearranging the above, we get $\phi=-u/\{\log[1-(V(u)-\tau^{2})/\sigma^{2}]\}$ with $u=500$ and $V(u)=56000$ , this gives $\phi=326$ (the exact answer is 293).