[R-sig-Geo] Moran.I function in "ape" library yielding unexpected values

[Roger Bivand]

On Fri, 5 Dec 2014, Steve Pickering wrote:

> Dear r-sig-geo community,
>
> I'm running the Moran.I function in the "ape" library and getting some
> values which are not what I'm expecting.  Consider the code below, in
> which a simple 8 x 8 grid is set up, and then populated with a
> chessboard pattern.  My understanding is that a chessboard is an
> example of perfect negative spatial autocorrelation, and that it
> should therefore yield a Moran's I value of -1.  The result it yields,
> while negative, is not close to -1:
>
> $observed
> [1] -0.07775248
>
> $expected
> [1] -0.01587302
>
> $sd
> [1] 0.01499786
>
> $p.value
> [1] 3.693094e-05
>
> The code is as follows.  Can anyone tell me what I'm doing wrong?

Your definition of all observations as neighbours of each other is causing 
this result. It isn't wrong, but very often the choice of weights will 
have an effect. Try:

library(spdep)
nb8r <- cell2nb(8, 8, type="rook", torus=FALSE)
nb8q <- cell2nb(8, 8, type="queen", torus=FALSE)
nb8rt <- cell2nb(8, 8, type="rook", torus=TRUE)
nb8qt <- cell2nb(8, 8, type="queen", torus=TRUE)

moran.test(testGrid$chess, nb2listw(nb8r, style="B"),
   alternative="two.sided")
# -1.000000000
moran.test(testGrid$chess, nb2listw(nb8q, style="B"),
   alternative="two.sided")
# -0.066666667
moran.test(testGrid$chess, nb2listw(nb8rt, style="B"),
   alternative="two.sided")
# -1.000000000
moran.test(testGrid$chess, nb2listw(nb8qt, style="B"),
   alternative="two.sided")
# 0.000000000

rdists <- nbdists(nb8r, cbind(testGrid$lng, testGrid$lat))
irdists <- lapply(rdists, function(x) 1/x)
qdists <- nbdists(nb8q, cbind(testGrid$lng, testGrid$lat))
iqdists <- lapply(qdists, function(x) 1/x)

moran.test(testGrid$chess, nb2listw(nb8r, glist=irdists, style="B"),
   alternative="two.sided")
# -1.000000000
moran.test(testGrid$chess, nb2listw(nb8q, glist=iqdists, style="B"),
   alternative="two.sided")
# -0.235545329

nb8all <- dnearneigh(cbind(testGrid$lng, testGrid$lat), 0,
   sqrt(2*(7^2)))
alldist <- nbdists(nb8all, cbind(testGrid$lng, testGrid$lat))
ialldist <- lapply(alldist, function(x) 1/x)
summary(rowSums(testGrid.dists.inv))
summary(sapply(ialldist, sum))
moran.test(testGrid$chess, nb2listw(nb8all, glist=ialldist, style="B"),
   randomisation=TRUE, alternative="two.sided")
# -0.0769782624

So the design of the spatial weights is driving the outcome. In addition, 
the 0/1 nature of the data is not strictly appropriate for Moran's I, 
which could take the residuals of a glm():

mod <- glm(chess ~ 1, data=testGrid, family=binomial)
lm.morantest(mod, nb2listw(nb8rt, style="B"), alternative="two.sided")
# -1.00000000

or join-count tests:

joincount.multi(as.factor(testGrid$chess), nb2listw(nb8r, style="B"))
#      Joincount Expected Variance z-value
# 0:0      0.000   27.556    8.325 -9.5503
# 1:1      0.000   27.556    8.325 -9.5503
# 1:0    112.000   56.889   27.205 10.5662
# Jtot   112.000   56.889   27.205 10.5662

joincount.multi(as.factor(testGrid$chess), nb2listw(nb8q, style="B"))
#      Joincount Expected Variance z-value
# 0:0     49.000   51.667   23.616 -0.5487
# 1:1     49.000   51.667   23.616 -0.5487
# 1:0    112.000  106.667   47.796  0.7714
# Jtot   112.000  106.667   47.796  0.7714

joincount.multi(as.factor(testGrid$chess), nb2listw(nb8all, 
glist=ialldist, style="B"))
#      Joincount Expected Variance z-value
# 0:0    148.864  158.719   34.009 -1.6899
# 1:1    148.864  158.719   34.009 -1.6899
# 1:0    347.388  327.678   22.506  4.1547
# Jtot   347.388  327.678   22.506  4.1547

In the rook case, there are no 0:0 or 1:1 joins, but there are in the 
Queen and all-in inverted distance cases.

I'm not sure that this clarifies, but it does account for your result.

Roger

>
> As always, many thanks!
>
> - Steve.
>
>
>
> library(ape)
>
> # first, build an 8 x 8 grid
> # lats and lngs will go from 0 to 7
>
> lng <- rep(seq(0, 7, by=1), 8)
>
> counter = 1
> subCounter = 0
> startNum = 0
> lat = NULL
>
> while (counter < 65) {
>
> if (subCounter == 8) {
> startNum = startNum + 1
> subCounter = 0
> }
>
> lat = c(lat, startNum)
> subCounter = subCounter + 1
> counter = counter + 1
>
> }
>
> # now add the black/ white chessboard pattern
>
> chess <- rep(c(0,1,0,1,0,1,0,1,1,0,1,0,1,0,1,0), 4)
>
> gridNum <- seq(1:64)
>
> testGrid <- data.frame(gridNum, lat, lng, chess)
>
> # simple Euclidean distance measure
>
> testGrid.dists <- as.matrix(dist(cbind(testGrid$lng, testGrid$lat)))
>
> # invert the distances, then replace the infinities with 0
>
> testGrid.dists.inv <- 1 / testGrid.dists
>
> diag(testGrid.dists.inv) <- 0
>
> # now test the chessboard pattern
>
> Moran.I(testGrid$chess, testGrid.dists.inv)
>
>
>
> Steve Pickering
> Assistant Professor
> Graduate School of Law
> Kobe University
> 2-1 Rokkodai-cho
> Nada-ku, Kobe
> 657-8501 Japan
> pickering at penguin.kobe-u.ac.jp
> http://www.stevepickering.net
>
> _______________________________________________
> 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