---
title: "Demonstrate GLS regression"
author: "D G Rossiter"
date: "28-December-2017"
output: html_document
---

# Load data

Mercer & Hall uniformity trial: 500 plots with measured wheat grain and straw yield, arranged in a grid.

Reference: Mercer, W. B., & Hall, A. D. (1911). The experimental error of field trials. Journal of Agricultural Science (Cambridge), 4, 107–132.

```{r}
mhw <- read.csv("mhw.csv")
summary(mhw)
```

Four variables: row & column in the square field plot; grain and straw yield in pounds per plot.

Build a model for straw yield as a function of grain yield.

# OLS solution in matrix form

Here the model is $\mathbf{y} = \mathbf{X} \mathbf{\beta} + \mathbf{\varepsilon}$.

```{r}
model.straw.grain <- lm(straw ~ grain, data=mhw)
coefficients(model.straw.grain)
```

The solution is $\mathbf{\widehat{\beta}}_{\mathrm{OLS}} = (\mathbf{X}^T \mathbf{X})^{-1}  \mathbf{X}^T \mathbf{y}$.

Let's see how this is computed.

Design matrix *X*:

```{r}
head(X <-model.matrix(model.straw.grain))
```

Quadratic form (variance-covariance matrix):

```{r}
(X2 <- t(X)%*%X)
```

The inverse of the quadratic form:

```{r}
(X2.inv <- solve(X2))
```

The $y$ vector multiplied by the transpose of the design matrix:

```{r}
str(y <- mhw$straw)
(X2y <- t(X)%*%y)
```

Now we can solve by OLS:

```{r}
(beta.ols <- X2.inv%*%X2y)
```

This is the same as computed directly by `lm`.

# Convert to a spatial object

Compute the  plot length and width in meters, then make a grid covering the centre of each plot.

```{r}
ha2ac <- 1/2.471054  # hectares per acre
ft2m <- .3048      # metres per foot
(field.area <- 10000*ha2ac); (plot.len <- sqrt(field.area)/20); (plot.wid <- sqrt(field.area)/25); (nrow <- length(unique(mhw$r))); (ncol <- length(unique(mhw$c)))
sx <- seq(plot.wid/2, plot.wid/2+(ncol-1)*plot.wid, length=ncol)
sy <- seq(plot.len/2+(nrow-1)*plot.len, plot.len/2, length=nrow)
xy <- expand.grid(x=sx, y=sy)
str(xy)
```

Make a `SpatialPointsDataFrame` version of the dataset; use the grid as spatial coordinates.

```{r}
require(sp)
mhw.sp <- mhw
coordinates(mhw.sp) <- xy
summary(mhw.sp)
spplot(mhw.sp, zcol="grain")
spplot(mhw.sp, zcol="straw", col.regions=topo.colors(64))
summary(mhw.sp$sgr <- mhw.sp$straw/mhw.sp$grain)
spplot(mhw.sp, zcol="sgr", col.regions=terrain.colors(64))
```


# Spatial structure of OLS residuals

Estimate the covariance by the variogram:

Convert the `mhw` object to a `SpatialPointsDataFrame`, add the model residuals as a data field.

```{r}
require(gstat)
mhw.sp$model.resid <- residuals(model.straw.grain)
summary(mhw.sp$model.resid)
bubble(mhw.sp, zcol="model.resid", pch=1)
```

There is obvious spatial dependence in the model residuals, and they seem to be well-separated into regions!

Now compute the residual variogram; there are two ways to express this, exactly the same.

```{r}
bbox(mhw.sp)
(vr <- variogram(model.resid ~ 1, loc=mhw.sp))
vr.alt <- variogram(straw ~ grain, loc=mhw.sp)
round(vr$gamma - vr.alt$gamma, 6)
plot(vr, pl=T)
```

There seems to be dependence to about 10 or 12 m.

```{r}
(vr <- variogram(model.resid ~ 1, loc=mhw.sp, cutoff=15, width=2))
plot(vr, pl=T)
```

Model the residual variogram:

```{r}
(vrmf <- fit.variogram(vr, model=vgm(0.2, "Exp", 3, 0.2)))
plot(vr, pl=T, model=vrmf)
```

