2-Geostatistical_models.Rmd

---
title: "Spatial statistics in ecology, Part 2"
author: "Philippe Marchand, Université du Québec en Abitibi-Témiscamingue"
date: "January 14, 2021"
output: 
    html_document:
        self_contained: false
        lib_dir: libs
        theme: united
        toc: true
        toc_float: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(cowplot)
theme_set(theme_cowplot())
```

# Spatial correlation of a variable

Correlation between measurements of a variable taken at nearby points often occurs in environmental data. This principle is sometimes referred to as the "first law of geography" and is expressed in the following quote from Waldo Tobler: "Everything is related to everything else, but near things are more related than distant things".

In statistics, we often refer to *autocorrelation* as the correlation between measurements of the same variable taken at different times (temporal autocorrelation) or places (spatial autocorrelation).

## Intrinsic or induced dependence

There are two basic types of spatial dependence on a measured variable $y$: an *intrinsic* dependence on $y$, or a dependence *induced* by external variables influencing $y$, which are themselves spatially correlated.

For example, suppose that the abundance of a species is correlated between two sites located near each other:

- this spatial dependence can be induced if it is due to a spatial correlation of habitat factors that are favorable or unfavorable to the species;

- or it can be intrinsic if it is due to the dispersion of individuals to nearby sites.

In many cases, both types of dependence affect a given variable.

If the dependence is simply induced and the external variables that cause it are included in the model explaining $y$, then the model residuals will be independent and we can use all the methods already seen that ignore spatial correlation.

However, if the dependence is intrinsic or due to unmeasured external factors, then the spatial correlation of the residuals in the model will have to be taken into account.

## Different ways to model spatial effects

In this training, we will directly model the spatial correlations of our data. It is useful to compare this approach to other ways of including spatial aspects in a statistical model.

First, we could include predictors in the model that represent position (e.g., longitude, latitude). Such predictors may be useful for detecting a systematic large-scale trend or gradient, whether or not the trend is linear (e.g., with a generalized additive model).

In contrast to this approach, the models we will see now serve to model a spatial correlation in the random fluctuations of a variable (i.e., in the residuals after removing any systematic effect).

Mixed models use random effects to represent the non-independence of data on the basis of their grouping, i.e., after accounting for systematic fixed effects, data from the same group are more similar (their residual variation is correlated) than data from different groups. These groups were sometimes defined according to spatial criteria (observations grouped into sites).

However, in the context of a random group effect, all groups are as different from each other, e.g., two sites within 100 km of each other are no more or less similar than two sites 2 km apart.

The methods we will see here and in the next parts of the training therefore allow us to model non-independence on a continuous scale (closer = more correlated) rather than just discrete (hierarchy of groups).


# Geostatistical models

Geostatistics refers to a group of techniques that originated in the earth sciences. Geostatistics is concerned with variables that are continuously distributed in space and where a number of points are sampled to estimate this distribution. A classic example of these techniques comes from the mining field, where the aim was to create a map of the concentration of ore at a site from samples taken at different points on the site.

For these models, we will assume that $z(x, y)$ is a stationary spatial variable measured at points with coordinates $x$ and $y$.

## Variogram

A central aspect of geostatistics is the estimation of the variogram $\gamma_z$ . The variogram is equal to half the mean square difference between the values of $z$ for two points $(x_i, y_i)$ and $(x_j, y_j)$ separated by a distance $h$.

$$\gamma_z(h) = \frac{1}{2} \text{E} \left[ \left( z(x_i, y_i) - z(x_j, y_j) \right)^2 \right]_{d_{ij} = h}$$

In this equation, the $\text{E}$ function with the index $d_{ij}=h$ designates the statistical expectation (i.e., the mean) of the squared deviation between the values of $z$ for points separated by a distance $h$.

If we want instead to express the autocorrelation $\rho_z(h)$ between measures of $z$ separated by a distance $h$, it is related to the variogram by the equation:

$$\gamma_z = \sigma_z^2(1 - \rho_z)$$ ,

where $\sigma_z^2$ is the global variance of $z$. 

Note that $\gamma_z = \sigma_z^2$ when we reach a distance where the measurements of $z$ are independent, so $\rho_z = 0$. In this case, we can see that $\gamma_z$ is similar to a variance, although it is sometimes called "semivariogram" or "semivariance" because of the 1/2 factor in the above equation.

## Theoretical models for the variogram

Several parametric models have been proposed to represent the spatial correlation as a function of the distance between sampling points. Let us first consider a correlation that decreases exponentially:

$$\rho_z(h) = e^{-h/r}$$

Here, $\rho_z = 1$ for $h = 0$ and the correlation is multiplied by $1/e \approx 0.37$ each time the distance increases by $r$. In this context, $r$ is called the *range* of the correlation.

From the above equation, we can calculate the corresponding variogram.

$$\gamma_z(h) = \sigma_z^2 (1 - e^{-h/r})$$

Here is a graphical representation of this variogram.

```{r, echo = FALSE}
ggplot(NULL) + xlim(0, 10) +
    labs(x = "h", y = expression(gamma)) +
    stat_function(fun = function(x) 5 * (1 - exp(-x/3)),
                  geom = "line", color = "#b3452c", size = 1) +
    geom_hline(yintercept = 5, linetype = "dotted") +
    geom_segment(aes(x = 0, xend = 2.8, y = 5*(1-exp(-1)),
                     yend = 5*(1-exp(-1))), 
                arrow = arrow(length = unit(0.05, "inches"), ends = "both", type = "closed")) +
    annotate("text", x = 1.5, y = 3.5, label = "range") +
    annotate("text", x = 9, y = 5.5, label = "sill") +
    scale_y_continuous(limits = c(0, 6))
