Chapter 15 Measures of spatial autocorrelation

When analysing areal data, it has long been recognised that, if present, spatial autocorrelation changes how we may infer, relative to the default position of independent observations. In the presence of spatial autocorrelation, we can predict the values of observation \(i\) from the values observed at \(j \in N_i\), the set of its proximate neighbours. Early results (Moran 1948; Geary 1954), entered into research practice gradually, for example the social sciences (Duncan, Cuzzort, and Duncan 1961). These results were then collated and extended to yield a set of basic tools of analysis (Cliff and Ord 1973, 1981).

Cliff and Ord (1973) generalised and extended the expression of the spatial weights matrix representation as part of the framework for establishing the distribution theory for join count, Moran’s \(I\) and Geary’s \(C\) statistics. This development of what have become known as global measures, returning a single value of autocorrelation for the total study area, has been supplemented by local measures returning values for each areal unit (Getis and Ord 1992; Anselin 1995).

15.1 Measures and process mis-specification

It is not and has never been the case that Tobler’s first law of geography: “Everything is related to everything else, but near things are more related than distant things” always holds absolutely. This is and has always been an oversimplification, disguising the underlying entitation, support and other mis-specification problems. Are the units of observation appropriate for the scale of the underlying spatial process? Could the spatial patterning of the variable of interest for the chosen entitation be accounted for by another variable?

Tobler (1970) was published in the same special issue of Economic Geography as Olsson (1970), but Olsson does grasp the important point that spatial autocorrelation is not inherent in spatial phenomena, but often is engendered by inappropriate entitation, by omitted variables and/or inappropriate functional form. The key quote from Olsson is on p. 228:

The existence of such autocorrelations makes it tempting to agree with Tobler (1970, 236 [the original refers to the pagination of a conference paper]) that ‘everything is related to everything else, but near things are more related than distant things.’ On the other hand, the fact that the autocorrelations seem to hide systematic specification errors suggests that the elevation of this statement to the status of ‘the first law of geography’ is at best premature. At worst, the statement may represent the spatial variant of the post hoc fallacy, which would mean that coincidence has been mistaken for a causal relation.

The status of the “first law” is very similar to the belief that John Snow induced the cause of cholera as water-borne from a map. It may be a good way of selling GIS, but it is inaccurate; Snow had a strong working hypothesis prior to visiting Soho, and the map was prepared after the Broad street pump was disabled as documentation that the hypothesis held (Brody et al. 2000).

Measures of spatial autocorrelation unfortunately pick up other mis-specifications in the way that we model data (Schabenberger and Gotway 2005; McMillen 2003). For reference, Moran’s \(I\) is given as (Cliff and Ord 1981, 17):

\[ I = \frac{n \sum_{(2)} w_{ij} z_i z_j}{S_0 \sum_{i=1}^{n} z_i^2} \] where \(x_i, i=1, \ldots, n\) are \(n\) observations on the numeric variable of interest, \(z_i = x_i - \bar{x}\), \(\bar{x} = \sum_{i=1}^{n} x_i / n\), \(\sum_{(2)} = \stackrel{\sum_{i=1}^{n} \sum_{j=1}^{n}}{i \neq j}\), \(w_{ij}\) are the spatial weights, and \(S_0 = \sum_{(2)} w_{ij}\). First we test a random variable using the Moran test, here under the normality assumption (argument randomisation=FALSE, default TRUE). Inference is made on the statistic \(Z(I) = \frac{I - E(I)}{\sqrt{\mathrm{Var}(I)}}\), the z-value compared with the Normal distribution for \(E(I)\) and \(\mathrm{Var}(I)\) for the chosen assumptions; this x does not show spatial autocorrelation with these spatial weights:

# Moran I statistic       Expectation          Variance       Std deviate 
#         -0.004772         -0.000401          0.000140         -0.369320 
#           p.value 
#          0.711889

