-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathsi4.qmd
169 lines (139 loc) · 5.92 KB
/
si4.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
---
title: Unknown, varying mean
format: html
engine: knitr
filters:
- webr
---
$$
\newcommand{\E}{{\rm E}} % E expectation operator
\newcommand{\Var}{{\rm Var}} % Var variance operator
\newcommand{\Cov}{{\rm Cov}} % Cov covariance operator
\newcommand{\Cor}{{\rm Corr}}
$$
For this, we need to know how the mean varies. Suppose we model this as
a linear regression model in $p$ known predictors:
$$Z(s_i) = \sum_{j=0}^p \beta_j X_j(s_i) + e(s_i)$$
$$Z(s) = \sum_{j=0}^p \beta_j X_j(s) + e(s) = X(s)\beta + e(s)$$
with $X(s)$ the matrix with predictors,
and row $i$ and column $j$ containing $X_j(s_i)$, and with $\beta =
(\beta_0,...\beta_p)$. Usually, the first column of $X$ contains zeroes
in which case $\beta_0$ is an intercept.
Predictor:
$$\hat{Z}(s_0) = x(s_0)\hat{\beta} + v'V^{-1} (Z-X\hat{\beta}) $$
with $x(s_0) = (X_0(s_0),...,X_p(s_0))$ and $\hat{\beta} = (X'V^{-1}X)^{-1} X'V^{-1}Z$
it has prediction error variance
$$\sigma^2(s_0) = \sigma^2_0 - v'V^{-1}v + Q$$
with $Q = (x(s_0) - X'V^{-1}v)'(X'V^{-1}X)^{-1}(x(s_0) - X'V^{-1}v)$
This form is called external drift kriging, universal kriging or sometimes
regression kriging.
Example in `meuse` data set: `log(zinc)` depending on
`sqrt(meuse)`
```{webr-r}
webr::install("sf")
library(sf) # st_distance
cov = function(h) exp(-h)
calc_beta = function(X, V, z) {
XtVinv = t(solve(V, X))
solve(XtVinv %*% X, XtVinv %*% z)
}
uk = function(data, newdata, X, x0, cov, beta) {
V = cov(st_distance(data))
v = cov(st_distance(data, newdata))
z = data[[1]]
if (missing(beta))
beta = calc_beta(X, V, z)
mu = X %*% beta
x0 %*% beta + t(v) %*% solve(V, z - mu)
}
# prediction location at (0,1):
newdata = st_as_sf(data.frame(x = 0, y = 1), coords = c("x", "y"))
# observation location at (1,1), with attribute value (y) 3:
data = st_as_sf(data.frame(x = 1, y = 1, z = 3), coords = c("x", "y"))
x0 = matrix(1, 1, 1)
X = x0
uk(data, newdata, X, x0, cov)
# three observations location, with attribute values (y) 3,2,5:
data = st_as_sf(data.frame(x = c(1,2,3), y = c(1,1,1), z = c(3, 2, 5)),
coords = c("x", "y"))
newdata = st_as_sf(data.frame(x = .1 * 0:20, y = 1), coords = c("x", "y"))
X = matrix(1,3,1)
x0 = matrix(1, nrow(newdata), 1)
uk(data, newdata, X, x0, cov) # mu = unknown
```
Plotting them:
```{webr-r}
newdata = st_as_sf(data.frame(x = seq(-4,6,by=.1), y = 1), coords = c("x", "y"))
x0 = matrix(1, nrow(newdata), 1)
z_hat = uk(data, newdata, X, x0, cov) # mu unknown
plot(st_coordinates(newdata)[,1], z_hat, type='l', ylim = c(1,5))
points(st_coordinates(data)[,1], data$z, col='red', pch=16)
beta = calc_beta(X, cov(st_distance(data)), data[[1]])
beta
abline(beta, 0, col = 'blue', lty = 2)
```
### Linear trend:
```{webr-r}
X = cbind(matrix(1,3,1), 1:3)
X
xcoord = seq(-4,6,by=.1)
newdata = st_as_sf(data.frame(x = xcoord, y = 1), coords = c("x", "y"))
x0 = cbind(matrix(1, nrow(newdata), 1), xcoord)
z_hat = uk(data, newdata, X, x0, cov) # mu unknown
plot(st_coordinates(newdata)[,1], z_hat, type='l')
points(st_coordinates(data)[,1], data$z, col='red', pch=16)
beta = calc_beta(X, cov(st_distance(data)), data[[1]])
beta
abline(beta, col = 'blue', lty = 2)
```
### With Gaussian covariance:
```{webr-r}
cov = function(h) exp(-(h^2))
z_hat = uk(data, newdata, X, x0, cov) # mu unknown
plot(st_coordinates(newdata)[,1], z_hat, type='l')
points(st_coordinates(data)[,1], data$z, col='red', pch=16)
beta = calc_beta(X, cov(st_distance(data)), data[[1]])
beta
abline(beta, col = 'blue', lty = 2)
```
### Estimating spatial correlation under the UK model
As opposed to the ordinary kriging model, the universal kriging
model needs knowledge of the mean vector in order to estimate
the semivariance (or covariance) from the residual vector:
$$\hat{e}(s) = Z(s) - X\hat\beta$$
but how to get $\hat\beta$ without knowing $V$? This is a chicken-egg
problem. The simplest, but not best, solution is to plug $\hat{\beta}_{OLS}$
in, and from the $e_{OLS}(s)$, estimate $V$ (i.e., the variogram of $Z$)
### Spatial Prediction
... involves errors, uncertainties
### Kriging varieties
* Simple kriging: $Z(s)=\mu+e(s)$, $\mu$ known
* Ordinary kriging: $Z(s)=m+e(s)$, $m$ unknown
* Universal kriging: $Z(s)=X\beta+e(s)$, $\beta$ unknown
* SK: linear predictor $\lambda'Z$ with $\lambda$ such that
$\sigma^2(s_0) = \E(Z(s_0)-\lambda'Z)^2$ is minimized
* OK: linear predictor $\lambda'Z$ with $\lambda$ such that it
* has minimum variance $\sigma^2(s_0) = \E(Z(s_0)-\lambda'Z)^2$, and
* is unbiased $\E(\lambda'Z) = m$
* second constraint: $\sum_{i=1}^n \lambda_i = 1$, weights sum to one.
* UK:
$$\hat{Z}(s_0) = x(s_0)\hat{\beta} + v'V^{-1} (Z-X\hat{\beta}) $$
with $x(s_0) = (X_0(s_0),...,X_p(s_0))$ and $\hat{\beta} = (X'V^{-1}X)^{-1} X'V^{-1}Z$
$$\sigma^2(s_0) = \sigma^2_0 - v'V^{-1}v + Q$$
with $Q = (x(s_0) - X'V^{-1}v)'(X'V^{-1}X)^{-1}(x(s_0) - X'V^{-1}v)$
* OK: fill in a column vector with ones for $X$: $X=(1,1,...,1)'$ and $X_0=1$
* SK: take out the trend/unknown mean
### UK and linear regression
If $Z$ has no spatial correlation, all covariances are zero and
$v=0$ and $V=\mbox{diag}(\sigma^2)$. This implies that
$$\hat{Z}(s_0) = x(s_0)\hat{\beta} + v'V^{-1} (Z-X\hat{\beta}) $$
with $\hat{\beta} = (X'V^{-1}X)^{-1} X'V^{-1}Z$ reduces to
$$\hat{Z}(s_0) = x(s_0)\hat{\beta}$$
with $\hat{\beta} = (X'X)^{-1} X'Z$, i.e., ordinary least squares regression.
Note that
* under this model the residual does not carry information, as it is white noise
* in spatial prediction, UK can not be worse than linear regression, as linear regression is a limiting case of a more general model.
## Application
For an application of universal kriging in mapping of air quality variables, and a comparison
to linear regression, see
* Rob Beelen, Gerard Hoek, Edzer Pebesma, Danielle Vienneau, Kees de Hoogh, David J. Briggs, Mapping of background air pollution at a fine spatial scale across the European Union, Science of The Total Environment, Volume 407, Issue 6, 2009, Pages 1852-1867, https://doi.org/10.1016/j.scitotenv.2008.11.048