```

Because of the exponential function, the value of $\gamma$ at large distances approaches the global variance $\sigma_z^2$ without exactly reaching it. This asymptote is called a *sill* in the geostatistical context and is represented by the symbol $s$.

Finally, it is sometimes unrealistic to assume a perfect correlation when the distance tends towards 0, because of a possible variation of $z$ at a very small scale. A *nugget* effect, denoted $n$, can be added to the model so that $\gamma$ approaches $n$ (rather than 0) if $h$ tends towards 0. The term nugget comes from the mining origin of these techniques, where a nugget could be the source of a sudden small-scale variation in the concentration of a mineral.

By adding the nugget effect, the remainder of the variogram is "compressed" to keep the same sill, resulting in the following equation.

$$\gamma_z(h) = n + (s - n) (1 - e^{-h/r})$$

In the *gstat* package that we use below, the term $(s-n)$ is called a *partial sill* or `psill` for the exponential portion of the variogram.

```{r, echo = FALSE}
ggplot(NULL) + xlim(0, 10) +
    labs(x = "h", y = expression(gamma)) +
    stat_function(fun = function(x) 4 * (1 - exp(-x/3)) + 1,
                  geom = "line", color = "#b3452c", size = 1) +
    geom_hline(yintercept = 5, linetype = "dotted") +
    geom_segment(aes(x = 0, xend = 2.8, y = 4*(1-exp(-1)) + 1,
                     yend = 4*(1-exp(-1)) + 1), 
                arrow = arrow(length = unit(0.05, "inches"), 
                              ends = "both", type = "closed")) +
    geom_segment(aes(x = 0, xend = 0, y = 0, yend = 0.9),
                 arrow = arrow(length = unit(0.05, "inches"),
                               ends = "both", type = "closed")) +
    annotate("text", x = 1.5, y = 4, label = "range") +
    annotate("text", x = 9, y = 5.5, label = "sill") +
    annotate("text", x = 0.5, y = 0.5, label = "nugget") +
    scale_y_continuous(limits = c(0, 6))