The test however detects quite strong positive spatial autocorrelation when we insert a gentle trend into the data, but omit to include it in the mean model, thus creating a missing variable problem but finding spatial autocorrelation instead:

# Moran I statistic       Expectation          Variance       Std deviate 
#          0.043403         -0.000401          0.000140          3.701491 
#           p.value 
#          0.000214

If we test the residuals of a linear model including the trend, the apparent spatial autocorrelation disappears:

# Observed Moran I      Expectation         Variance      Std deviate 
#        -0.004777        -0.000789         0.000140        -0.337306 
#          p.value 
#         0.735886

A comparison of implementations of measures of spatial autocorrelation shows that a wide range of measures is available in R in a number of packages, chiefly in the spdep package (R. Bivand 2021b), and that differences from other implementations can be attributed to design decisions (Bivand and Wong 2018). The spdep package also includes the only implementations of exact and saddlepoint approximations to global and local Moran’s I for regression residuals (Tiefelsdorf 2002; Bivand, Müller, and Reder 2009).

15.2 Global measures

Global measures consider the average level of spatial autocorrelation across all observations; they can of course be biassed (as most spatial statistics) by edge effects where important spatial process components fall outside the study area.

15.2.1 Join-count tests for categorical data

We will begin by examining join count statistics, where joincount.test() takes a "factor" vector of values fx= and a listw object, and returns a list of htest (hypothesis test) objects defined in the stats package, one htest object for each level of the fx= argument. The observed counts are of neighbours with the same factor levels, known as same-colour joins.

# function (fx, listw, zero.policy = NULL, alternative = "greater", 
#     sampling = "nonfree", spChk = NULL, adjust.n = TRUE) 

The function takes an alternative= argument for hypothesis testing, a sampling= argument showing the basis for the construction of the variance of the measure, where the default "nonfree" choice corresponds to analytical permutation; the spChk= argument is retained for backward compatibility. For reference, the counts of factor levels for the type of municipality or Warsaw borough are:

# .
#          Rural          Urban    Urban/rural Warsaw Borough 
#           1563            303            611             18

Since there are four levels, we re-arrange the list of htest objects to give a matrix of estimated results. The observed same-colour join counts are tabulated with their expectations based on the counts of levels of the input factor, so that few joins would be expected between for example Warsaw boroughs, because there are very few of them. The variance calculation uses the underlying constants of the chosen listw object and the counts of levels of the input factor. The z-value is obtained in the usual way by dividing the difference between the observed and expected join counts by the square root of the variance.

The join count test was subsequently adapted for multi-colour join counts (Upton and Fingleton 1985). The implementation as joincount.mult() in spdep returns a table based on nonfree sampling, and does not report p-values.

#                               Joincount Expected Variance z-value
# Rural:Rural                    3087.000 2793.920 1126.534    8.73
# Urban:Urban                     110.000  104.719   93.299    0.55
# Urban/rural:Urban/rural         656.000  426.526  331.759   12.60
# Warsaw Borough:Warsaw Borough    41.000    0.350    0.347   68.96
# Urban:Rural                     668.000 1083.941  708.209  -15.63
# Urban/rural:Rural              2359.000 2185.769 1267.131    4.87
# Urban/rural:Urban               171.000  423.729  352.190  -13.47
# Warsaw Borough:Rural             12.000   64.393   46.460   -7.69
# Warsaw Borough:Urban              9.000   12.483   11.758   -1.02
# Warsaw Borough:Urban/rural        8.000   25.172   22.354   -3.63
# Jtot                           3227.000 3795.486 1496.398  -14.70

So far, we have used binary weights, so the sum of join counts multiplied by the weight on that join remains integer. If we change to row standardised weights, where the weights almost always fractions of 1, the counts, expectations and variances change, but there are few major changes in the z-values.

Using an inverse distance based listw object does, however, change the z-values markedly, because closer centroids are upweighted relatively strongly:

