Classification Trees and Spatial Autocorrelation

[This article was first published on theBioBucket*, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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.


# 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)

# 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")

# 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))

To leave a comment for the author, please follow the link and comment on their blog: theBioBucket*. offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)