```

In addition to the exponential model, two other common theoretical models for the variogram are the Gaussian model (where the correlation follows a half-normal curve), and the spherical model (where the variogram increases linearly at the start and then curves and reaches the plateau at a distance equal to its range $r$). The spherical model thus allows the correlation to be exactly 0 at large distances, rather than gradually approaching zero in the case of the other models.

 Model | $\rho(h)$ | $\gamma(h)$
-------|-----------|-------------
Exponential | $\exp\left(-\frac{h}{r}\right)$ | $s \left(1 - \exp\left(-\frac{h}{r}\right)\right)$
Gaussian | $\exp\left(-\frac{h^2}{r^2}\right)$ | $s \left(1 - \exp\left(-\frac{h^2}{r^2}\right)\right)$
Spherical $(h < r)$ * | $1 - \frac{3}{2}\frac{h}{r} + \frac{1}{2}\frac{h^3}{r^3}$ | $s \left(\frac{3}{2}\frac{h}{r} - \frac{1}{2}\frac{h^3}{r^3} \right)$

\* For the spherical model, $\rho = 0$ and $\gamma = s$ if $h \ge r$.


```{r, echo = FALSE, fig.dim = c(9, 4)}
vexp <- ggplot(NULL) + xlim(0, 10) +
    labs(x = "h", y = expression(gamma), title = "Exponential") +
    stat_function(fun = function(x) 5 * (1 - exp(-x/3)),
                  geom = "line", color = "#b3452c", size = 1) +
    geom_hline(yintercept = 5, linetype = "dotted")
    

vgau <- ggplot(NULL) + xlim(0, 10) +
    labs(x = "h", y = expression(gamma), title = "Gaussian") +
    stat_function(fun = function(x) 5 * (1 - exp(-x^2/4^2)),
                  geom = "line", color = "#b3452c", size = 1) +
    geom_hline(yintercept = 5, linetype = "dotted")

vsph <- ggplot(NULL) + xlim(0, 10) +
    labs(x = "h", y = expression(gamma), title = "Spherical") +
    stat_function(fun = function(x) ifelse(x < 8, 5 * (1.5*x/8 - 0.5*x^3/8^3), 5),
                  geom = "line", color = "#b3452c", size = 1) +
    geom_hline(yintercept = 5, linetype = "dotted")

plot_grid(vexp, vgau, vsph, nrow = 1)
```

## Empirical variogram

To estimate $\gamma_z(h)$ from empirical data, we need to define distance classes, thus grouping different distances within a margin of $\pm \delta$ around a distance $h$, then calculating the mean square deviation for the pairs of points in that distance class.

$$\hat{\gamma_z}(h) = \frac{1}{2 N_{\text{paires}}} \sum \left[ \left( z(x_i, y_i) - z(x_j, y_j) \right)^2 \right]_{d_{ij} = h \pm \delta}$$

We will see in the next section how to estimate a variogram in R.

## Regression model with spatial correlation

The following equation represents a multiple linear regression including residual spatial correlation:

$$v = \beta_0 + \sum_i \beta_i u_i + z + \epsilon$$

Here, $v$ designates the response variable and $u$ the predictors, to avoid confusion with the spatial coordinates $x$ and $y$. 

In addition to the residual $\epsilon$ that is independent between observations, the model includes a term $z$ that represents the spatially correlated portion of the residual variance.

Here are suggested steps to apply this type of model:

1. Fit the regression model without spatial correlation.

2. Verify the presence of spatial correlation from the empirical variogram of the residuals.

3. Fit one or more regression models with spatial correlation and select the one that shows the best fit to the data. 

# Geostatistical models in R

The *gstat* package contains functions related to geostatistics. For this example, we will use the `oxford` dataset from this package, which contains measurements of physical and chemical properties for 126 soil samples from a site, along with their coordinates `XCOORD` and `YCOORD`.

```{r, message = FALSE, warning = FALSE}
library(gstat)

data(oxford)
str(oxford)
```

Suppose that we want to model the magnesium concentration (`MG1`), represented as a function of the spatial position in the following graph.

```{r}
library(ggplot2)
ggplot(oxford, aes(x = YCOORD, y = XCOORD, size = MG1)) +
    geom_point() +
    coord_fixed()
```

Note that the $x$ and $y$ axes have been inverted to save space. The `coord_fixed()` function of *ggplot2* ensures that the scale is the same on both axes, which is useful for representing spatial data.

We can immediately see that these measurements were taken on a 100 m grid. It seems that the magnesium concentration is spatially correlated, although it may be a correlation induced by another variable. In particular, we know that the concentration of magnesium is negatively related to the soil pH (`PH1`).

```{r}
ggplot(oxford, aes(x = PH1, y = MG1)) +
    geom_point()