#                               Joincount Expected Variance z-value
# Rural:Rural                    3.46e+02 3.61e+02 4.93e+01   -2.10
# Urban:Urban                    2.90e+01 1.35e+01 2.23e+00   10.39
# Urban/rural:Urban/rural        4.65e+01 5.51e+01 9.61e+00   -2.79
# Warsaw Borough:Warsaw Borough  1.68e+01 4.53e-02 6.61e-03  206.38
# Urban:Rural                    2.02e+02 1.40e+02 2.36e+01   12.73
# Urban/rural:Rural              2.25e+02 2.83e+02 3.59e+01   -9.59
# Urban/rural:Urban              3.65e+01 5.48e+01 8.86e+00   -6.14
# Warsaw Borough:Rural           5.65e+00 8.33e+00 1.73e+00   -2.04
# Warsaw Borough:Urban           9.18e+00 1.61e+00 2.54e-01   15.01
# Warsaw Borough:Urban/rural     3.27e+00 3.25e+00 5.52e-01    0.02
# Jtot                           4.82e+02 4.91e+02 4.16e+01   -1.38

15.2.2 Moran’s \(I\)

The implementation of Moran’s \(I\) in spdep in the moran.test() function has similar arguments to those of joincount.test(), but sampling= is replaced by randomisation= to indicate the underlying analytical approach used for calculating the variance of the measure. It is also possible to use ranks rather than numerical values (Cliff and Ord 1981, 46). The drop.EI2= argument may be used to reproduce results where the final component of the variance term is omitted as found in some legacy software implementations.

# function (x, listw, randomisation = TRUE, zero.policy = NULL, 
#     alternative = "greater", rank = FALSE, na.action =, 
#     spChk = NULL, adjust.n = TRUE, drop.EI2 = FALSE) 

The default for the randomisation= argument is TRUE, but here we will simply show that the test under normality is the same as a test of least squares residuals with only the intercept used in the mean model; the analysed variable is first round turnout proportion of registered voters in municipalities and Warsaw burroughs in the 2015 Polish presidential election. The spelling of randomisation is that of Cliff and Ord (1973).

# Moran I statistic       Expectation          Variance       Std deviate 
#          0.691434         -0.000401          0.000140         58.461349 
#           p.value 
#          0.000000

The lm.morantest() function also takes a resfun= argument to set the function used to extract the residuals used for testing, and clearly lets us model other salient features of the response variable (Cliff and Ord 1981, 203). To compare with the standard test, we are only using the intercept here and, as can be seen, the results are the same.

# Observed Moran I      Expectation         Variance      Std deviate 
#         0.691434        -0.000401         0.000140        58.461349 
#          p.value 
#         0.000000

The only difference between tests under normality and randomisation is that an extra term is added if the kurtosis of the variable of interest indicates a flatter or more peaked distribution, where the measure used is the classical measure of kurtosis. Under the default randomisation assumption of analytical randomisation, the results are largely unchanged.

# Moran I statistic       Expectation          Variance       Std deviate 
#          0.691434         -0.000401          0.000140         58.459835 
#           p.value 
#          0.000000

From the very beginning in the early 1970s, interest was shown in Monte Carlo tests, also known as Hope-type tests and as permutation bootstrap. By default, returns a "htest" object, but may simply use boot::boot() internally and return a "boot" object when return_boot=TRUE. In addition the number of simulations needs to be given as nsim=; that is the number of times the values of the observations are shuffled at random.

The bootstrap permutation retains the outcomes of each of the random permutations, reporting the observed value of the statistic, here Moran’s \(I\), the difference between this value and the mean of the simulations under randomisation (equivalent to \(E(I)\)), and the standard deviation of the simulations under randomisation.

If we compare the Monte Carlo and analytical variances of \(I\) under randomisation, we typically see few differences, arguably rendering Monte Carlo testing unnecessary.

#    Permutation bootstrap Analytical randomisation 
#                 0.000144                 0.000140

