# Chapter 13 Area Data and Spatial Autocorrelation

Areal units of observation are very often used when simultaneous observations are aggregated within non-overlapping boundaries. The boundaries may be those of administrative entities, and may be related to underlying spatial processes, such as commuting flows, but are usually arbitrary. If they do not match the underlying and unobserved spatial processes in one or more variables of interest, proximate areal units will contain parts of the underlying processes, engendering spatial autocorrelation. This is at least in part because the aggregated observations are driven by autocorrelated factors which may or may not themselves have been observed.

With support of data we mean the physical size (lenth, area, volume) associated with an individual observational unit (measurement). It is possible to represent the support of areal data by a point, despite the fact that the data have polygonal support. The centroid of the polygon may be taken as a representative point, or the centroid of the largest polygon in a multi-polygon object. When data with intrinsic point support are treated as areal data, the change of support goes the other way, from the known point to a non-overlapping tesselation such as a Voronoi diagram or Dirichlet tessellation or Thiessen polygons often through a Delaunay triangulation and using a Euclidean plane (projected coordinates). Here, different metrics may also be chosen, or distances measured on a network rather than on the plane. There is also a literature using weighted Voronoi diagrams in local spatial analysis (see for example Boots and Okabe 2007; Okabe et al. 2008; She et al. 2015).

When the intrinsic support of the data is as points, but the underlying process is between proximate observations rather than driven chiefly by distance however measured between observations, the data may be aggregate counts or totals (polling stations, retail turnover) or represent a directly observed characteristic of the observation (opening hours of the polling station). Obviously, the risk of mis-representing the footprint of the underlying spatial processes remains in all of these cases, not least because the observations are taken as encompassing the entirety of the underlying process in the case of tesselation of the whole area of interest. This is distinct from the geostatistical setting in which observations are rather samples taken using some scheme within the area of interest. It is also partly distinct from the practice of taking areal sample plots within the area of interest but covering only a small proportion of the area, typically used in ecological and environmental research.

This chapter then considers a subset of the methods potentially available for exploring spatial autocorrelation in areal data, or data being handled as areal, where the spatial processes are considered as working through proximity understood in the first instance as contiguity, as a graph linking observations taken as neighbours. This graph is typically undirected and unweighted, but may be directed and/or weighted in certain settings, which then leads to further issues with regard to symmetry. In principle, proximity would be expected to operate symmetrically in space, that is that the influence of \(i\) on \(j\) and of \(j\) on \(i\) based on their relative positions should be equivalent. Edge effects are not considered in standard treatments.

## 13.1 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).

## 13.2 Spatial weight matrices

Handling spatial autocorrelation using relationships to neighbours on a graph takes the graph as given, chosen by the analyst. This differs from the geostatistical approach in which the analyst chooses the binning of the empirical variogram and function used, and then the way the fitted variogram is fitted. Both involve a priori choices, but represent the underlying correlation in different ways (Wall 2004). In Bavaud (1998) and work citing his contribution, attempts have been made to place graph-based neighbours in a broader context.

One issue arising in the creation of objects representing neighbourhood relationships is that of no-neighbour areal units (Bivand and Portnov 2004). Islands or units separated by rivers may not be recognised as neighbours when the units have areal support and when using topological relationships such as shared boundaries. In some settings, for example `mrf`

(Markov Random Field) terms in `mgcv::gam()`

and similar model fitting functions that require undirected connected graphs, a requirement is violated when there are disconnected subgraphs.

No-neighbour observations can also occur when a distance threshold is used between points, where the threshold is smaller than the maximum nearest neighbour distance. Shared boundary contiguities are not affected by using geographical, unprojected coordinates, but all point-based approaches use distance in one way or another, and need to calculate distances in an appropriate way.

The **spdep** package provides an `nb`

class for neighbours, a list of length equal to the number of observations, with integer vector components. No-neighbours are encoded as an integer vector with a single element `0L`

, and observations with neighbours as sorted integer vectors containing values in `1L:n`

pointing to the neighbouring observations. This is a typical row-oriented sparse representation of neighbours. **spdep** provides many ways of constructing `nb`

objects, and the representation and construction functions are widely used in other packages.

**spdep** builds on the `nb`

representation (undirected or directed graphs) with the `listw`

object, a list with three components, an `nb`

