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?
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
Moran.I function in "ape" library yielding unexpected values
3 messages · Roger Bivand, Steve Pickering
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
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