Geary’s global \(C\) is implemented in geary.test() largely following the same argument structure as moran.test(). The Getis-Ord \(G\) test includes extra arguments to accommodate differences between implementations, as Bivand and Wong (2018) found multiple divergences from the original definitions, often to omit no-neighbour observations generated when using distance band neighbours. It is given by (Getis and Ord 1992, 194). For \(G^*\), the \(\sum_{(2)}\) constraint is relaxed by including \(i\) as a neighbour of itself (thereby also removing the no-neighbour problem, because all observations have at least one neighbour).

Finally, the empirical Bayes Moran’s \(I\) takes account of the denominator in assessing spatial autocorrelation in rates data (Assunção and Reis 1999). Until now, we have considered the proportion of valid votes cast in relation to the numbers entitled to vote by spatial entity, but using we can try to accommodate uncertainty in extreme rates in entities with small numbers entitled to vote. There is, however, little impact on the outcome in this case.

Global measures of spatial autocorrelation using spatial weights objects based on graphs of neighbours are, as we have seen, rather blunt tools, which for interpretation depend critically on a reasoned mean model of the variable in question. If the mean model is just the intercept, the global measures will respond to all kinds of mis-specification, not only spatial autocorrelation. The choice of entities for aggregation of data will typically be a key source of mis-specification.

15.3 Local measures

Building on insights from the weaknesses of global measures, local indicators of spatial association began to appear in the first half of the 1990s (Anselin 1995; Getis and Ord 1992, 1996).

In addition, the Moran plot was introduced, plotting the values of the variable of interest against their spatially lagged values, typically using row-standardised weights to make the axes more directly comparable (Anselin 1996). The moran.plot() function also returns an influence measures object used to label observations exerting more than propotional influence on the slope of the line representing global Moran’s \(I\). In figure 15.1, we can see that there are many spatial entities exerting such influence. These pairs of observed and lagged observed values make up in aggregate the global measure, but can also be explored in detail. The quadrants of the Moran plot also show low-low pairs in the lower left quadrant, high-high in the upper right quadrant, and fewer low-high and high-low pairs in the upper left and lower right quadrants.

Moran plot of I round turnout, row standardised weights

Figure 15.1: Moran plot of I round turnout, row standardised weights

If we extract the hat value influence measure from the returned object, Figure 15.2 suggests that some edge entities exert more than proportional influence (perhaps because of row standardisation), as do entities in or near larger urban areas.

Moran plot hat values, row standardised neighbours

Figure 15.2: Moran plot hat values, row standardised neighbours

15.3.1 Local Moran’s \(I_i\)

Bivand and Wong (2018) discuss issues impacting the use of local indicators, such as local Moran’s \(I_i\) and local Getis-Ord \(G_i\). Some issues affect the calculation of the local indicators, others inference from their values. Because \(n\) statistics may be being calculated from the same number of observations, there are multiple comparison problems that need to be addressed. Although the apparent detection of hotspots from values of local indicators has been quite widely adopted, it remains fraught with difficulty because adjustment of the inferential basis to accommodate multiple comparisons is not often chosen, and as in the global case, mis-specification also remains a source of confusion. Further, interpreting local spatial autocorrelation in the presence of global spatial autocorrelation is challenging (Ord and Getis 2001; Tiefelsdorf 2002; Bivand, Müller, and Reder 2009).

# function (x, listw, zero.policy = NULL, na.action =, 
#     conditional = FALSE, alternative = "greater", p.adjust.method = "none", 
#     mlvar = TRUE, spChk = NULL, adjust.x = FALSE) 

The mlvar= and adjust.x= arguments to localmoran() are discussed in Bivand and Wong (2018), and permit matching with other implementations. The p.adjust.method= argument uses an untested speculation implemented in p.adjustSP() that adjustment should only take into account the cardinality of the neighbour set of each observation when adjusting for multiple comparisons; using stats::p.adjust() is preferable.

