famuvie / breedR

Statistical methods for forest genetic resources analysts
http://famuvie.github.io/breedR/
GNU General Public License v3.0
31 stars 24 forks source link

Provide automatically a spatial autocorrelation test on residuals #33

Open famuvie opened 9 years ago

famuvie commented 9 years ago

e.g. Moran's I On behalf of Bruno FADY, PhD, Directeur de recherches, INRA.

famuvie commented 9 years ago

There are so many different tests and flavours for testing spatial autocorrelation, that initially I had made the decision to give the user easy access to the spatial coordinates and residual values and let him apply the method(s) of choice.

This could be something along the lines of

library(breedR)
library(spdep)   # neighbours, distances and Moran's test
res <- breedR:::breedR.result()

## distance between two diagonal neighbours
coord <- coordinates(res)
coordinates(coord) <- c('x', 'y') # cast to SpatialPointsDataFrame
gridded(coord) <- TRUE            # cast to SpatialPixelsDataFrame
dst <- sqrt(sum((slot(slot(coord, "grid"), "cellsize"))**2))

## list of neighbours for each individual
res_nb <- dnearneigh(coordinates(coord), 0, dst)

## Moran's I test
moran.test(residuals(res), nb2listw(res_nb))

Note that it's only two lines of relevant code (apart from loading spdep, and computing the diagonal distance which can be provided directly).

However, I recognize that it doesn't hurt to provide the result of some test in the summary. The interested user can always improve upon that.

For more details on this topic, I recommend reading Chapter 9, and in particular Section 9.4 of Bivand, R. S., E. Pebesma, and V. Gómez-Rubio (2013). Applied Spatial Data Analysis with R (Second ed.). Use R! Springer.

famuvie commented 9 years ago

By the way, I find more interesting and useful the approach described in:

Li, H., C. A. Calder, and N. Cressie (2007). Beyond moran's I: Testing for spatial dependence based on the spatial autoregressive model. Geographical Analysis 39 (4), 357-375. DOI: 10.1111/j.1538-4632.2007.00708.x

Turns out that spdep has it implemented. It can go like this:

## APLE value
aple(residuals(res), nb2listw(res_nb, style="W"))

## APLE test
(boot_res <- aple.mc(residuals(res), nb2listw(res_nb, style="W"), nsim=500))
plot(boot_res)
boot::boot.ci(boot_res, type = 'norm')

Note that there also are plotting methods for visual assessment of both Moran's I and APLE.