object, a matching list of numerical weights, and a single element character vector containing the single letter name of the way in which the weights were calculated. The most frequently used approach in the social sciences is calculating weights by row standardization, so that all the non-zero weights for one observation will be the inverse of the cardinality of its set of neighbours (`1/card(nb[[i]])`

).

We will be using election data from the 2015 Polish Presidential election in this chapter, with 2495 municipalities and Warsaw boroughs (see Figure 13.1 for a **tmap** map (section 9.5) of the municipality types) , and complete count data from polling stations aggregated to these areal units. The data are an **sf** `sf`

object:

```
# Simple feature collection with 6 features and 3 fields
# Geometry type: MULTIPOLYGON
# Dimension: XY
# Bounding box: xmin: 235000 ymin: 367000 xmax: 281000 ymax: 413000
# Projected CRS: ETRS89 / Poland CS92
# TERYT name types geometry
# 1 020101 BOLESŁAWIEC Urban MULTIPOLYGON (((261089 3855...
# 2 020102 BOLESŁAWIEC Rural MULTIPOLYGON (((254150 3837...
# 3 020103 GROMADKA Rural MULTIPOLYGON (((275346 3846...
# 4 020104 NOWOGRODZIEC Urban/rural MULTIPOLYGON (((251770 3770...
# 5 020105 OSIECZNICA Rural MULTIPOLYGON (((263424 4060...
# 6 020106 WARTA BOLESŁAWIECKA Rural MULTIPOLYGON (((267031 3870...
```

Between early 2002 and April 2019, **spdep** contained functions for constructing and handling neighbour and spatial weights objects, tests for spatial autocorrelation, and model fitting functions. The latter have been split out into **spatialreg**, and will be discussed in the next chapter. **spdep** now accommodates objects represented using **sf** classes and **sp** classes directly.

`# Loading required package: sp`

`# Loading required package: spData`

### 13.2.1 Contiguous neighbours

The `poly2nb()`

function in **spdep** takes the boundary points making up the polygon boundaries in the object passed as the `pl=`

argument, and for each observation checks whether at least one (`queen=TRUE`

, default), or at least two (rook, `queen=FALSE`

) points are within `snap=`

distance units of each other. The distances are planar in the raw coordinate units, ignoring geographical projections. Once the required number of sufficiently close points is found, the search is stopped.

```
# function (pl, row.names = NULL, snap = sqrt(.Machine$double.eps),
# queen = TRUE, useC = TRUE, foundInBox = NULL, small_n = 500)
# NULL
```

Two other arguments are also worth discussing, `foundInBox=`

and `small_n=`

. The first, `foundInBox=`

, accepted the output of the **rgeos** `gUnarySTRtreeQuery()`

function to list candidate neighbours, that is polygons whose bounding boxes intersect the bounding boxes of other polygons. From **spdep** 1.1-7, the GEOS interface of the **sf** package is used within `poly2nb()`

if `foundInBox=NULL`

and the number of observations is greater than `small_n=`

, to find the candidate neighbours and populate `foundInBox`

internally. In this case, this use of spatial indexing (STRtree queries) in GEOS through **sf** is the default, as the number of observations is greater than `small_n`

:

```
# user system elapsed
# 0.859 0.000 0.860
```

The print method shows the summary structure of the neighbour object:

```
# Neighbour list object:
# Number of regions: 2495
# Number of nonzero links: 14242
# Percentage nonzero weights: 0.229
# Average number of links: 5.71
```

Raising `small_n=`

above the observation count, we see that the processing time is increased a little; the benefits of indexing are more apparent with larger data sets.

```
# user system elapsed
# 0.998 0.000 0.998
```

The output objects are identical:

`# [1] TRUE`

Much of the work involved in finding contiguous neighbours is spent on finding candidate neighbours with intersecting bounding boxes. Note that `nb`

objects record both symmetric neighbour relationships, because these objects admit asymmetric relationships as well, but these duplications are not needed for object construction.

Most of the **spdep** functions for constructing neighbour objects take a `row.names=`

argument, the value of which is stored as a `region.id`

attribute. If not given, the values are taken from `row.names()`

of the first argument. These can be used to check that the neighbours object is in the same order as data. If `nb`

objects are subsetted, the indices change to continue to be within `1:length(subsetted_nb)`

, but the `region.id`

attribute values point back to the object from which it was constructed.

We can also check that this undirected graph is connected using the `n.comp.nb()`