Taking "two.sided" p-values, we obtain:

The \(I_i\) local indicators when summed and divided by the sum of the spatial weights equal global Moran’s \(I\), showing the possible presence of positive and negative local spatial autocorrelation:

# [1] TRUE

Using stats::p.adjust() to adjust for multiple comparisons, we see that almost 29% of the 2495 local measures have p-values < 0.05 if no adjustment is applied, but only 12% using Bonferroni adjustment to control the familywise error rate, with two other choices shown: fdr is the Benjamini and Hochberg (1995) false discovery rate and BY (Benjamini and Yekutieli 2001), another false discovery rate adjustment:

#       none bonferroni        fdr         BY 
#        715        297        576        424

In the global measure case, bootstrap permutations may be used as an alternative to analytical methods for possible inference, where both the theoretical development of the analytical variance of the measure, and the permutation scheme, shuffle all of the observed values. In the local case, conditional permutation may be used, fixing the value at observation \(i\) and randomly sampling from the remaining \(n-1\) values to find randomised values at neighbours, and is provided as localmoran_perm(), which may use multiple nodes to sample in parallel if provided, and permits the setting of a seed for the random number generator across the compute nodes:

#    user  system elapsed 
#   0.629   1.513   0.695

The outcome is that almost 32% of observations have two sided p-values < 0.05 without multiple comparison adjustment, and under 3% with Bonferroni adjustment.

#       none bonferroni        fdr         BY 
#        797         76        463        161

We can see what is happening by tabulating counts of the standard deviate of local Moran’s \(I\), where the two-sided \(\alpha=0.05\) bounds would be \(0.025\) and \(0.975\), but Bonferroni adjustment is close to \(0.00001\) and \(0.99999\). Without adjustment, almost 800 observations are significant, with Bonferroni adjustment, only almost 70 in the conditional permutation case:

# .
#  (-Inf,-4.26] (-4.26,-3.72] (-3.72,-3.09] (-3.09,-2.33] (-2.33,-1.96] 
#             0             0             1             4             5 
#     (-1.96,0]      (0,1.96]   (1.96,2.33]   (2.33,3.09]   (3.09,3.72] 
#           459          1239           195           316           145 
#   (3.72,4.26]   (4.26, Inf] 
#            55            76
# [1] 797
# [1] 76

In an important clarification, Sauer et al. (2021) show that the comparison of standard deviates for local Moran’s \(I_i\) based on analytical formulae and conditional permutation in Bivand and Wong (2018) was based on a misunderstanding. Sokal, Oden, and Thomson (1998) provide alternative analytical formulae for standard deviates of local Moran’s \(I_i\) based either on total or conditional permutation, but the analytical formulae used in Bivand and Wong (2018), based on earlier practice, only use total permutation, and consequently do not match the simulation conditional permutations. Thanks to a timely pull request, localmoran() now has a conditional= argument using alternative formulae from the appendix of Sokal, Oden, and Thomson (1998):

yielding standard deviates that correspond closely to those from simulation:

#       none bonferroni        fdr         BY 
#        789         69        468        156
Analytical total and conditional permutation, and bootstrap conditional permutation standard deviates of local Moran’s I for first round turnout, row-standardised neighbours

Figure 15.3: Analytical total and conditional permutation, and bootstrap conditional permutation standard deviates of local Moran’s I for first round turnout, row-standardised neighbours

Figure 15.3 shows that conditional permutation scales back the proportion of standard deviate values taking extreme values, especially positive values. The analytical total standard deviates of local Moran’s \(I\) should probably not be used since alternatives are available, not least thanks to the clarification by Sauer et al. (2021).

