Classification Trees and Spatial Autocorrelation

I'm currently trying to model species presence / absence data (N = 523) that were collected over a geographic area and are possibly spatially autocorrelated. Samples come from preferential sites (sea level > 1200 m, obligatory presence of permanent waterbodies, etc). My main goal is to infere on environmental factors determining the occurrence rate of several (amphibian) species and to rule out spatial autocorrelation.

A fellow R-sig-eco list member assisted me in working out a method for dealing with spatial autocorrelation in classification trees (conditional inference trees in R-package party).

The whole idea is actually plain simple: compute the classification tree, calculate residuals and use it for a Mantel-test and Mantel correlograms.
The Mantel correlograms test differrences in dissimilarities of
the residuals across several spatial distances and thus enable you to detect lag-distances where possible spatial autocorrelation vanishes.

I did the mantel tests and correlograms with the package ecodist.

The below toy example is obviously not spatially autocorrelated - as is shown by the Mantel tests and correlograms. If I'd encounter autocorrelation with real-world data, I'd try to use subsamples of the data avoiding resampling within the lag-distance..

For an example with this scheme in action see this paper by Bálint Czúcz et al.

library(party)

# for simplicity only one env-factor:
env <- data.frame(F1 = as.factor(rep(LETTERS[1:4], each = 40)))

n <- nrow(env)

# probs for simulated sampling:
p <- seq(0, 1, 1/n)

# simulated binary response:
binresp <- numeric()

for(i in 1:n) {
   binresp[i] <- sample(0:1, 1, replace = T, c(p[i], 1-p[i]))
   }

binresp <- as.factor(binresp)

ct <- ctree(binresp ~ env$F1)
plot(ct)

# predicted nodes:
pred <- names(predict(ct))

# predicted probabilties:
ct_resp <- treeresponse(ct)

# pull second elements of list items:
ct_resp <- sapply(ct_resp, "[[", 2)

# convert response from factor back to numeric:
binresp <- as.numeric(ifelse(binresp == 1, 1, 0))

# calculate residuals:
Residuals <- ct_resp - binresp
cbind(ct_resp, binresp, Residuals)


# simulate XY-coords:
Y <- jitter(sample(39:42, n, replace = T))
X <- jitter(sample(10:12, n, replace = T))
XY <- cbind(X, Y)

# install.packages("ecodist")
library(ecodist)

# make distance matrices:
d_XY <- lower(dist(XY))
d_Resid <- lower(dist(Residuals))

# check distribution of residuals:
table(d_Resid); hist(d_Resid)

# mantel tests (with Pearson's and Spearman's Corr. Coeff.):
mantel(d_Resid ~ d_XY, nperm = 10000)
mantel(d_Resid ~ d_XY, nperm = 10000, mrank = T)

# mantel correlograms:
plot(mgram(d_Resid, d_XY))
plot(mgram(d_Resid, d_XY, mrank = T))