function; while some model estimation techniques do not support graphs that are not connected, it is helpful to be aware of possible problems (Freni-Sterrantino, Ventrucci, and Rue 2018):

`# [1] 1`

Neighbour objects may be exported and imported in GAL format for exchange with other software, using `write.nb.gal()`

and `read.gal()`

:

### 13.2.2 Graph-based neighbours

If areal units are an appropriate representation, but only points have been observed, contiguity relationships may be approximated using graph-based neighbours. In this case, the imputed boundaries tesselate the plane such that points closer to one observation than any other fall within its polygon. The simplest form is by using triangulation, here using the `deldir()`

function in the **deldir** package. Because the function returns from and to identifiers, it is easy to construct a long representation of a `listw`

object, as used in the S-Plus SpatialStats module and the `sn2listw()`

function internally to construct an `nb`

object (ragged wide representation). Alternatives often fail to return sufficient information to permit the neighbours to be identified.

The `soi.graph()`

function takes triangulated neighbours and prunes off neighbour relationships represented by unusually long edges, especially around the convex hull, but may render the output object asymmetric. Other graph-based approaches include `relativeneigh()`

and `gabrielneigh()`

.

The output of these functions is then converted to the `nb`

representation using `graph2nb()`

, with the possible use of the `sym=`

argument to coerce to symmetry. We take the centroids of the largest component polygon for each observation as the point representation; population-weighted centroids might have been a better choice if they were available:

```
coords <- st_centroid(st_geometry(pol_pres15), of_largest_polygon=TRUE)
suppressMessages(nb_tri <- tri2nb(coords))
nb_tri
```

```
# Neighbour list object:
# Number of regions: 2495
# Number of nonzero links: 14930
# Percentage nonzero weights: 0.24
# Average number of links: 5.98
```

The average number of neighbours is similar to the Queen boundary contiguity case, but if we look at the distribution of edge lengths using `nbdists()`

, we can see that although the upper quartile is about 15 km, the maximum is almost 300 km, an edge along much of one side of the convex hull. The short minimum distance is also of interest, as many centroids of urban municipalities are very close to the centroids of their surrounding rural counterparts.

```
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 247 9847 12151 13485 14994 296974
```

Triangulated neighbours also yield a connected graph:

`# [1] 1`

The sphere of influence graph trims a neighbour object such as `nb_tri`

to remove edges that seem long in relation to typical neighbours (Avis and Horton 1985).

```
# Neighbour list object:
# Number of regions: 2495
# Number of nonzero links: 12792
# Percentage nonzero weights: 0.205
# Average number of links: 5.13
```

Unpicking the triangulated neighbours does however remove the connected character of the underlying graph:

`# [1] 16`

The SoI algorithm has stripped out longer edges leading to urban and rural municipality pairs where their centroids are very close to each other because the rural ones completely surround the urban, giving 15 pairs of neighbours unconnected to the main graph:

```
#
# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
# 2465 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
```

The largest length edges along the convex hull have been removed, but “holes” have appeared where the unconnected pairs of neighbours have appeared. The differences between `nb_tri`

and `nb_soi`

are shown in orange in Figure 13.2.

```
opar <- par(mar=c(0,0,0,0)+0.5)
plot(st_geometry(pol_pres15), border="grey", lwd=0.5)
plot(nb_soi, coords=st_coordinates(coords), add=TRUE, points=FALSE, lwd=0.5)
plot(diffnb(nb_tri, nb_soi), coords=st_coordinates(coords), col="orange",
add=TRUE, points=FALSE, lwd=0.5)
```

### 13.2.3 Distance-based neighbours

Distance-based neighbours can be constructed using `dnearneigh()`

, with a distance band with lower `d1=`

and upper `d2=`

bounds controlled by the `bounds=`

argument. If spherical coordinates are used and either specified in the coordinates object `x`

or with `x`

as a two column matrix and `longlat=TRUE`

, great circle distances in km will be calculated assuming the WGS84 reference ellipsoid. From **spdep** 1.1-7, two arguments have been added, to use functionality in the **dbscan** package for finding neighbours using a kd-tree in two or three dimensions by default, and not to test the symmetry of the output neighbour object.

The `knearneigh()`

function for \(k\)-nearest neighbours returns a `knn`

object, converted to an `nb`

object using `knn2nb()`

. It can also use great circle distances, not least because nearest neighbours may differ when uprojected coordinates are treated as planar. `k=`