In presenting local Moran’s \(I\), use is often made of “hotspot” maps. Because \(I_i\) takes high values both for strong positive autocorrelation of low and high values of the input variable, it is hard to show where “clusters” of similar neighbours with low or high values of the input variable occur. The quadrants of the Moran plot are used, by creating a categorical quadrant variable interacting the input variable and its spatial lag split at their means. The quadrant categories are then set to NA if, for the chosen standard deviate and adjustment, \(I_i\) would be considered insignificant. Here, for the Bonferroni adjusted conditional analytical standard deviates, 14 observations belong to “Low-Low clusters”, and 55 to “High-High clusters”:

#                  Moran plot quadrants Unadjusted analytical total
# Low X : Low WX                   1040                         370
# High X : Low WX                   264                           3
# Low X : High WX                   213                           0
# High X : High WX                  978                         342
#                  Bonferroni analytical cond. Bonferroni cond. perm.
# Low X : Low WX                            14                     18
# High X : Low WX                            0                      0
# Low X : High WX                            0                      0
# High X : High WX                          55                     58
Local Moran’s I hotspot maps \(\alpha = 0.05\): left panel upper: unadjusted analytical standard deviates; right upper panel: Bonferroni adjusted analytical conditional standard deviates; left lower panel: Bonferroni adjusted bootstrap conditional permutation standard deviates, first round turnout, row-standardised neighbours

Figure 15.4: Local Moran’s I hotspot maps \(\alpha = 0.05\): left panel upper: unadjusted analytical standard deviates; right upper panel: Bonferroni adjusted analytical conditional standard deviates; left lower panel: Bonferroni adjusted bootstrap conditional permutation standard deviates, first round turnout, row-standardised neighbours

Figure 15.4 shows the impact of using analytical or conditional permutation standard deviates, and no or Bonferroni adjustment, reducing the counts of observations in “Low-Low clusters” from 370 to 14, and “High-High clusters” from 342 to 54; the “High-High clusters” are metropolitan areas.

Tiefelsdorf (2002) argues that standard approaches to the calculation of the standard deviates of local Moran’s \(I_i\) should be supplemented by numerical estimates, and shows that Saddlepoint approximations are a computationally efficient way of achieving this goal. The localmoran.sad() function takes a fitted linear model as its first argument, so we first fit a null (intercept only) model, but use case weights because the numbers entitled to vote vary greatly between observations:

Saddlepoint approximation is much more computationally intensive than conditional permutation:

#    user  system elapsed 
#  10.550   0.084  10.643

However, standard approaches do not permit richer mean models with covariates or case weights. Next we add the categorical variable distinguishing between rural, urban and other types of observational unit:

#    user  system elapsed 
#    11.4     0.0    11.4

To conclude, we add the spatially lagged categories, although the spatial lag of a categorical variable, represented by dummies in the model matrix, is not well defined:

#    user  system elapsed 
#  11.904   0.004  11.908
Saddlepoint weighted null model, weighted types model, weighted Durbin types model, and analytical conditional standard deviates of local Moran’s I for first round turnout, row-standardised neighbours

Figure 15.5: Saddlepoint weighted null model, weighted types model, weighted Durbin types model, and analytical conditional standard deviates of local Moran’s I for first round turnout, row-standardised neighbours

Figure 15.5 includes the analytical conditional standard deviates for comparison (lower right panel), but in general it can be seen that the Saddlepoint approximation standard deviates lie closer to zero, possibly because some of the mis-specification in the mean model has been removed by using richer versions, and possibly because the approximation approach is inherently local, relating regression residual values at \(i\) to those of its neighbours. It is also possible to use Saddlepoint approximation where the global spatial process has been incorporated, removing the conflation of global and local spatial autocorrelation in standard approaches. The same can also be accomplished using exact methods (Bivand, Müller, and Reder 2009).

15.3.2 Local Getis-Ord \(G_i\)

The local Getis-Ord \(G\) measure is reported as a standard deviate, and may also take the \(G^*_i\) form where self-neighbours are inserted into the neighbour object using include.self(). The observed and expected values of local \(G\) with their analytical variances may also be returned if return_internals=TRUE.

#    user  system elapsed 
#   0.027   0.000   0.036
#    user  system elapsed 
#   0.519   1.317   0.627

