--- title: "Spatial Statistics Notation" output: pdf_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = F, message = F) library(tidyverse) library(knitr) library(gstat) library(sp) set.seed(02192021) ``` \vfill \vfill \vfill \vfill \vfill \vfill \vfill \newpage Given that we assume $$\boldsymbol{Y}|\mu, \underline{\theta} \sim MVN_n(\underline{\mu}, \Sigma(\underline{\theta})),$$ \vfill \vfill \vfill ## Stationarity Assume the spatial process, $\boldsymbol{Y}(\boldsymbol{s})$ has a mean, $\underline{\mu}(\boldsymbol{s})$, and that the variance of $\boldsymbol{Y}(\boldsymbol{s})$ exists everywhere. \vfill ## Strict Stationarity Is this strictly stationary? ```{r} # Data a <- data.frame( x=rnorm(20000, 0, .5), y=rnorm(20000, 1, .5) ) b <- data.frame( x=rnorm(20000, 2, 1), y=rnorm(20000, -.5, 1) ) c <- data.frame( x=rnorm(20000, -.5, 1), y=rnorm(20000, 0, 1) ) data <- rbind(a,b,c) ggplot(data, aes(x=x, y=y) ) + stat_density_2d(aes(fill = ..density..), geom = "raster", contour = FALSE) + scale_fill_distiller(palette= "Spectral", direction=1) + scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0)) + theme( legend.position='none' ) ``` \newpage ## Weak Stationarity 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}),$ \vfill Typically with spatial models, the process is assumed to be mean zero as covariates explain the mean structure. ## Weak Stationarity Is this weakly stationary? ```{r} # Data dat <- rep(0,100) dat[1:100 %%2 == 1] <- sample(x = c(1,-1), size = 50, replace =T) dat[1:100 %%2 == 0] <- rnorm(50) comb <- data.frame(dat = dat, n = 1:100) ggplot(comb, aes(x=n, y=dat) ) + geom_point() + geom_line() + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) + theme_bw() ``` \newpage ## Intrinsic Stationarity A third kind of stationarity, known as intrinsic stationarity, describes the behavior of differences in the spatial process, rather than the data. \vfill 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}$. \newpage ## Variograms Intrinsic stationary justifies the use of variograms $2\gamma( \boldsymbol{h})$ \vfill Variograms are often used to visualize spatial patterns: \vfill \vfill \vfill There is a mathematical link between the covariance function $C(\boldsymbol{h})$ and the variogram $2\gamma(\boldsymbol{h})$. ## Isotropy \vfill \vfill ## Variogram to the Covariance Function \begin{eqnarray*} 2 \gamma (\boldsymbol{h}) &=& Var(\boldsymbol{Y}(\boldsymbol{s+h}) - \boldsymbol{Y}(\boldsymbol{s}))\\ &=& Var(\boldsymbol{Y}(\boldsymbol{s+h})) + Var(\boldsymbol{Y}(\boldsymbol{s})) - 2 Cov(Var(\boldsymbol{Y}(\boldsymbol{s+h}),\boldsymbol{Y}(\boldsymbol{s})))\\ &=&C(\boldsymbol{0}) + C(\boldsymbol{0}) - 2C(\boldsymbol{h})\\ &=& 2 \left[C(\boldsymbol{0}) - C(\boldsymbol{h})\right] \end{eqnarray*} Given $C()$, the variogram can easily be recovered \newpage ## 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 \vfill \vfill \vfill \newpage ## 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. ```{r, eval = T} tau.sq <- 1 sigma.sq <- 2 range <- 2 phi <- 1/range d <- seq(0,3, by =.01) sph.gam <- tau.sq + sigma.sq * (3 * phi * d / 2 - .5 * (phi*d)^3) sph.gam[d >= range] <- tau.sq + sigma.sq sph.var <- data.frame(d=d, sph.gam = sph.gam ) # ggplot(data = sph.var, aes(x=d, y=sph.gam)) + geom_line() + ylim(0,4) + ylab('spherical variogram') + annotate('text', x=0.1, y = 0.9, label = 'nugget') + annotate('text', x=2.5, y = 3.1, 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) + annotate("segment", x = 2, xend = 3, y = 1, yend = 1,colour = "red", linetype =3, size = .5) + annotate("segment", x = 2, xend = 3, y = 3, yend = 3,colour = "red", linetype =3, size = .5) + annotate("segment", x = 2.5, xend = 2.5, y = 1, yend = 3,colour = "red", linetype =3, size = .5) + annotate('text', x=2.2, y = 2, label = 'partial sill', color = 'red') ``` \newpage ## 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} 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 \vfill \vfill ```{r} d <- seq(0,8, by = .01) 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=7.9, y = 3.1, label = 'sill') + annotate('text', x=6.6, y = .1, label = 'effective range', color = 'blue') + annotate("segment", x = 6, xend = 6, y = 0, yend = 4, colour = "blue", linetype =3, size = .5) + theme_bw() ``` \vfill \newpage ## 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 \newpage ## Variogram Creation: How? ```{r, echo=T,} data(meuse) meuse.small <- meuse %>% select(x, y, copper) %>% as_tibble() meuse.small %>% head(15) %>% kable() ``` ## 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 ## Variogram Creation: R function ```{r, eval = T} coordinates(meuse) = ~x+y variogram(copper~1, meuse) %>% plot() ``` ## Variogram Creation: How - Step 1 Calculate Distances between sampling locations ```{r, echo= T} #dist(meuse.small) dist.mat <- dist(meuse.small %>% select(x,y)) ``` ## Variogram Creation: How - Step 2 Choose grid for distance calculations ```{r, echo= T} cutoff <- max(dist.mat) / 3 # default maximum distance num.bins <- 15 bin.width <- cutoff / 15 ``` ## Variogram Creation: How - Step 3 Calculate empirical semivariogram ```{r, echo = T} dist.sq <- dist(meuse.small$copper)^2 vario.dat <- data.frame(dist = as.numeric(dist.mat), diff = as.numeric(dist.sq)) %>% mutate(bin = floor(dist / bin.width) + 1) %>% filter(bin < 16) %>% group_by(bin) %>% summarize(emp.sv = .5 * mean(diff)) vario.dat ``` ## Variogram Creation: How - Step 4 Plot empirical semivariogram ```{r} vario.dat %>% ggplot(aes(x=bin, y= emp.sv)) + geom_point() + xlim(0,15) + ylim(0,800) + scale_x_continuous(breaks=c(5, 10, 15), labels =c('500','1000','1500')) + ggtitle('Empirical Semivariogram for Copper in Meuse Data Set') + xlab('distance (m)') + ylab('Empirical Semivariogram') + theme_bw() ```