should be a small number. For projected coordinates, the **dbscan** package is used to compute nearest neighbours more efficiently. Note that `nb`

objects constructed in this way are most unlikely to be symmetric, hence `knn2nb()`

has a `sym=`

argument to permit the imposition of symmetry, which will mean that all units have at least `k=`

neighbours, not that all units will have exactly `k=`

neighbours.

The `nbdists()`

function returns the length of neighbour relationship edges in the units of the coordinates if the coordinates are projected, in km otherwise. In order to set the upper limit for distance bands, one may first find the maximum first nearest neighbour distance, using `unlist()`

to remove the list structure of the returned object.

```
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 247 6663 8538 8275 10124 17979
```

Here the largest first nearest neighbour distance is just under 18 km, so using this as the upper threshold gives certainty that all units will have at least one neighbour:

```
# user system elapsed
# 0.168 0.000 0.168
```

```
# user system elapsed
# 0.152 0.000 0.152
```

`# [1] TRUE`

```
# Neighbour list object:
# Number of regions: 2495
# Number of nonzero links: 20358
# Percentage nonzero weights: 0.327
# Average number of links: 8.16
```

However, even though there are no no-neighbour observations (their presence is reported by the print method for `nb`

objects), the graph is not connected, as a pair of observations are each others’ only neighbours.

`# [1] 2`

```
#
# 1 2
# 2493 2
```

Adding 300 m to the threshold gives us a neighbour object with no no-neighbour units, and all units can be reached from all others across the graph.

```
# Neighbour list object:
# Number of regions: 2495
# Number of nonzero links: 21086
# Percentage nonzero weights: 0.339
# Average number of links: 8.45
```

`# [1] 1`

One characteristic of distance-based neighbours is that more densely settled areas, with units which are smaller in terms of area (Warsaw boroughs are much smaller on average, but have almost 30 neighbours). Having many neighbours smooths the neighbour relationship across more neighbours. For use later, we also construct a neighbour object with no-neighbour units, using a threshold of 16 km:

```
# Neighbour list object:
# Number of regions: 2495
# Number of nonzero links: 15850
# Percentage nonzero weights: 0.255
# Average number of links: 6.35
# 7 regions with no links:
# 569 1371 1522 2374 2385 2473 2474
```

It is possible to control the numbers of neighbours directly using \(k\)-nearest neighbours, either accepting asymmetric neighbours or imposing symmetry:

```
# Neighbour list object:
# Number of regions: 2495
# Number of nonzero links: 14970
# Percentage nonzero weights: 0.24
# Average number of links: 6
# Non-symmetric neighbours list
```

```
# Neighbour list object:
# Number of regions: 2495
# Number of nonzero links: 16810
# Percentage nonzero weights: 0.27
# Average number of links: 6.74
```

Here the size of `k=`

is sufficient to ensure connectedness, although the graph is not planar as edges cross at locations other than nodes, which is not the case for contiguous or graph-based neighbours.

`# [1] 1`

### 13.2.4 Weights specification

Once neighbour objects are available, further choices need to made in specifying the weights objects. The `nb2listw()`

function is used to create a `listw`

weights object with an `nb`

object, a matching list of weights vectors, and a style specification. Because handling no-neighbour observations now begins to matter, the `zero.policy=`

argument is introduced. By default, this is `FALSE`

, indicating that no-neighbour observations will cause an error, as the spatially lagged value for an observation with no neighbours is not available. By convention, zero is substituted for the lagged value, as the cross product of a vector of zero-valued weights and a data vector, hence the name of `zero.policy`

.

```
# function (neighbours, glist = NULL, style = "W", zero.policy = NULL)
# NULL
```

We will be using the helper function `spweights.constants()`

below to show some consequences of varing style choices. It returns constants for a `listw`

object, \(n\) is the number of observations, `n1`

to `n3`

are \(n-1, \ldots\), `nn`

is \(n^2\) and \(S_0\), \(S_1\) and \(S_2\) are constants, \(S_0\) being the sum of the weights. There is a full discussion of the constants in Bivand and Wong (2018).

```
# function (listw, zero.policy = NULL, adjust.n = TRUE)
# NULL
```

The `"B"`

binary style gives a weight of unity to each neighbour relationship, and typically upweights units with no boundaries on the edge of the study area.