Once again we face the problem of multiple comparisons, with the count of areal unit p-values < 0.05 being reduced by an order of magnitude when employing Bonferroni correction:

#       none bonferroni        fdr         BY 
#        789         69        468        156
Plots of analytical conditional against bootstrap standard deviates; left: local Moran’s I; right: local G; first round turnout, row-standardised neighbours

Figure 15.6: Plots of analytical conditional against bootstrap standard deviates; left: local Moran’s I; right: local G; first round turnout, row-standardised neighbours

Figure 15.6 shows that, when using analytical conditional standard deviates, the values from analytical and bootstrap estimates coincide for both \(I_i\) and \(G_i\). In both cases, one may argue that the bootstrap approach is superfluous in exploratory spatial data analysis.

#             Low Not significant            High 
#              14            2426              55
Left: analytical standard deviates of local G; right: Bonferroni hotspots; first round turnout, row-standardised neighbours

Figure 15.7: Left: analytical standard deviates of local G; right: Bonferroni hotspots; first round turnout, row-standardised neighbours

As can be seen from Figure 15.6, we do not need to contrast the two estimation methods, and showing the mapped standard deviate is as informative as the “hotspot” status for the chosen adjustment (Figure 15.7). In the case of \(G_i\), the values taken by the measure reflect the values of the input variable, so a “High cluster” is found for observations with high values of the input variable, here high turnout in metropolitan areas.

Very recently, Geoda has been wrapped for R as rgeoda (Li and Anselin 2021), and will provide very similar functionalities for the exploration of spatial autocorrelation in areal data as spdep. The active objects are kept as pointers to a compiled code workspace; using compiled code for all operations (as in Geoda itself) makes rgeoda perform fast, but leaves less flexible when modifications or enhancements are desired.

The contiguity neighbours it constructs are the same as those found by poly2nb(), as almost are the \(I_i\) measures. The difference is as established by Bivand and Wong (2018), that localmoran() calculates the input variable variance divinding by \(n\), but Geoda uses \((n-1)\), as can be reproduced by setting mlvar=FALSE:

# here
#    user  system elapsed 
#   0.105   0.000   0.114
#                      name               value
# 1 number of observations:                2495
# 2          is symmetric:                 TRUE
# 3               sparsity: 0.00228786229774178
# 4        # min neighbors:                   1
# 5        # max neighbors:                  13
# 6       # mean neighbors:    5.70821643286573
# 7     # median neighbors:                   6
# 8           has isolates:               FALSE
#    user  system elapsed 
#   0.273   0.008   0.071
# [1] TRUE
# [1] TRUE

15.4 Exercises

  1. Why are join-count measures on a chessboard so different between rook and queen neighbours?
  2. Please repeat the simulation shown in section 15.1 using the chessboard polygons and the row-standardized queen contiguity neighbours. Why is it important to understand that spatial autocorrelation usually signals (unavoidable) mis-specification in our data?
  3. Do we need to use conditional permutation in inference for local measures of spatial autocorrelation?
  4. Why is false discovery rate adjustment recommended for local measures of spatial autocorrelation?
  5. Compare the local Moran’s \(I_i\) standard deviate values for the simulated data from exercise 15.2 for the analytical conditional approach, and Saddlepoint approximation. Consider the advantages and disadvantages of the Saddlepoint approximation approach.


Anselin, L. 1995. “Local indicators of spatial association - LISA.” Geographical Analysis 27 (2): 93–115.

Anselin, L. 1996. “The Moran Scatterplot as an ESDA Tool to Assess Local Instability in Spatial Association.” In Spatial Analytical Perspectives on GIS, edited by M. M. Fischer, H. J. Scholten, and D. Unwin, 111–25. London: Taylor & Francis.

Assunção, R.M., and E. A. Reis. 1999. “A New Proposal to Adjust Moran’s \(I\) for Population Density.” Statistics in Medicine 18: 2147–62.

