Classification Trees and Spatial Autocorrelation

March 25, 2012
By

(This article was first published on theBioBucket*, and kindly contributed to R-bloggers)

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

To leave a comment for the author, please follow the link and comment on his blog: theBioBucket*.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , , ,

Comments are closed.