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