```
# n n1 n2 n3 nn S0 S1 S2
# 2495 2494 2493 2492 6225025 14242 28484 357280
```

The `"W"`

row-standardized style upweights units around the edge of the study area that necessarily have fewer neighbours. This style first gives a weight of unity to each neighbour relationship, then divides these weights by the per unit sums of weights. Naturally this leads to division by zero where there are no neighbours, a not-a-number result, unless the chosen policy is to permit no-neighbour observations. We can see that \(S_0\) is now equal to \(n\).

```
# n S0 S1 S2
# 2495 2495 958 10406
```

Inverse distance weights are used in a number of scientific fields. Some use dense inverse distance matrices, but many of the inverse distances are close to zero, so have little practical contribution, especially as the spatial process matrix is itself dense. Inverse distance weights may be constructed by taking the lengths of edges, changing units to avoid most weights being too large or small (here from m to km), taking the inverse, and passing through the `glist=`

argument to `nb2listw()`

:

```
gwts <- lapply(nbdists(nb_d183, coords), function(x) 1/(x/1000))
lw_d183_idw_B <- nb2listw(nb_d183, glist=gwts, style="B")
unlist(spweights.constants(lw_d183_idw_B))[c(1,6:8)]
```

```
# n S0 S1 S2
# 2495 1841 534 7265
```

No-neighbour handling is by default to prevent the construction of a weights object, making the analyst take a position on how to proceed.

`# Error in nb2listw(nb_d16, style = "B") : Empty neighbour sets found`

Use can be made of the `zero.policy=`

argument to many functions used with `nb`

and `listw`

objects.

```
lw_d16_B <- nb2listw(nb_d16, style="B", zero.policy=TRUE)
unlist(spweights.constants(lw_d16_B, zero.policy=TRUE))[c(1,6:8)]
```

```
# n S0 S1 S2
# 2488 15850 31700 506480
```

Note that by default the `adjust.n=`

argument to `spweights.constants()`

is set by default to `TRUE`

, subtracting the count of no-neighbour observations from the observation count, so \(n\) is smaller with possible consequences for inference. The complete count can be retrieved by changing the argument.

## 13.3 Measures of spatial autocorrelation

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`

):

```
set.seed(1)
x <- rnorm(nrow(pol_pres15))
mt <- moran.test(x, lw_q_B, randomisation=FALSE, alternative="two.sided")
```

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:

```
glance_htest <- function(ht) c(ht$estimate, "Std deviate"=unname(ht$statistic), "p.value"=unname(ht$p.value))
glance_htest(mt)
```

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

The test however fails to detect a missing trend in the data as a missing variable problem, finding spatial autocorrelation instead:

```
beta <- 0.15e-02
t <- st_coordinates(coords)[,1]/1000
x_t <- x + beta*t
mt <- moran.test(x_t, lw_q_B, randomisation=FALSE, alternative="two.sided")
glance_htest(mt)
```

```
# 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, 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).

### 13.3.1 Global measures

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)
# NULL
```

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 are not unity in all cases, 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
```

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=`

agrument 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 = na.fail,
# spChk = NULL, adjust.n = TRUE, drop.EI2 = FALSE)
# NULL
```

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 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
```

Of course, from the very beginning, interest was shown in Monte Carlo testing, also known as a Hope-type test and as a permutation bootstrap. By default, `moran.mc()`

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 of the variable of interest by permutation, that is shuffling the values across the observations at random, needs to be given as `nsim=`

.

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 `EBImoran.mc()`

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. A key source of mis-specification will typically also include the choice of entities for aggregation of data.

### 13.3.2 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 13.3, 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.

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

Bivand and Wong (2018) discuss issues impacting the use of local indicators, such as local Moran’s \(I\) and local Getis-Ord \(G\). 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). 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 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 because these local indicators when summed and divided by the sum of the spatial weights, and thus positive and negative local spatial autocorrelation may be present, we obtain:

```
all.equal(sum(locm[,1])/Szero(lw_q_W),
unname(moran.test(pol_pres15$I_turnout, lw_q_W)$estimate[1]))
```

`# [1] TRUE`

Using `stats::p.adjust()`

to adjust for multiple comparisons, we see that
almost 29% of the local measures have p-values < 0.05 if no adjustment is applied, but only 12% using Bonferroni adjustment, with two other choices also shown:

