Spatial autocorrelation in GAM model residuals in R

Question

I have difficulties with calculating Morans I values and a correlogram. I did a GAM analysis (script below) and I don't know how to find out if there's spatial autocorrelation in model residuals. Maybe you could help with the next steps?

formula <- as.formula(paste('Asustus ~ s(Elupaiga.pindala, k=3) + 
           s(Haava.tagavara, k=3) + 
           s(Lehtpuu.tagavara, k=3) + s(Kuuse.tagavara, k=3) + 
           s(Männi.tagavara, k=3) + 
           s(Kase.tagavara, k=3) + s(Sanglepa.tagavara, k=3) + 
           s(Hall.lepa.tagavara) + s(Aasta, bs=\'re\')'))
model <- gam(formula, data=table, family=binomial())

Edti: I tried different approaches which were recommended but I always got the error at some point. Something to do with the weight dimensions. I tried to solve this but I couldn't get it right.

Answer 1

It might be a good idea to start with a simpler model in order to find potential spatial autocorrelation in your data. The R package 'spdep' provides some spatial econometric models like the CAR or the SAR model. Hence you could extract the residuals from your gam model by applying 'residuals.gam {mgcv}' and then try a spatially autocorrelated model of the form $$Y=c+\rho W Y + u $$ where $Y$ represents your residuals from the gam estimation and $W$ gives a matrix with spatial weights. Then do a simple significance test on $\rho$ in order to find out about spatial autocorrelation. However you could also apply some measure as Morans $I$ which is also implemented in 'spdep' to the residuals. As your predictor seems to be binary it might be also helpful to use the raw residuals. Furthermore you could search for a spatial logit or probit which surely exists for R (see for example https://journal.r-project.org/archive/2013-1/wilhelm-matos.pdf).

Some R code for this idea:

#load the needed packages 
library(Matrix)
library(mgcv)
library(spdep)

# generate a response variable y where the first 100 observations are more likely to be one than the following 100 observations
y<-c(rbinom(100,1,0.8),rbinom(100,1,0.1))

# generate some random covariate
x<-rnorm(200,0,1)

# generate a spatial weigh matrix for the first 100 observations with a lot of neighbours
W1<-matrix(rbinom(10000,1,0.9),ncol=100)

# generate a spatial weigh matrix for the first 100 observations with not very much neighbours
W2<-matrix(rbinom(10000,1,0.1),ncol=100)

# put them together
W<-bdiag(W1,W2)

# do a gam
mod<-gam(y~x,family=binomial())

# extract the residuals
Y<-residuals(mod)

# do a CAR model
lag<-lagsarlm(Y~1,listw = mat2listw(W/rowSums(W)),zero.policy = T)

# check the output of the car model in order to see how big rho is (should be around 0.8-0.9 and highly significant)
summary(lag)

If you are equipped with longitude and latitude in your data you could include it in your gam model via a smooth term interaction in order to control for the spatial structure. (See here Why does including latitude and longitude in a GAM account for spatial autocorrelation?)