[R-sig-Geo] Moran.I function in "ape" library yielding unexpected values
Steve Pickering
pickering at penguin.kobe-u.ac.jp
Fri Dec 5 15:51:20 CET 2014
Certainly does clarify, and gives me the results I was expecting with
simulated and real data: thanks, Roger!
- Steve.
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
On 5 December 2014 at 23:03, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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
>
More information about the R-sig-Geo
mailing list