```
pvs <- cbind("none"=locm[,5], "bonferroni"=p.adjust(locm[,5], "bonferroni"),
"fdr"=p.adjust(locm[,5],"fdr"), "BY"=p.adjust(locm[,5], "BY"))
apply(pvs, 2, function(x) sum(x < 0.05))
```

```
# none bonferroni fdr BY
# 715 297 576 424
```

In the global measure case, bootstrap permutations could be used as an alternative to analytical methods for possible inference. In the local case, conditional permutation may be used, retaining 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 will 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:

`# NULL`

```
# user system elapsed
# 0.599 1.506 0.767
```

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.

```
pv <- locm_p[,5]
pvsp <- cbind("none"=pv, "bonferroni"=p.adjust(pv, "bonferroni"), "fdr"=p.adjust(pv, "fdr"), "BY"=p.adjust(pv, "BY"))
apply(pvsp, 2, function(x) sum(x <= 0.05))
```

```
# 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 68 in the conditional permutation case:

```
brks <- qnorm(c(0, 0.00001, 0.0001, 0.001, 0.01, 0.025, 0.5, 0.975, 0.99, 0.999, 0.9999, 0.99999, 1))
(tab <- table(cut(locm_p[,4], brks)))
```

```
#
# (-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`

```
pol_pres15$locm_Z <- locm[,4]
pol_pres15$locm_p_Z <- locm_p[,4]
tm_shape(pol_pres15) + tm_fill(c("locm_Z", "locm_p_Z"), breaks=brks, midpoint=0, title="Standard deviates of\nLocal Moran's I") + tm_facets(free.scales=FALSE) + tm_layout(panel.labels=c("Analytical", "Conditional permutation"))
```

Figure 13.5 shows that conditional permutation scales back the proportion of standard deviate values taking extreme values, especially positive values. As we will see below, the analytical standard deviates of local Moran’s \(I\) should probably not be used if alternatives are available.

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 conditional permutation standard deviates, Bonferroni adjusted, 14 observations belong to “Low-Low clusters”, and 54 to “High-High clusters”:

```
quadr <- interaction(cut(infl_W$x, c(-Inf, mean(infl_W$x), Inf), labels=c("Low X", "High X")), cut(infl_W$wx, c(-Inf, mean(infl_W$wx), Inf), labels=c("Low WX", "High WX")), sep=" : ")
a <- table(quadr)
pol_pres15$hs_an_q <- quadr
is.na(pol_pres15$hs_an_q) <- !(pol_pres15$locm_Z < brks[6] | pol_pres15$locm_Z > brks[8])
b <- table(pol_pres15$hs_an_q)
pol_pres15$hs_cp_q <- quadr
is.na(pol_pres15$hs_cp_q) <- !(pol_pres15$locm_p_Z < brks[2] | pol_pres15$locm_p_Z > brks[12])
c <- table(pol_pres15$hs_cp_q)
t(rbind("Moran plot quadrants"=a, "Unadjusted analytical"=b, "Bonferroni cond. perm."=c))
```

```
# Moran plot quadrants Unadjusted analytical
# 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 cond. perm.
# Low X : Low WX 18
# High X : Low WX 0
# Low X : High WX 0
# High X : High WX 58
```

Figure 13.6 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.

