Account for spatial autocorrelation with a binomial regression model

I am using a binomial regression model for presence/absence, with 20 independent variables to test. The data has x and y coordinates and I would like to understand how can I take into account the spatial autocorrelation.

I already studied the correlation between the variables and run the same model for 1000 different samples (I have a big dataset that allows me to do this) to understand the distribution of each parameter and check for variables that might be introducing problems in my model.

glm_model <- glm(PA ~ Var1 + Var2 + Var3 + Var4 + Var5,family=binomial(link=logit))

However I believe I also need to account for spatial autocorrelation. I saw that there is a package that might help me (spdep), however I am not sure I completely understand if I can use my model or not. My question is what are my options ?

If you are happy to assume your binomial responses are coming from a spatially correlated gaussian random field via a logit link, and your non-spatial covariates have the usual log-linear form, then stuff it all into geoRglm:

http://cran.r-project.org/web/packages/geoRglm/vignettes/geoRglmintro.pdf

and once you've got your MCMC all tuned, out pops the parameter estimates.

• +1 This package also has non-Bayesian solutions if you want to do some classical variography yourself without specifying priors. – whuber Jul 1 '13 at 15:50

It sounds to me like x, y are potential independent variables in your model. The issue is that the binomial regression model you mention assumes that the independent variables are not correlated. Some people add interaction terms to their model to deal with this, but a lot of model interpretability is lost when you do this.

You have several options for your classification problem. You could use for example k-means clustering. You can find a nice cheat sheet for classification methods in R here.

I found this tutorial for the spdep package to be quite useful.

In the end, you would want to create spatial weights like this:

us.nb4 <- knearneigh(coordinates(data[,1:2]), k=4)
us.nb4 <- knn2nb(us.nb4)
us.nb4 <- make.sym.nb(us.nb4)
us.wt4 <- nb2listw(us.nb4, style="W")

... which you could use into a Moran eigenvector filtering function to remove spatial autocorrelation from the residuals of (generalised) linear models (if you see the tutorial it will make sense).