```

The `variogram` function of *gstat* is used to estimate a variogram from empirical data. Here is the result obtained for the variable `MG1`.

```{r}
var_mg <- variogram(MG1 ~ 1, locations = ~ XCOORD + YCOORD, data = oxford)
var_mg
```

The formula `MG1 ~ 1` indicates that no linear predictor is included in this model, while the argument `locations` indicates which variables in the data frame correspond to the spatial coordinates. 

In the resulting table, `gamma` is the value of the variogram for the distance class centered on `dist`, while `np` is the number of pairs of points in that class. Here, since the points are located on a grid, we obtain regular distance classes (e.g.: 100 m for neighboring points on the grid, 141 m for diagonal neighbors, etc.).

Here, we limit ourselves to the estimation of isotropic variograms, i.e. the variogram depends only on the distance between the two points and not on the direction. Although we do not have time to see it today, it is possible with *gstat* to estimate the variogram separately in different directions.

We can illustrate the variogram with `plot`.

```{r}
plot(var_mg, col = "black")
```

If we want to estimate the residual spatial correlation of `MG1` after including the effect of `PH1`, we can add that predictor to the formula.

```{r}
var_mg <- variogram(MG1 ~ PH1, locations = ~ XCOORD + YCOORD, data = oxford)
plot(var_mg, col = "black")
```

Including the effect of pH, the range of the spatial correlation seems to decrease, while the plateau is reached around 300 m. It even seems that the variogram decreases beyond 400 m. In general, we assume that the variance between two points does not decrease with distance, unless there is a periodic spatial pattern.

The function `fit.variogram` accepts as arguments a variogram estimated from the data, as well as a theoretical model described in a `vgm` function, and then estimates the parameters of that model according to the data. The fitting is done by the method of least squares.

For example, `vgm("Exp")` means we want to fit an exponential model.

```{r}
vfit <- fit.variogram(var_mg, vgm("Exp"))
vfit
```

There is no nugget effect, because `psill = 0` for the `Nug` (nugget) part of the model. The exponential part has a sill at 1951 and a range of 95 m.

To compare different models, a vector of model names can be given to `vgm`. In the following example, we include the exponential, gaussian ("Gau") and spherical ("Sph") models.

```{r, warning = FALSE, message = FALSE}
vfit <- fit.variogram(var_mg, vgm(c("Exp", "Gau", "Sph")))
vfit
```

The function gives us the result of the model with the best fit (lowest sum of squared deviations), which here is the same exponential model.

Finally, we can superimpose the theoretical model and the empirical variogram on the same graph.

```{r}
plot(var_mg, vfit, col = "black")
```

## Regression with spatial correlation

We have seen above that the *gstat* package allows us to estimate the variogram of the residuals of a linear model. In our example, the magnesium concentration was modeled as a function of pH, with spatially correlated residuals.

Another tool to fit this same type of model is the `gls` function of the *nlme* package, which is included with the installation of R. 

This function applies the *generalized least squares* method to fit linear regression models when the residuals are not independent or when the residual variance is not the same for all observations. Since the estimates of the coefficients depend on the estimated correlations between the residuals and the residuals themselves depend on the coefficients, the model is fitted by an iterative algorithm:

1. A classical linear regression model (without correlation) is fitted to obtain residuals.

2. The spatial correlation model (variogram) is fitted with those residuals.

3. The regression coefficients are re-estimated, now taking into account the correlations.

Steps 2 and 3 are repeated until the estimates are stable at a desired precision.

Here is the application of this method to the same model for the magnesium concentration in the `oxford` dataset. In the `correlation` argument of `gls`, we specify an exponential correlation model as a function of our spatial coordinates and we include a possible nugget effect.

In addition to the exponential correlation `corExp`, the `gls` function can also estimate a Gaussian (`corGaus`) or spherical (`corSpher`) model.

```{r}
library(nlme)
gls_mg <- gls(MG1 ~ PH1, oxford, 
              correlation = corExp(form = ~ XCOORD + YCOORD, nugget = TRUE))