The local Getis-Ord \(G\) measure is reported as a standard deviate, and may also take the \(G^*\) 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.014 0.000 0.013
```

```
# user system elapsed
# 0.123 0.762 0.788
```

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:

```
pv <- 2 * pnorm(abs(c(locG)), lower.tail = FALSE)
pvsp <- cbind(pv, p.adjust(pv, "bonferroni"), p.adjust(pv, "fdr"),
p.adjust(pv, "BY"))
colnames(pvsp) <- c("none", "bonferroni", "fdr", "BY")
apply(pvsp, 2, function(x) sum(x < 0.05))
```

```
# none bonferroni fdr BY
# 789 69 468 156
```

In the \(G_i\) case, however, there is no systematic differencce between the analytical and conditional permutation standard deviates.

```
library(ggplot2)
p1 <- ggplot(data.frame(Zi=locm[,4], Zi_perm=locm_p[,4])) + geom_point(aes(x=Zi, y=Zi_perm), alpha=0.2) + xlab("Analytical") + ylab("Conditional permutation") + coord_fixed() + ggtitle("Local Moran's I")
p2 <- ggplot(data.frame(Zi=c(locG), Zi_perm=c(locG_p))) + geom_point(aes(x=Zi, y=Zi_perm), alpha=0.2) + xlab("Analytical") + ylab("Conditional permutation") + coord_fixed() + ggtitle("Local G")
gridExtra::grid.arrange(p1, p2, nrow=1)
```

Figure 13.7 shows that, keeping fixed aspect in both panels, conditional permutation changes the range and distribution of the standard deviate values for \(I_i\), but that for \(G_i\), the two sets of standard deviates are equivalent.

```
pol_pres15$locG_Z <- c(locG)
pol_pres15$hs_G <- cut(c(locG), c(-Inf, brks[2], brks[12], Inf), labels=c("Low", "Not significant", "High"))
table(pol_pres15$hs_G)
```

```
#
# Low Not significant High
# 14 2426 55
```

```
m1 <- tm_shape(pol_pres15) + tm_fill(c("locG_Z"), midpoint=0, title="Standard\ndeviate")
m2 <- tm_shape(pol_pres15) + tm_fill(c("hs_G"), title="Bonferroni\nhotspot status")
tmap_arrange(m1, m2, nrow=1)
```

As can be seen from Figure 13.7, 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 13.8). 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`

:

`# Loading required package: digest`

```
#
# Attaching package: 'rgeoda'
```

```
# The following object is masked from 'package:spdep':
#
# skater
```

`# here`

```
# user system elapsed
# 0.111 0.004 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
```

```
system.time(lisa <- local_moran(Geoda_w, pol_pres15["I_turnout"],
cpu_threads=ifelse(parallel::detectCores() == 1, 1, parallel::detectCores()-1L), permutations=499, seed=1))
```

```
# user system elapsed
# 0.288 0.009 0.079
```

`# [1] TRUE`

`# [1] TRUE`

### References

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.

Avis, D., and J. Horton. 1985. “Remarks on the Sphere of Influence Graph.” In *Discrete Geometry and Convexity*, edited by J. E. Goodman, 323–27. New York: New York Academy of Sciences, New York.

Bavaud, F. 1998. “Models for Spatial Weights: A Systematic Look.” *Geographical Analysis* 30: 153–71. https://doi.org/10.1111/j.1538-4632.1998.tb00394.x.

Bivand, Roger S., and David W. S. Wong. 2018. “Comparing Implementations of Global and Local Indicators of Spatial Association.” *TEST* 27 (3): 716–48. https://doi.org/10.1007/s11749-018-0599-x.

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.

Bivand, R. S., and B. A. Portnov. 2004. “Exploring Spatial Data Analysis Techniques Using R: The Case of Observations with No Neighbours.” In *Advances in Spatial Econometrics: Methodology, Tools, Applications*, edited by L. Anselin, R. J. G. M. Florax, and S. J. Rey, 121–42. Berlin: Springer.

Boots, B., and A. Okabe. 2007. “Local Statistical Spatial Analysis: Inventory and Prospect.” *International Journal of Geographical Information Science* 21 (4): 355–75. https://doi.org/10.1080/13658810601034267.

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.

Freni-Sterrantino, Anna, Massimo Ventrucci, and Håvard Rue. 2018. “A Note on Intrinsic Conditional Autoregressive Models for Disconnected Graphs.” *Spatial and Spatio-Temporal Epidemiology* 26: 25–34. https://doi.org/https://doi.org/10.1016/j.sste.2018.04.002.

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.

*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*. https://CRAN.R-project.org/package=rgeoda.

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.

Okabe, A., T. Satoh, T. Furuta, A. Suzuki, and K. Okano. 2008. “Generalized Network Voronoi Diagrams: Concepts, Computational Methods, and Applications.” *International Journal of Geographical Information Science* 22 (9): 965–94. https://doi.org/10.1080/13658810701587891.

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.

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

She, Bing, Xinyan Zhu, Xinyue Ye, Wei Guo, Kehua Su, and Jay Lee. 2015. “Weighted Network Voronoi Diagrams for Local Spatial Analysis.” *Computers, Environment and Urban Systems* 52: 70–80. https://doi.org/https://doi.org/10.1016/j.compenvurbsys.2015.03.005.

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.

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

Wall, M. M. 2004. “A Close Look at the Spatial Structure Implied by the CAR and SAR Models.” *Journal of Statistical Planning and Inference* 121: 311–24.