GLS regression for raster r

Question

there are a number of posts (1, 2) showing how to conduct OLS regression on raster stacks in R - but I can't find anything about using GLS regression with a specified spatial correlation structure.

Some code below

library(nlme)
library(raster)

# example raster
ras <- brick(nrow = 72, ncol = 144, nl = 5)
values(ras) <- runif(ncell(ras) * nlayers(ras), -1, 1)

time <- ts(seq(1:5))
ras <- setZ(ras, time, "Time")
ras

# define the function to calculate/extract lm slope
calculateSlope <- function(x) {if (is.na(x[1])) {NA} else{ lm(x ~ time)$coefficients[2] }}

# lm slope
slopeRaster <- calc(ras, fun = calculateSlope)
slopeRaster

# define the function to calculate/extract gls slope
calculateGLSSlope <- function(x) {if (is.na(x[1])) {NA} else{ nlme::gls(x ~ time)$coefficients[2] }}

# GLS slope
slopeRasterGLS <- calc(ras, fun = calculateGLSSlope)
slopeRasterGLS

# use all.equal to allow for small differences (sqrt(.Machine$double.eps))
# gls and ols are equal when no correlation structre is defined
all.equal(values(slopeRaster), values(slopeRasterGLS))
[1] TRUE

If I then re-write the slope function to something like

calculateGLSSlope <- function(x) {if (is.na(x[1])) {NA} else{ nlme::gls(x ~ time, correlation = corSpher(form = ~ time)))$coefficients[2] }}

but I'm not 100 % sure if this accounts for the X and Y coords of the data, only the correlation between time. How is it possible to specify the x~y coords of each of the raster cells when calc works on a pixel-by-pixel basis? I thought something along the lines of

calculateSlope <- function(j) {if (is.na(j[1])) {NA} else{
  coords <- xyFromCell(ras, j)
  nlme::gls(j ~ Year, correlation = corSpher(form = ~ coords$x + coords$y + time))$coefficients[2] }}

But in this case j is represented by the cell-value and there is no information about the cell location, so xyFromCell fails

Edit 11/05/2017

Following on from @geoCompR answer below, I've attempted to follow their steps to get the lat/long coordinates in a model. The additions to the code above are below.

# create rasters for latitude and longitude
rLat <- raster(ncols = ncol(ras), nrows = nrow(ras))
values(rLat) <- values(init(ras, "y"))
rLon <- raster(ncols = ncol(ras), nrows = nrow(ras))
values(rLon) <- values(init(ras, "x"))

# Stack the original raster with the lon and lat rasters
newStack <- stack(ras, rLon, rLat)
newStack

# define the function to calculate/extract lm slope with coords as var
calculateSlopeWCoords <- function(x) {if (is.na(x[1])) {NA} else{
# layers x[1:5] is the vector of dependent variables, time is the time vector from above
# x[6] is the longitude, x[7] is the latitude
lm(x[1:5] ~ time + x[6] + x[7])$coefficients[2] }}
slopeRaster2 <- calc(newStack, fun = calculateSlopeWCoords)

However, when attempting to run slopeRaster2 I receive the following error

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
  cannot use this function

How can I include the lat and long variables in a function that can be passed to calc?

If I try to create a raster stack where I have 5 layers of dependent variables, 5 layers for time, and then 5 each for latitude and longitude; I could then write a function like this

slopeFun = function(x) {
  # In case of any NA return four NA values for coefficient estimates
  # Four variables = four coefficients (1 intercept, 3 indep)
  if (any(is.na(x))) return(c(NA, NA, NA, NA))
  model <- lm(x[1:5] ~ x[6:10] + x[11:15] + x[16:20], na.action = na.omit)
  coefficients(model)}

This can be passed to calc like so: slopeRas <- calc(rasStack, fun = slopeFun), which returns a raster stack of model coefficients, but the coefficients for the Latitude and Longitude layers all come back as NA as each of the layers for latitude, and each of the layers for longitude, are perfectly correlated with themselves and therefore can't be modelled with lm. The code below shows this approach.