summary(gls_mg)
```

To compare this result with the adjusted variogram above, the parameters given by `gls` must be transformed. The range has the same meaning in both cases and corresponds to 478 m for the result of `gls`. The global variance of the residuals is the square of `Residual standard error`. The nugget effect here (0.294) is expressed as a fraction of that variance. Finally, to obtain the partial sill of the exponential part, the nugget effect must be subtracted from the total variance. 

After performing these calculations, we can give these parameters to the `vgm` function of *gstat* to superimpose this variogram estimated by `gls` on our variogram of the residuals of the classical linear model.

```{r}
gls_range <- 478
gls_var <- 53.823^2
gls_nugget <- 0.294 * gls_var
gls_psill <- gls_var - gls_nugget

gls_vgm <- vgm("Exp", psill = gls_psill, range = gls_range, nugget = gls_nugget)

plot(var_mg, gls_vgm, col = "black", ylim = c(0, 4000))
```

Does the model fit the data less well here? In fact, this empirical variogram represented by the points was obtained from the residuals of the linear model ignoring the spatial correlation, so it is a biased estimate of the actual spatial correlations. The method is still adequate to quickly check if spatial correlations are present. However, to simultaneously fit the regression coefficients and the spatial correlation parameters, the generalized least squares (GLS) approach is preferable and will produce more accurate estimates.

Finally, note that the result of the `gls` model also gives the AIC, which we can use to compare the fit of different models (with different predictors or different forms of spatial correlation).

## Exercise

The [bryo_belg.csv](data/bryo_belg.csv) dataset is adapted from the data of this study: 

> Neyens, T., Diggle, P.J., Faes, C., Beenaerts, N., Artois, T. et Giorgi, E. (2019) Mapping species richness using opportunistic samples: a case study on ground-floor bryophyte species richness in the Belgian province of Limburg. *Scientific Reports* 9, 19122. https://doi.org/10.1038/s41598-019-55593-x

This data frame shows the specific richness of ground bryophytes (*richness*) for different sampling points in the Belgian province of Limburg, with their position *(x, y)* in km, in addition to information on the proportion of forest (*forest*) and wetlands (*wetland*) in a 1 km^2$ cell containing the sampling point.

```{r}
bryo_belg <- read.csv("data/bryo_belg.csv")
head(bryo_belg)
```

For this exercise, we will use the square root of the specific richness as the response variable. The square root transformation often allows to homogenize the variance of the count data in order to apply a linear regression. 

a) Fit a linear model of the transformed species richness to the proportion of forest and wetlands, without taking into account spatial correlations. What is the effect of the two predictors in this model?

b) Calculate the empirical variogram of the model residuals in (a). Does there appear to be a spatial correlation between the points?

*Note*: The `cutoff` argument to the `variogram` function specifies the maximum distance at which the variogram is calculated. You can manually adjust this value to get a good view of the sill.

c) Re-fit the linear model in (a) with the `gls` function in the *nlme* package, trying different types of spatial correlations (exponential, Gaussian, spherical). Compare the models (including the one without spatial correlation) with the AIC.

d) What is the effect of the proportion of forests and wetlands according to the model in (c)? Explain the differences between the conclusions of this model and the model in (a). 


# Kriging

As mentioned before, a common application of geostatistical models is to predict the value of the response variable at unsampled locations, a form of spatial interpolation called kriging (pronounced with a hard "g").

There are three basic types of kriging based on the assumptions made about the response variable:

- Ordinary kriging: Stationary variable with an unknown mean.

- Simple kriging: Stationary variable with a known mean.

- Universal kriging: Variable with a trend given by a linear or non-linear model.

For all kriging methods, the predictions at a new point are a weighted mean of the values at known points. These weights are chosen so that kriging provides the best linear unbiased prediction of the response variable, if the model assumptions (in particular the variogram) are correct. That is, among all possible unbiased predictions, the weights are chosen to give the minimum mean square error. Kriging also provides an estimate of the uncertainty of each prediction.

While we will not present the detailed kriging equations here, the weights depend on both the correlations (estimated by the variogram) between the sampled points and the new point, as well of the correlations between the sampled points themselves. In other words, sampled points near the new point are given more weight, but isolated sampled points are also given more weight, because sample points close to each other provide redundant information.

Kriging is an interpolation method, so the prediction at a sampled point will always be equal to the measured value (the measurement is supposed to have no error, just spatial variation). However, in the presence of a nugget effect, any small displacement from the sampled location will show variability according to the nugget.

In the example below, we generate a new dataset composed of randomly-generated (x, y) coordinates within the study area as well as randomly-generated pH values based on the `oxford` data. We then apply the function `krige` to predict the magnesium values at these new points. Note that we specify the variogram derived from the GLS results in the `model` argument to `krige`.

```{r}
set.seed(14)
new_points <- data.frame(
    XCOORD = runif(100, min(oxford$XCOORD), max(oxford$XCOORD)),
    YCOORD = runif(100, min(oxford$YCOORD), max(oxford$YCOORD)),
    PH1 = rnorm(100, mean(oxford$PH1), sd(oxford$PH1))
)

