---
title: "Week 10b"
format: pdf
editor: source
editor_options:
chunk_output_type: console
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
library(tidyverse)
library(sp)
library(knitr)
library(gstat)
```
### Last Time
- Fitting GP models
- Kriging with plug-in estimates
### This week
- Stationarity
- Locations close are more similar, but how? -> Variograms
---
## Stationarity
Assume the spatial process, $\boldsymbol{Y}(\boldsymbol{s})$ has a mean, ${\mu}(\boldsymbol{s})$, and that the variance of $\boldsymbol{Y}(\boldsymbol{s})$ exists everywhere.
The process is __strictly stationary__ if: for any $n \geq 1$, any set of $n$ sites $\{\boldsymbol{s}_1, \dots, \boldsymbol{s}_n\}$ and any $\boldsymbol{h} \in \mathcal{R}^r$ (typically r =2), the distribution of $(\boldsymbol{Y}(\boldsymbol{s}_1). \dots, \boldsymbol{Y}(\boldsymbol{s}_n))$ is the same as $(\boldsymbol{Y}(\boldsymbol{s_1+h}) \dots, \boldsymbol{Y}(\boldsymbol{s_n + h})$
Weak stationarity, or second-order stationarity of a spatial process, requires a constant mean and covariance that is a function of $\boldsymbol{h}$, $Cov(\boldsymbol{Y}(\boldsymbol{s}),\boldsymbol{Y}(\boldsymbol{s +h})) = C(\boldsymbol{h}),$ where $C(\boldsymbol{h})$ is a covariance function that only requires the separation vector $\boldsymbol{h}.$
Typically with spatial models, the process is assumed to be mean zero as covariates explain the mean structure. So second-order stationarity is primarily focused on the covariance structure.
A third kind of stationarity, known as intrinsic stationarity, assumes $E[\boldsymbol{Y}(\boldsymbol{s + h}) -\boldsymbol{Y}(\boldsymbol{s}) ] =0$ then define
$$E[\boldsymbol{Y}(\boldsymbol{s + h}) - \boldsymbol{Y}(\boldsymbol{s})]^2 = Var(\boldsymbol{Y}(\boldsymbol{s + h}) - \boldsymbol{Y}(\boldsymbol{s})) = 2 \gamma(\boldsymbol{h}),$$
which only works, and satisfies intrinsic stationarity, if the equation only depends on $\boldsymbol{h}$.
## Variograms
Variograms, defined as $2\gamma( \boldsymbol{h})$ are often used to visualize spatial patterns:
\vfill
- If the distance between points is small, the variogram is expected
\vfill
- As the distance between points increases,
\vfill
There is a mathematical link between the covariance function $C(\boldsymbol{h})$ and the variogram $2\gamma(\boldsymbol{h})$.
## Linear semivariogram
$$\gamma(d)=
\begin{cases}
\tau^2 + \sigma^2d \; \; \text{if } d > 0\\
0 \; \; \text{otherwise}
\end{cases}$$
```{r}
tau.sq <- 1
sigma.sq <- 1
d <- seq(0,3, by =.01)
lin.gam <- tau.sq + sigma.sq * d
lin.var <- data.frame(d=d, lin.gam = lin.gam )
ggplot(data = lin.var, aes(x=d, y=lin.gam)) +
geom_line() + ylim(0,4) + ylab('linear variogram') + theme_bw()
```
### Nugget, sill, partial-sill, and range
Nugget is defined as $\gamma(d): \; d \rightarrow 0_+$
\vfill
Sill is defined as $\gamma(d): \; d \rightarrow \infty$
\vfill
Partial sill is defined as the $\text{sill} - \text{nugget}$
\vfill
Range is defined as the first point where $\gamma(d)$ reaches the sill. Range is sometimes (1/$\phi$) and sometimes $\phi$ depending on the paremeterization.
## Spherical semivariogram:
$$\gamma(d)=
\begin{cases}
\tau^2 + \sigma^2 \; \; \text{if } d \geq 1/ \phi\\
\tau^2 + \sigma^2\left[\frac{3\phi d}{2} - \frac{1}{2} (\phi d)^3 \right]\\
0 \; \; \text{otherwise}
\end{cases}$$
- Sketch, or generate in R, a spherical semivariogram
- On this figure label the nugget, sill, partial sill, and range.
## Exponential
We saw the exponential covariance earlier in class, what is the mathematical form of the covariance?
\vfill
\vfill
The variogram is
$$\gamma(d)=
\begin{cases}
\tau^2 + \sigma^2(1 - \exp(- d / \phi)) \; \; \text{if } d > 0\\
0 \; \; \text{otherwise}
\end{cases}$$
```{r}
#| echo: false
range <- 2
phi <- 1/range
exp.gam <- tau.sq + sigma.sq * (1 - exp( - phi * d))
exp.var <- data.frame(d=d, exp.gam = exp.gam )
ggplot(data = exp.var, aes(x=d, y=exp.gam)) + geom_line() + ylim(0,4) +
ylab('exponential variogram') +
annotate('text', x=0.1, y = 0.9, label = 'nugget') +
# annotate('text', x=2.9, y = 2.65, label = 'sill?') +
# annotate('text', x=2.2, y = .1, label = 'range?', color = 'blue') +
annotate("segment", x = 2, xend = 2, y = 0, yend = 4,
colour = "blue", linetype =3, size = .5) + theme_bw()
```
\vfill
The *effective range* is defined as the distance where there is *effectively* no spatial structure. Generally this is determined when the correlation is .05.
\vfill
\begin{eqnarray*}
\exp(-d_o / \phi) &=&.05 \\
d_0 / \phi &=& - \log(.05) \\
d_0 &\approx& 3 \phi
\end{eqnarray*}
\vfill
## More Semivariograms: Equations
Gaussian:
$$\gamma(d)=
\begin{cases}
\tau^2 + \sigma^2(1 - \exp(-\phi^2 d^2)) \; \; \text{if } d > 0\\
0 \; \; \text{otherwise}
\end{cases}$$
\vfill
Powered Exponential:
$$\gamma(d)=
\begin{cases}
\tau^2 + \sigma^2(1 - \exp(-|\phi d|^p)) \; \; \text{if } d > 0\\
0 \; \; \text{otherwise}
\end{cases}$$
\vfill
Mat$\acute{\text{e}}$rn:
$$\gamma(d)=
\begin{cases}
\tau^2 + \sigma^2\left[1 - \frac{(2\sqrt{\nu}d \phi)^\nu}{2^{\nu-1}\Gamma(\nu)} K_\nu (2 \sqrt{\nu} d \phi)\right] \; \; \text{if } d > 0\\
0 \; \; \text{otherwise}
\end{cases},$$
where $K_\nu$ is a modified Bessel function and $\Gamma()$ is a Gamma function.
\vfill
### Variogram Creation: How?
```{r, echo=T,}
data(meuse)
meuse.small <- meuse %>% select(x, y, copper) %>% as_tibble()
meuse.small %>% head(15) %>% kable()
meuse |>
ggplot(aes(x = x , y=y, size = copper, color = copper)) +
geom_point() +
scale_color_viridis_b() +
theme_minimal()
```
### Variogram Creation: Steps
1. Calculate distances between sampling locations
2. Choose grid for distance calculations
3. Calculate empirical semivariogram
$$\hat{\gamma}(d_k) = \frac{1}{2N(d_k)} \sum_{\boldsymbol{s}_i,\boldsymbol{s}_i, \in N(d_k)}\left[ \boldsymbol{Y}(\boldsymbol{s}_i) - \boldsymbol{Y}(\boldsymbol{s}_j) \right]^2,$$
where $$N(d_k) = \{(\boldsymbol{s}_i,\boldsymbol{s}_j): ||\boldsymbol{s}_i -\boldsymbol{s}_j|| \in I_k\}$$ and $I_k$ is the $k^{th}$ interval.
4. Plot the semivariogram