library(raster)
library(nlme)
# Generate an example raster brick
set.seed(1)
ras <- raster(res = 2.5)
brRas <- brick(ncols = ncol(ras), nrows = nrow(ras), nl = 5)
values(brRas) <- rnorm(ncell(brRas) * nlayers(brRas), 0, 1)
brRas
# generate a brick to represent time
timeRas <- brick(ncols = ncol(ras), nrows = nrow(ras), nl = 5)
for (r in 1:nlayers(timeRas)){
  values(timeRas[[r]]) <- seq(1:5)[r]
}
timeRas
time <- ts(seq(1:5))
# Generate Lat/Long raster bricks
lat <- brick(ncols = ncol(ras), nrows = nrow(ras), nl = 5)
values(lat) <- values(init(brRas, "y"))
long <- brick(ncols = ncol(ras), nrows = nrow(ras), nl = 5)
values(long) <- values(init(brRas, "x"))
# Set the layer names so they are easily identifiable
names(brRas) <- paste0("Layer", seq(1:5))
names(timeRas) <- paste0("Time", seq(1:5))
names(lat) <- paste0("Lat", seq(1:5))
names(long)<- paste0("Long", seq(1:5))
# Stack it all together and check
rasStack <- stack(brRas, timeRas, long, lat)
rasStack
names(rasStack)
# Define the function for the coefficients
slopeFun <- function(x) {
      # In case of any NA return four NA values for coefficient estimates
      # Four variables in the model = four coefficients (1 intercept, 3 indep)
      if (any(is.na(x))) return(c(NA, NA, NA, NA))
      model <- lm(x[1:5] ~ x[6:10] + x[11:15] + x[16:20], na.action = na.omit)
      coefficients(model)}
# Pass the function to calc
slopeRas <- calc(rasStack, fun = slopeFun)
slopeRas

class       : RasterBrick 
dimensions  : 72, 144, 10368, 4  (nrow, ncol, ncell, nlayers)
resolution  : 2.5, 2.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : X.Intercept.,   x.6.10.,  x.11.15.,  x.16.20. 
min values  :    -4.236375, -1.066759,        NA,        NA 
max values  :     4.112934,  1.228611,        NA,        NA 

The only other thing I can think of doing, is creating a dataframe of my raster stack containing xy centroids, the value of the dependent variable, the z-index of the raster (time), and the CellID, and then fitting the models row-wise, and then summarising the results by CellID to avoid the perfect collinearity issue.

Answer 1

Here's a way to include Lat/Long in your equation. You'll still need to check for spatial autocorrelation in your residuals to ensure that, in fact, this global fix is enough to ameliorate the effects of S.A.!

You can convert it to a matrix of cell coordinates, create two new rasters from the old one, and replace their values with the point coordinates. May not be great for a huge raster... for that you'll probably have to do some math and convert the bounding box/nrow/ncell into something to speed it up. As far as the first method goes, this worked lightning fast on a raster containing ~1 million cells:

require(raster)
require(rgdal)
r<-raster("C:/your/raster.tif")

x<-xyFromCell(r,1:ncell(r)) # get cell coordinates as matrix
lon<-x[,1] # get all longitude as vector
lat<-x[,2] # get all latitude as vector

rLon<-r # make a copy of the original raster (easier than creating a blank one)
rLat<-r # make a copy of the original raster (easier than creating a blank one)

values(rLon)<-lon # replace values of the raster
values(rLat)<-lat # replace values of the raster

This will leave you with two rasters (rLon and rLat) that you can then add to your stack and use as independent variables in your model. Happy modeling, glad to see people doing raster modeling in R (basically all I do all day...)! ... keep in mind: You don't have to stick with lm() for this. Any model can be applied to the stack - this includes SVM, Neural Nets, or Random Forest (my favorite... it plays well with spatial data, has no assumptions of normality, eliminates the need to check for multicollinearity), just to name a few!