--- title: Linear Prediction for Spatial and Spatio-Temporal Fields date: 10 September 2026 (Lecture 6) author: 36-740, Fall 2026 output: slidy_presentation: math_method: engine: mathjax url: https://bactra.org/mathjax/tex-svg.js bibliography: locusts.bib --- ```{r, include=FALSE} # General set-up options library(knitr) opts_chunk$set(size="small", background="white", highlight=FALSE, cache=TRUE, autodep=TRUE, tidy=TRUE, tidy.opts=list(comment=FALSE), warning=FALSE, message=FALSE, echo=FALSE) ``` \[ \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? ## Optimal linear prediction for spatial data {.smaller} - Given: $X(r_1), X(r_2), \ldots X(r_n)$ - Desired: prediction of $X(r_0)$ \begin{eqnarray} \EstRegFunc(r_0) & = & \alpha + \vec{\beta} \cdot \left[\begin{array}{c} X(r_1) \\ X(r_2) \\ \vdots \\ X(r_n) \end{array}\right]\\ \alpha & = & \Expect{X(r_0)} - \vec{\beta} \cdot \left[\begin{array}{c} \Expect{X(r_1)}\\ \Expect{X(r_2)} \\ \vdots \\ \Expect{X(r_n)}\end{array}\right] ~ \text{(goes away if everything's centered)}\\ \vec{\beta} & = & {\left[\begin{array}{cccc} \Var{X(r_1)} & \Cov{X(r_1), X(r_2)} & \ldots & \Cov{X(r_1), X(r_n)}\\ \Cov{X(r_1), X(r_2)} & \Var{X(r_2)} & \ldots & \Cov{X(r_2), X(r_n)}\\ \vdots & \vdots & \ldots & \vdots\\ \Cov{X(r_1), X(r_n)} & \Cov{X(r_2), X(r_n)} & \ldots & \Var{X(r_n)}\end{array}\right]}^{-1} \left[\begin{array}{c} \Cov{X(r_0), X(r_1)}\\ \Cov{X(r_0), X(r_2)}\\ \vdots \\ \Cov{X(r_0), X(r_n)}\end{array}\right] \end{eqnarray} ## 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 ## Smoothing for covariance estimation {.smaller} - Two (equivalent) definitions of the covariance: \[ \Cov{X(r), X(s)} = \Expect{X(r) X(s)} - \Expect{X(r)}\Expect{X(s)} = \Expect{(X(r) - \Expect{X(t)}) (X(s) - \Expect{X(s)})} \] - Use the 2nd form: _define_ \[ \Gamma(r,s) \equiv (X(r) - \Expect{X(r)})(X(s) - \Expect{X(s)}) \] so $\Cov{X(r), X(s)} = \Expect{\Gamma(r,s)}$ - Now assume constant mean so we can estimate \[ \widehat{\Gamma}(r,s) = (X(r) - \overline{x}) (X(s) - \overline{x}) \] + Or substitute in trend estimates of your choice - Now assume stationarity so $\Cov{X(r), X(s)} = \Expect{\Gamma(r,s)} = \AutoCov(r-s)$ - Finally assume $\AutoCov(h)$ changes slowly in $h$ - We can _estimate_ $\AutoCov(h)$ by finding pairs $r, s$ with $r-s \approx h$ and averaging $\widehat{\Gamma}(r,s)$ + Or kernel smoothing or spline smoothing or... - **Warning**: not _necessarily_ a valid covariance function ## 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 ```{r} # Setting up the wind data for nice display, after help(wind) library(gstat) data(wind) library(sp) # char2dms wind.loc$y = as.numeric(char2dms(as.character(wind.loc[["Latitude"]]))) wind.loc$x = as.numeric(char2dms(as.character(wind.loc[["Longitude"]]))) coordinates(wind.loc) = ~x+y library(mapdata) map("worldHires", xlim = c(-11,-5.4), ylim = c(51,55.5)) plot(wind.loc, add=TRUE, pch=16, cex=2, col=topo.colors(5)[cut(wind.loc$MeanWind, breaks=5)]) text(coordinates(wind.loc), pos=1, label=wind.loc$Station) legend("topleft", legend=levels(cut(wind.loc$MeanWind, breaks=5)), col=topo.colors(5), pch=16, cex=0.5) ``` ## What about Belfast? ```{r} map("worldHires", xlim = c(-11,-5.4), ylim = c(51,55.5)) plot(wind.loc, add=TRUE, pch=16, cex=2, col=topo.colors(5)[cut(wind.loc$MeanWind, breaks=5)]) text(coordinates(wind.loc), pos=1, label=wind.loc$Station) legend("topleft", legend=levels(cut(wind.loc$MeanWind, breaks=5)), col=topo.colors(5), pch=16, cex=0.5) belfast.longlat <- data.frame(x=-5.926437, y=54.607868) points(belfast.longlat, cex=2, pch=16) text(belfast.longlat, pos=1, label="Belfast") ``` (Why is Belfast not in the data set?) ## Try a stationary, isotropic covariance - Assume $\TrueRegFunc =$ constant so $\EstRegFunc= `r signif(mean(wind.loc$MeanWind), 3)`$. - 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$ ```{r, echo=TRUE} # Get the over-all mean of the wind speed global.mean <- mean(wind.loc$MeanWind) # Calculate actual distances (in kilometers) between map points distances <- spDists(coordinates(wind.loc), longlat=TRUE) # Find (x(r) - mean(x)) * (x(q) - mean(x)) for all combinations of points # in the data # expectations of these would equal the covariance so we'll smooth them c.rq <- outer(wind.loc$MeanWind - global.mean, wind.loc$MeanWind - global.mean, "*") # scale down to see the terms going in to correlations # presumes constant variance over space, but we're already committed to that # by stationarity rho.rq <- c.rq/var(wind.loc$MeanWind) ``` ## Try a stationary, isotropic covariance ```{r, echo=TRUE} # Plot distance vs. correlation 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)) # Run a spline smoother through this lines(smooth.spline(x=as.vector(distances), y=as.vector(rho.rq), cv=TRUE)) # exponential correlations often work well --- eye-ball how much the correlation # has fallen in 100 km and back out the exponent L.rough <- -log(0.14)/100 # add the curve for the roughly-fitted exponential correlation function 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 `r signif(1/L.rough, 2)` \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 ```{r, echo=TRUE} # Inputs: vector for point to make prediction at (r) # Inverse correlation length (L) # Output: Predicted value for the wind at r # Presumes: r is a two-entry vector, with first entry longitude, second latitude # wind.loc data frame is loaded # using an exponential correlation function, without nugget wind.prediction <- function(r, L) { # Prepare a little vector to store prediction and its standard error fit <- c(fit=NA, se=NA) global.mean <- mean(wind.loc$MeanWind) global.var <- var(wind.loc$MeanWind) # Find distance from r to the data points distances.to.r <- spDistsN1(pts=coordinates(wind.loc), pt=r, longlat=TRUE) # Ditto among data points distances.among.q <- spDists(coordinates(wind.loc), longlat=TRUE) # Find correlations # All we need under assumption of spatial stationarity corYZ <- exp(-L*distances.to.r) corZ <- exp(-L*distances.among.q) # The data... Z <- wind.loc$MeanWind # Make the prediction fit["fit"] <- global.mean + (Z-global.mean) %*% solve(corZ) %*% corYZ # the standard error fit["se"] <- sqrt(global.var - t(corYZ) %*% solve(corZ) %*% corYZ) return(fit) } wind.prediction(r=as.numeric(belfast.longlat), L=L.rough) ``` ## 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-Guttorp]: for some velocity $v$, $\gamma(0,t) = \gamma(vt,0)$ * (after G. I. Taylor, see @Batchelor-on-Taylor; 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-Cressie-prediction-error-of-spatial-models-with-estimated-parameters ## 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-on-kriging + Krige does not appear to have known about Wiener, Kolmogorov, etc. ## Some figures from @Krige-on-kriging ![](krige-kriging-surface.jpg) "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-on-kriging ![](krige-acf.jpg) "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-regionalised] - "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-spatial-statistics, 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: Estimating the Mean of a Stationary Random Field {.smaller} - I have been using the time- or space- average to estimate the expectation value of stationary processes: $\Expect{X(r)} = \mu$ has been estimated by \[ \overline{x} = \frac{1}{n}\sum_{i=1}^{n}{X(r_i)} = \argmin_{m\in\mathbb{R}}{\frac{1}{n} (\mathbf{X} - m\mathbf{1})^T(\mathbf{X} - m\mathbf{1})} \] + where $\mathbf{1} =$ $n\times 1$ matrix of all 1s - This is _consistent_ (assuming weak stationarity), but not necessarily _efficient_ - If the variance isn't constant, we shouldn't try as hard to match the high-variance observations as the low-variance observations - If $X(r_1)$ and $X(r_2)$ are positively correlated, we shouldn't try as hard to match _both_ of them (compared to their being uncorrelated) - **Generalized least squares**: if $\mathbf{X}$ has variance matrix $\mathbf{v}$, then minimize \[ \frac{1}{n}(\mathbf{X} - m\mathbf{1})^T\mathbf{v}^{-1}(\mathbf{X}-m\mathbf{1}) \] Solution: \[ \hat{m} = (\mathbf{1}^T \mathbf{v}^{-1} \mathbf{1})^{-1} \mathbf{1}^T\mathbf{v}^{-1} \mathbf{X} \] + Check: this is $1\times 1$, as it should be ## 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-Guttorp 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: Parametric covariance (or correlation) estimation {.smaller} - $\Cov{X(r), X(q)} = \Expect{(X(r)-\TrueRegFunc(r))(X(q)-\TrueRegFunc(q))}$ - For each pair of points $r$ and $q$, create the product of deviations $C(r,q) = (X(r) - \EstRegFunc(r))(X(q) - \EstRegFunc(q))$ - If we had many realizations of $X(r), X(q)$, we could average $C(r,q)$ to estimate $\Cov{X(r), X(q)}$ - Assume instead that we have a **parametric** form for the covariance, $\Cov{X(r), X(q)} = f(r,q,\theta)$ (parameter $\theta$ unknown) + E.g., $f(r,q,\lambda, \sigma^2) = \sigma^2 e^{-\|r-q\|\lambda}$ for a stationary, isotropic, exponentially-decaying covariance, with correlation length $1/\lambda$ and variance $\sigma^2$ - Then $C(r,q) = f(r,q,\theta)+$ noise, with noise having mean 0 - Use **nonlinear least squares** (NLS): \[ \hat{\theta} = \argmin_{\theta}{\frac{1}{n^2}\sum_{i=1}^{n}{\sum_{j=1}^{n}{{\left( C(r_i, r_j) - f(r_i, r_j, \theta)\right)}^2}}} \] + Can apply this to correlations rather than covariances if the variances are known, or assumed constant ## Backup: Nonlinear least squares {.smaller} - General setting: $Y=f(Z,\theta) +$ noise + data $(z_1, y_1), (z_2, y_2), \ldots (z_m, y_m)$ - NLS estimate is \[ \hat{\theta} = \argmin_{\theta}{\frac{1}{m}{\sum_{i=1}^{m}{{\left( y_i - f(z_i, \theta)\right)}^2}}} \] - Statistical properties very similar to ordinary least squares, but computationally more annoying + Statistically, optimization methods are all pretty similar, and we'll see the theory in a few weeks + Computationally, no closed-form solution (unlike OLS), so you need to iteratively approach the optimum - In R: either `nls` (specifically for nonlinear least squares) or `optim` (for general-purpose optimization) - Usually needs an initial guess about $\theta$ ## Backup: Nonlinear least squares illustrated Using `optim()`: ```{r, echo=TRUE} # 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 ``` so estimated correlation length is $`r signif(1/L.optim$par, 2)` \mathrm{km}$ rather than initial guess of $`r signif(1/L.rough, 2)` \mathrm{km}$ ## Backup: Nonlinear least squares illustrated Using `nls()`: ```{r, echo=TRUE} # Need to use as.vector() here because both rho.rq and distances are # set up as matrices, but this makes nls() unhappy nls(as.vector(rho.rq) ~ exp(-L*as.vector(distances)), start=list(L=L.rough)) ``` ## Backup: "Variograms" {.smaller} - A lot of the geostatistical literature emphasizes idea of the **variogram** \[ v(r, s) = \frac{1}{2}\Var{X(r) - X(s)} \] - Notice that \[ v(r,s) = \frac{1}{2}\left(\Var{X(r)} + \Var{X(s)} - 2\Cov{X(r), X(s)}\right) \] so if you know the variances, the variogram and the covariance function are equivalent - If variance is constant, $\Var{X(r)} = \sigma^2$, \[ v(r,s) = \sigma^2 - \Cov{X(r), X(s)} \] - Any assumption that constrains the covariance function constrains the variogram and vice-versa, e.g., stationarity, \[ v(r,s) = v(r-s) \] so we can pool data to estimate (There are some situations where $v(r,s)$ is well-defined but $\Cov{X(r), X(s)}$ is not, but they are _weird_) ## Backup: Generalized Least Squares and Gauss-Markov {.smaller} - Given: $\mathbf{y}$ ($n\times 1$), $\mathbf{x}$ ($n \times p$), positive-definite $\mathbf{v}$ ($n\times n$) - We decided to minimize \[ \frac{1}{n}(\mathbf{y} -\mathbf{x} b)^T \mathbf{v}^{-1} (\mathbf{y} - \mathbf{x} b) \] - The solution is \[ \hat{\beta}_{GLS} = (\mathbf{x}^T \mathbf{v}^{-1} \mathbf{x})^{-1} \mathbf{x}^T \mathbf{v}^{-1} \mathbf{y} \] - Assume $\Expect{Y|X} = \beta \mathbf{X} + \epsilon$, and $\Expect{\epsilon} = 0$ - Then $\Expect{\hat{\beta}_{GLS}} = \beta$ (unbiased) - If $\Var{\epsilon} = \mathbf{v}$, then $\hat{\beta}_{GLS}$ is the minimum variance unbiased linear estimator of $\beta$ (Gauss-Markov theorem) + Despite the name, does not require the noise to be Gaussian + Proof: Not given here, but see @CRS-ADAfaEPoV, sec. 10.2.2 ## References