\[ \newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]} \newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]} \newcommand{\SampleVar}[1]{\widehat{\mathrm{Var}}\left[ #1 \right]} \newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]} \newcommand{\TrueRegFunc}{\mu} \newcommand{\EstRegFunc}{\widehat{\TrueRegFunc}} \DeclareMathOperator{\tr}{tr} \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator{\det}{det} \newcommand{\TrueNoise}{\epsilon} \newcommand{\EstNoise}{\widehat{\TrueNoise}} \newcommand{\AutoCov}{\gamma} \newcommand{\LinearExpect}[1]{\mathbb{L}\left[ #1 \right]} \]

In our previous episode

• General approach to optimal linear prediction
  • Predict \(Y\) from \(\vec{Z} = [Z_1, Z_2, \ldots Z_p]\)
  • Best linear prediction is \(\LinearExpect{Y|Z} = \alpha + \vec{\beta} \cdot \vec{Z}\)
  • Best \(\alpha = \Expect{Y} - \vec{\beta} \cdot \Expect{\vec{Z}}\)
  • Best \(\beta = \Var{\vec{Z}}^{-1} \Cov{\vec{Z}, Y}\)
• Today: application to spatial and spatio-temporal data
  • What can this do for us?
  • How do we find the covariances?

This should look familiar

• It’s exactly the same as optimal linear prediction for time series
• “Using linear prediction on spatial data” is called kriging
  • We’ll come back to the name later

What do we need to know?

• \(\Expect{X(r)}\)
• \(\Var{X(r)}\)
• \(\Cov{X(r), X(q)}\)
• Problem: we only have one observation of each \(X(r)\)

Expectation values

• We can estimate \(\Expect{X(r)}\) by smoothing
  • Assuming \(\TrueRegFunc(r)\) doesn’t change very rapidly
• Need to pick a smoother
  • Assuming \(\Expect{X(r)} =\) constant is one choice of smoother…
  • Moving averages can work, splines are nice…
• Or we can assume spatial stationarity and use the global mean
  • See back-up about (potentially) improving on this

What about variance and covariance?

\[ \Cov{X(r), X(q)} = \Expect{X(r)X(q)} -\TrueRegFunc(r) \TrueRegFunc(q) = \Expect{(X(r)-\TrueRegFunc(r))(X(q) - \TrueRegFunc(q))} \]

• If we saw \(m\) realizations of \(X(r)\) and \(X(q)\), we could use \[ \Cov{X(r), X(q)} \approx \frac{1}{m}\sum_{i=1}^{m}{(X^{(i)}(r)-\TrueRegFunc(r))(X^{(i)}(q) - \TrueRegFunc(q))} \]

• but we see \(X(r)\) and \(X(q)\) only once

• To get anywhere, we need to assume something about how the covariance function works
  • Useful assumptions are restrictions on \(\Cov{X(r), X(q)}\)
  • We could assume some covariances are exactly equal, so we can pool the corresponding pairs of points in estimating the covariance
  • Or we could assume \(\Cov{X(r), X(q)}\) is a smooth function of \(r\) and \(q\), and use smoothing
  • Or both

Some restrictions on \(\Cov{X(r), X(q)}\)

• Stationarity: for any displacement vector \(h\), \[ \Cov{X(r), X(r+h)} = \Cov{X(q), X(q+h)} \equiv \gamma(h) \]
  • Often also assume \(\TrueRegFunc(q) =\) constant
  • Lets us pool pairs of points separated by the same displacement
• Can be checked:
  • Sub-divide data in space
  • Re-estimate the covariance function in each spatial region
  • Do they match?

Some restrictions on \(\Cov{X(r), X(q)}\)

• Isotropy (in 2 or 3D): if \(\|h\| = \|k\|\), then \[ \Cov{X(r), X(r+h)} = \Cov{X(r), X(r+k)} \]
  • Lets us pool pairs of points separated by the same distance from \(r\)
  • Can also be stationary and isotropic
  • That lets us pool pairs of points at the same distance
• Can be checked:
  • Sub-divide data into pairs of points separated in the same (or similar) directions
  • Re-estimate the covariance function along each direction
  • Do they match?

Some restrictions on \(\Cov{X(r), X(q)}\)

• Separability: In 2D, \(r=(r_1, r_2)\) and \(q=(q_1, q_2)\), \[ \Cov{X(r), X(q)} = C_1(r_1, q_1) C_2(r_2, q_2) \]
  • Similarly for 3D
  • Lets us pool pairs of points with same \(r_1, q_1\) (for \(C_1\))
  • Can also be stationary and separable
  • Lets us pool pairs of points with same \(r_1 - q_1\) (for \(C_1\))
  • Often unphysical but assumed for statistical convenience
  • Is separable and isotropic possible?
• Can be checked:
  • How?

Estimating a stationary covariance function

• For each \(h\), \(n_h =\) number of \(r_i, r_j\) pairs with \(r_i - r_j = h\) \[ \hat{\gamma}(h) = \frac{1}{n_h} \sum_{i, j ~ : ~ r_i - r_j = h}{(X(r_i) - \EstRegFunc(r_i)) (X(r_j) - \EstRegFunc(r_j))} \]

• Works best on a regular grid…
• … which we don’t always have

Think of the Irish wind data

What about Belfast?

(Why is Belfast not in the data set?)

Try a stationary, isotropic covariance

• Assume \(\TrueRegFunc =\) constant so \(\EstRegFunc= 5.26\).
• Assume \(\Cov{X(r), X(q)} = \gamma(\|r-q\|)\)
• Plot \((X(r_i) - \EstRegFunc(r_i)) (X(r_j) - \EstRegFunc(r_j))\) against \(r_i - r_j\)

global.mean <- mean(wind.loc$MeanWind)
distances <- spDists(coordinates(wind.loc), longlat = TRUE)
c.rq <- outer(wind.loc$MeanWind - global.mean, wind.loc$MeanWind - global.mean, "*")
rho.rq <- c.rq/var(wind.loc$MeanWind)

Try a stationary, isotropic covariance

plot(x = as.vector(distances), y = as.vector(rho.rq), xlab = "Distance (km)", ylab = expression((x(r) -
    bar(x))(x(q) - bar(x))/sigma^2))
lines(smooth.spline(x = as.vector(distances), y = as.vector(rho.rq), cv = TRUE))
L.rough <- -log(0.14)/100
curve(exp(-L.rough * x), add = TRUE, col = "blue")
legend("topright", legend = c("Data", "Spline", "Exponential"), pch = c(1, NA, NA),
    lty = c("blank", "solid", "solid"), col = c("black", "black", "blue"))

Try a stationary, isotropic covariance

• Fit the spline first to guide the eye
• Tried exponential decay because that often works
  • \(\rho(\|h\|) = e^{-\|h\| \lambda}\), \(1/\lambda=\) correlation length
  • Rough guess at \(\lambda\) from how much the spline had fallen in 100km
  • In this case, \(1/\lambda \approx 51 \mathrm{km}\)
  • Could improve this through fitting the exponential model (see back-up slides on nonlinear least squares)
  • Sometimes try \(\rho(\|h\|) = c\mathbf{1}(h=0) + (1-c)e^{-\|h\| \lambda}\), nugget effect
    • This looks weird and has a silly name but there’s a reason for it which we’ll come to next week

Kriging

wind.prediction <- function(r, L) {
    fit <- c(fit = NA, se = NA)
    global.mean <- mean(wind.loc$MeanWind)
    global.var <- var(wind.loc$MeanWind)
    distances.to.r <- spDistsN1(pts = coordinates(wind.loc), pt = r, longlat = TRUE)
    distances.among.q <- spDists(coordinates(wind.loc), longlat = TRUE)
    corYZ <- exp(-L * distances.to.r)
    corZ <- exp(-L * distances.among.q)
    Z <- wind.loc$MeanWind
    fit["fit"] <- global.mean + (Z - global.mean) %*% solve(corZ) %*% corYZ
    fit["se"] <- sqrt(global.var - t(corYZ) %*% solve(corZ) %*% corYZ)
    return(fit)
}

wind.prediction(r = as.numeric(belfast.longlat), L = L.rough)

##      fit       se 
## 5.362427 1.358965

Spatio-temporal prediction

• Once more, with feeling
  • The math is exactly the same
• Again: all about estimating covariance functions
• Stationarity: over space? over time? both?
  • If stationary over both space and time, \(\Cov{X(r,t), X(q,s)} = \Cov{X(0,0), X(r-q, s-t)} = \gamma(r-q, s-t)\)
  • Taylor hypothesis (Gneiting, Genton, and Guttorp 2007): for some velocity \(v\), \(\gamma(0,t) = \gamma(vt,0)\)
    • (after G. I. Taylor, see Batchelor (1996); often useful in fluid flows)
• Isotropy: usually only makes sense for space, not for space+time
  • \(\|h\| = \|k\|\) \(\Rightarrow\) \(\Cov{X(r,t), X(r+h, s)} = \Cov{X(q,t), X(q+k, s)}\)
• Separability: space from time? coordinates of space?
  • Space-time separability is often even more artificial than spatial separability

Prediction error

• General result for optimal linear prediction:

\[ Expect{(Y - \LinearExpect{Y|Z})^2} = \Var{Y-\LinearExpect{Y|Z}} = \Var{Y} - \Cov{Y,Z}^T (\Var{Z})^{-1} \Cov{Y,Z} \]

• This is fine if the Oracle tells us the covariance function
• We need to estimate the covariance \(\Rightarrow\) more uncertainty
  • Some approximations, e.g., Zimmerman and Cressie (1992)

Some historical notes: Krige

• The technology of optimal linear prediction we saw for time series was independently re-invented, for purely spatial data, in geology by D. G. Krige in 1951
  • Krige was a mining engineer in South Africa who wanted to estimate the mineral value of new sites from old sites
    • This was his MS thesis!
  • Reference: Krige (1981)
  • Krige does not appear to have known about Wiener, Kolmogorov, etc.

Some figures from Krige (1981)

“Typical gold inch-dwt [=concentration] trend surfaces in the Klerksdorp goldfield on the basis of two-dimensional moving averages [= kriging]… Moving averages of \(100 \times 100\) ft areas within a mined-out section of \(500\times 500\) ft” (p. 28, caption of fig. 21a)

Some figures from Krige (1981)

“Illustration of the agreement between auto-covariance patterns of gold values for 25, 50, 100, and 200-ft square ore units for a section of the Hartebeestfontein mine” (p. 23, caption of fig. 18)

Some historical notes: “kriging”, “geostatistics”

• Krige’s work was taken up by French mathematicians and engineers in the 1960s and 1970s
• The link was Georges Matheron, a mathematician and mining engineer who was the chief scientist for the mining operations in Algeria when it was a French colony (Matheron 2019)
• “Kriging” = French mathematicians and engineers inventing an English verb for “using optimal linear prediction on spatial data”
• “Geostatistics” = French mathematicians and engineers getting really into covariance functions on \(\mathbb{R}^2\) and/or \(\mathbb{S}^2\)
  • A peculiar jargon, and some distinctive ideas, developed in parallel to the rest of statistics
    • Ripley (1981), ch. 4, was an early, influential, and still very useful translation into mainstream statistical terms

Summary

• Spatial and spatio-temporal prediction works just like time series
• The only trick is estimating covariances

Backup: Valid covariance functions

• \(\Cov{X(r), X(q)}\) can’t be just any function of \(r\) and \(q\)
  • \(\Cov{X(r), X(r)} \geq 0\), all implied correlations \(\in [-1, 1]\), etc.
• For every valid \(\Cov{X(r), X(q)}\), for any set of points \(r_1, r_2, \ldots r_n\), the matrix \(c_{ij} = \Cov{X(r_i), X(r_j)}\) is non-negative definite
  • i.e., \(\mathbf{v}^T \mathbf{c} \mathbf{v} \geq 0\) for any vector \(\mathbf{v}\)
• Any function of \(r\) and \(q\) which always leads to non-negative definite matrices is the covariance function of some distribution
• The point of the French mathematicians’ work on geostatistics in the 1960s and 1970s was to realize this, and to use it to restrict the search for covariance functions to ones which are non-negative definite
  • Gneiting, Genton, and Guttorp (2007) has a good discussion, but needs Fourier analysis

Backup: Varieties of kriging

The geostatistics people distinguished a lot of varieties of kriging:

• Simple kriging: assumes weak stationarity, with known mean and covariance function
• Ordinary kriging: assumes weak stationarity, estimate the mean by generalized least squares, still taking the covariance function as known
• Universal kriging; assumes a translation-invariant covariance function but the mean varies spatially, and is a linear function of some measured spatial variable; estimate the linear model by generalized least squares
• etc., etc.

Backup: Nonlinear least squares illustrated

Using optim():

rho.mse <- function(L) {
    mean((rho.rq - exp(-L * distances))^2)
}
L.optim <- optim(par = L.rough, fn = rho.mse, method = "BFGS")
L.optim$par

## [1] 0.02432098

so estimated correlation length is \(41 \mathrm{km}\) rather than initial guess of \(51 \mathrm{km}\)

Backup: Nonlinear least squares illustrated

Using nls():

nls(as.vector(rho.rq) ~ exp(-L * as.vector(distances)), start = list(L = L.rough))

## Nonlinear regression model
##   model: as.vector(rho.rq) ~ exp(-L * as.vector(distances))
##    data: parent.frame()
##       L 
## 0.02428 
##  residual sum-of-squares: 109.6
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 3.649e-06