So we have an exponential covariance structure. In covariance terms, the structural sill, nugget, total sill, and range parameters are:

```{r}
(c.sill <- vrmf[2,"psill"]); (c.nugget <- vrmf[1,"psill"]); (a <- vrmf[2,"range"])
```


# GLS regression

The model for GLS regression is $\mathbf{y} = \mathbf{X} \mathbf{\beta} + \mathbf{\eta}, \; \mathbf{\eta} \sim \mathcal{N}(0, \mathbf{V})$.

The solution is $\mathbf{\widehat{\beta}}_{\mathrm{GLS}} = (\mathbf{X}^T \mathbf{V}^{-1}\mathbf{X})^{-1}  \mathbf{X}^T \mathbf{V}^{-1} \mathbf{y}$.

That is, we need a variance-covariance matrix of the residuals $V$. We have a model for this, based on the covariance model developed in the previous section.

First, a matrix of the distances between points:

```{r}
D <- spDists(mhw.sp)
dim(D)
D[1:5,1:5]
```

Convert these to covariances according to the exponential covariance function $C(h, a) = \sigma^2 \exp (-{|h/a|})$:

```{r}
# range parameter a, structural sill c, separation h
exp.cov <- function(h, c, a) {
  c * exp(-abs(h/a))
}
V <- c.nugget + exp.cov(D, c.sill, a)
V[1:5, 1:5]
```

So as points are further apart, the covariance decreases. On the diagonal it is the total sill.

Now solve using this:

```{r}
solve(t(X)%*%solve(V)%*%X) %*% t(X)%*%solve(V)%*%y
coefficients(model.straw.grain)
```

Notice the large change in slope: it is much lower, because of less influence of the (spatially-clustered) residuals.

See this also as a `nlme` fit:

```{r}
library(nlme)
# nugget proportion
(s <- vrmf[1,"psill"]/sum(vrmf[,"psill"]))
model.gls.straw.grain <-  gls(model=straw ~ grain,
         data=as.data.frame(mhw.sp),
         correlation=corSpher(
                       value=c(vrmf[2,"range"],s),
                       form=~x+y,
                       nugget=T))
summary(model.gls.straw.grain)
```

Somewhat different because of the re-estimation of the covariance structure, by REML.

# Spatial mean

This just a special case of GLS regression, where the $\mathbf{X}$ matrix is just a column vector of $1$'s. The resulting intercept $\mathbf{\beta}$ is the *spatial* mean, i.e., the mean value but corrected for spatial correlation. If extreme high or low values are clustered, the mean will not be so affected by these.

The $V$ matrix has to be re-computed, because it is not based on residuals, rather on original values.

```{r}
v <- variogram(straw ~ 1, loc=mhw.sp)
plot(v, pl=T)
(vmf <- fit.variogram(v, model=vgm(0.4, "Exp", 10, 0.4)))
plot(v, pl=T, model=vmf)
```

Convert to covariance parameters:

```{r}
(c.sill <- vmf[2,"psill"]); (c.nugget <- vmf[1,"psill"]); (a <- vmf[2,"range"])
```

Obviously, longer range and more variance when not using information about grain yield.

Build the variance-covariance matrix:

```{r}
V <- c.nugget + exp.cov(D, c.sill, a)
V[1:5, 1:5]
```

```{r}
X <- model.matrix(lm(straw ~ 1, data=mhw))
head(X)
solve(t(X)%*%solve(V)%*%X) %*% t(X)%*%solve(V)%*%y
mean(mhw$straw)
lm(straw ~ 1, data=mhw)
```


The *spatial* mean in this case slightly lower, for reasons explained above.

See this also as a `nlme` fit:

```{r}
# nugget proportion
(s <- vmf[1,"psill"]/sum(vmf[,"psill"]))
model.gls.straw.1 <-  gls(model=straw ~ 1,
         data=as.data.frame(mhw.sp),
         correlation=corSpher(
                       value=c(vmf[2,"range"],s ),
                       form=~x+y,
                       nugget=T))
summary(model.gls.straw.1)
coefficients(model.gls.straw.1)
```

Somewhat different because of the re-estimation of the covariance structure, by REML; in this case much closer to the non-spatial mean.