Testing for Spatial Autocorrelation in Areal Data Via Moran's and Geary's Tests: A Toy Example
First we read in the data from the upper left of the 5x5 grids on page 330 of the lecture notes, which are saved in our spatialstats directory as "acex1.dat":
> acex1 <- matrix(scan(".public-html/spatialstats/acex1.dat"),ncol=2,byrow=T)
> acex1
[,1] [,2]
[1,] 1 5
[2,] 2 4
[3,] 3 4
[4,] 4 4
[5,] 5 4
[6,] 6 5
[7,] 7 5
[8,] 8 5
[9,] 9 3
[10,] 10 3
[11,] 11 1
[12,] 12 5
[13,] 13 2
[14,] 14 3
[15,] 15 4
[16,] 16 1
[17,] 17 1
[18,] 18 2
[19,] 19 2
[20,] 20 3
[21,] 21 1
[22,] 22 2
[23,] 23 1
[24,] 24 2
[25,] 25 3
Next we create a spatial-neighbors object for a 5x5 square grid, with adjacent grid cells defined as neighbors (the "rook's" definition of neighbors):
> library(spdep)
> grid5by5nb <- cell2nb(5,5)
> grid5by5nb
Neighbour list object:
Number of regions: 25
Number of nonzero links: 80
Percentage nonzero weights: 12.8
Average number of links: 3.2
> summary.nb(grid5by5nb)
Neighbour list object:
Number of regions: 25
Number of nonzero links: 80
Percentage nonzero weights: 12.8
Average number of links: 3.2
Link number distribution:
2 3 4
4 12 9
4 least connected regions:
1:1 5:1 1:5 5:5 with 2 links
9 most connected regions:
2:2 3:2 4:2 2:3 3:3 4:3 2:4 3:4 4:4 with 4 links
> nb2listw(grid5by5nb)
Characteristics of weights list object:
Neighbour list object:
Number of regions: 25
Number of nonzero links: 80
Percentage nonzero weights: 12.8
Average number of links: 3.2
Weights style: W
Weights constants summary:
n nn S0 S1 S2
W 25 625 25 16.19444 100.6667
Now we are ready to perform Moran's test for the null hypothesis of no spatial autocorrelation. We first do this test using the p-value based on asymptotic normal theory. The alternative hypothesis is that of positive spatial autocorrelation. Furthermore, we use neighbor weights that are "row-standardized", i.e. made to sum to 1 for each datum:
> moran.test(acex1[,2],nb2listw(grid5by5nb))
Moran's I test under randomisation
data: acex1[, 2]
weights: nb2listw(grid5by5nb)
Moran I statistic standard deviate = 3.4384, p-value = 0.0002926
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.48833333 -0.04166667 0.02375926
Next we repeat this, but using weights that are binary, i.e. all 0 or 1, and not forced to sum to 1 for each datum:
> moran.test(acex1[,2],nb2listw(grid5by5nb,style="B"))
Moran's I test under randomisation
data: acex1[, 2]
weights: nb2listw(grid5by5nb, style = "B")
Moran I statistic standard deviate = 3.2769, p-value = 0.0005248
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.45000000 -0.04166667 0.02251249
Now we list two runs of Moran's test using permutation-based p-values, based on 10000 randomizations, and returning to the use of row-standardized weights:
> moran.mc(acex1[,2],nb2listw(grid5by5nb),10000)
Monte-Carlo simulation of Moran's I
data: acex1[, 2]
weights: nb2listw(grid5by5nb)
number of simulations + 1: 10001
statistic = 0.4883, observed rank = 9999, p-value = 2e-04
alternative hypothesis: greater
> moran.mc(acex1[,2],nb2listw(grid5by5nb),10000)
Monte-Carlo simulation of Moran's I
data: acex1[, 2]
weights: nb2listw(grid5by5nb)
number of simulations + 1: 10001
statistic = 0.4883, observed rank = 10000, p-value = 9.999e-05
alternative hypothesis: greater
Finally, we perform Geary's test, using a p-value corresponding to the asymptotic standard normal distribution.
> geary.test(acex1[,2],nb2listw(grid5by5nb))
Geary's C test under randomisation
data: acex1[, 2]
weights: nb2listw(grid5by5nb)
Geary C statistic standard deviate = 3.4597, p-value = 0.0002704
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic Expectation Variance
0.4856000 1.0000000 0.0221069