pred <- krige(MG1 ~ PH1, locations = ~ XCOORD + YCOORD, data = oxford,
              newdata = new_points, model = gls_vgm)
head(pred)
```

The result of `krige` includes the new point coordinates, the prediction of the variable `var1.pred` along with its estimated variance `var1.var`. In the graph below, we show the mean MG1 predictions from kriging (triangles) along with the measurements (circles).

```{r}
pred$MG1 <- pred$var1.pred

ggplot(oxford, aes(x = YCOORD, y = XCOORD, color = MG1)) +
    geom_point() +
    geom_point(data = pred, shape = 17, size = 2) +
    coord_fixed()
```

The estimated mean and variance from kriging can be used to simulate possible values of the variable at each new point, conditional on the sampled values. In the example below, we performed 4 conditional simulations by adding the argument `nsim = 4` to the same `krige` instruction.

```{r}
sim_mg <- krige(MG1 ~ PH1, locations = ~ XCOORD + YCOORD, data = oxford,
                newdata = new_points, model = gls_vgm, nsim = 4)
head(sim_mg)
```

```{r, message = FALSE, warning = FALSE, fig.dim = c(10, 5)}
library(tidyr)
sim_mg <- pivot_longer(sim_mg, cols = c(sim1, sim2, sim3, sim4), 
                       names_to = "sim", values_to = "MG1")
ggplot(sim_mg, aes(x = YCOORD, y = XCOORD, color = MG1)) +
    geom_point() +
    coord_fixed() +
    facet_wrap(~ sim)
```

# Solutions

```{r}
bryo_lm <- lm(sqrt(richness) ~ forest + wetland, data = bryo_belg)
summary(bryo_lm)
```

The proportion of forest has a significant positive effect and the proportion of wetlands has a significant negative effect on bryophyte richness.

```{r}
plot(variogram(sqrt(richness) ~ forest + wetland, locations = ~ x + y,
               data = bryo_belg, cutoff = 50), col = "black")
```

The variogram is increasing from 0 to at least 40 km, so there appears to be spatial correlations in the model residuals.

```{r}
bryo_exp <- gls(sqrt(richness) ~ forest + wetland, data = bryo_belg,
                correlation = corExp(form = ~ x + y, nugget = TRUE))
bryo_gaus <- gls(sqrt(richness) ~ forest + wetland, data = bryo_belg,
                correlation = corGaus(form = ~ x + y, nugget = TRUE))
bryo_spher <- gls(sqrt(richness) ~ forest + wetland, data = bryo_belg,
                  correlation = corSpher(form = ~ x + y, nugget = TRUE))
```

```{r}
AIC(bryo_lm)
AIC(bryo_exp)
AIC(bryo_gaus)
AIC(bryo_spher)
```

The spherical model has the smallest AIC.

```{r}
summary(bryo_spher)
```

Both effects are less important in magnitude and the effect of wetlands is not significant anymore. As is the case for other types of non-independent residuals, the "effective sample size" here is less than the number of points, since points close to each other provide redundant information. Therefore, the relationship between predictors and response is less clear than given by the model assuming all these points were independent.

Note that the results for all three `gls` models are quite similar, so the choice to include spatial correlations was more important than the exact shape assumed for the variogram.