Benjamini, Yoav, and Yosef Hochberg. 1995. “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society. Series B (Methodological) 57 (1): 289–300.

Benjamini, Yoav, and Daniel Yekutieli. 2001. “The control of the false discovery rate in multiple testing under dependency.” The Annals of Statistics 29 (4): 1165–88.

Bivand, Roger. 2021b. Spdep: Spatial Dependence: Weighting Schemes, Statistics.

Bivand, Roger S., and David W. S. Wong. 2018. “Comparing Implementations of Global and Local Indicators of Spatial Association.” TEST 27 (3): 716–48.

Bivand, R. S., W. Müller, and M. Reder. 2009. “Power Calculations for Global and Local Moran’s \(I\).” Computational Statistics and Data Analysis 53: 2859–72.

Brody, Howard, Michael Russell Rip, Peter Vinten-Johansen, Nigel Paneth, and Stephen Rachman. 2000. “Map-Making and Myth-Making in Broad Street: The London Cholera Epidemic, 1854.” The Lancet 356 (9223): 64–68.

Cliff, A. D., and J. K. Ord. 1973. Spatial Autocorrelation. London: Pion.

Cliff, A. 1981. Spatial Processes. London: Pion.

Duncan, O. D., R. P. Cuzzort, and B. Duncan. 1961. Statistical Geography: Problems in Analyzing Areal Data. Glencoe, IL: Free Press.

Geary, R. C. 1954. “The Contiguity Ratio and Statistical Mapping.” The Incorporated Statistician 5: 115–45.

Getis, A., and J. K. Ord. 1992. “The Analysis of Spatial Association by the Use of Distance Statistics.” Geographical Analysis 24 (2): 189–206.

Getis, A., and J. K. Ord. 1992. “The Analysis of Spatial Association by the Use of Distance Statistics.” Geographical Analysis 24 (2): 189–206.

1996. “Local Spatial Statistics: An Overview.” In Spatial Analysis: Modelling in a Gis Environment, edited by P. Longley and M Batty, 261–77. Cambridge: GeoInformation International.

Li, Xun, and Luc Anselin. 2021. Rgeoda: R Library for Spatial Data Analysis.

McMillen, Daniel P. 2003. “Spatial Autocorrelation or Model Misspecification?” International Regional Science Review 26: 208–17.

Moran, P. A. P. 1948. “The Interpretation of Statistical Maps.” Journal of the Royal Statistical Society, Series B (Methodological) 10 (2): 243–51.

Olsson, Gunnar. 1970. “Explanation, Prediction, and Meaning Variance: An Assessment of Distance Interaction Models.” Economic Geography 46: 223–33.

Ord, J. K., and A. Getis. 2001. “Testing for Local Spatial Autocorrelation in the Presence of Global Autocorrelation.” Journal of Regional Science 41 (3): 411–32.

Sauer, Jeffery, Taylor M Oshan, Sergio Rey, and Levi J Wolf. 2021. “On Null Hypotheses and Heteroskedasticity.” OSF Preprints.

Schabenberger, O., and C. A. Gotway. 2005. Statistical Methods for Spatial Data Analysis. Boca Raton/London: Chapman & Hall/CRC.

Sokal, R. R, N. L. Oden, and B. A. Thomson. 1998. “Local Spatial Autocorrelation in a Biological Model.” Geographical Analysis 30: 331–54.

Tiefelsdorf, M. 2002. “The Saddlepoint Approximation of Moran’s I and Local Moran’s \({I}_i\) Reference Distributions and Their Numerical Evaluation.” Geographical Analysis 34: 187–206.

Tobler, W. R. 1970. “A Computer Movie Simulating Urban Growth in the Detroit Region.” Economic Geography 46: 234–40.

Upton, G., and B. Fingleton. 1985. Spatial Data Analysis by Example: Point Pattern and Qualitative Data. New York: Wiley.