[R-sig-Geo] how to use the jouncount.test for spatial autocorrelation test

[Roger Bivand]

On Mon, 17 Mar 2014, Archer_xy wrote:

> I want to test a species' presence / absence records for spatial
> autocorrelation. My data contain >130,000 grids in GIS and with about 700
> species' presence records.

This is not clear at all. See also below for a discussion of whether this 
is a sensible thing to do. If you are going to do it, start with one 
species and one spatial object. Run the example in joincount.test() and 
see how the input objects are structured. In the example, after starting R 
do:

library(spdep)
data(oldcol)
ls()
HICRIME <- cut(COL.OLD$CRIME, breaks=c(0,35,80), labels=c("low","high"))
ls()
# this constructs a factor, a categorical variable
table(HICRIME)
str(HICRIME)

# the list of spatial neighbours is part of the oldcol dataset
print(COL.nb)
class(COL.nb)
joincount.test(HICRIME, nb2listw(COL.nb, style="B"))

# nb2listw converts COL.nb to a binary weights object

You must create the spatial neighbours object for your spatial object. 
Your statement "data contain >130,000 grids in GIS" tells us nothing. 
Having 700 species is not the issue - once you have carried our a 
joincount test for one species, you can do the same for each in turn. If 
you mean that you have a grid with >130K cells (eg. 600*230), then you can 
use for example dnearneigh() to create this object with a suitable 
threshold (below 1.5 with 1x1 cells to create 8 neighbours for each cell 
except on the edges).

library(sp)
GT <- GridTopology(c(0.5, 0.5), c(1, 1), c(600, 230))
SG <- SpatialGrid(GT)
library(spdep)
nb <- dnearneigh(coordinates(SG), 0, 1.5)
# this step takes most time, but the object can be re-used
nb
table(card(nb)) # numbers of neighbours
lwb <- nb2listw(nb, style="B") 
set.seed(140317)
pts <- spsample(SG, n=50000, type="random") # generate random pattern
o <- over(pts, SG)
pa <- rep(0L, prod(slot(GT, "cells.dim")))
pa[o] <- 1L
paf <- factor(pa)
joincount.test(paf, lwb)

is an example, but it should be possible to the raster package focal 
approach too. SG could be read into R using rgdal::readGDAL(), but you 
have problems in ensuring that your tabular species data match the 
ordering of the raster cells - consider creating a SpatialGridDataFrame 
object with multiple species, or a stack or brick also with multiple 
species.

BUT! is this sensible? If your grids are "too small" with respect to the 
scale that applies to one species, you will find strong positive 
autocorrelation resulting simply from the size/resolution of the raster 
cells; if too large, negative autocorrelation. In addition, much of the 
"patchiness" is driven by ecological and environmental factors, to which 
the species respond. Consequently, calculating joincount statistics by 
species may actually tell you nothing at all, perhaps unless you have 
addressed these inhomogeneity issues. If your resolution matches the 
scales for all species and the study area is homogeneous in drivers, the 
results may make sense, but not necessarily otherwise.

>
>
> I have read that the normal Moran's I can't deal with this kind of data, but
> that the join count method in package spdep can do it. However, I'm new to R
> and I still can't understand the information and code in the help for
> joincount.mc or joincount.test.
>
>
> My data is like this:
>
> gridnumber  species
>     1                 1
>     2                 0
>     3                 0
>     4                 1
>     …… 
> I know how to read a shp.file into R and I know I must calculate the weight
> of my data, but the following steps with spdep is beyond my ability.
>
> Any advice in code will be very appreciated.
>
>
>
> --
> View this message in context: http://r-sig-geo.2731867.n2.nabble.com/how-to-use-the-jouncount-test-for-spatial-autocorrelation-test-tp7585965.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no