Section 9.2 Exploring spatial autocorrelation in lattice data
Subsection 9.2.1 Burglary in Manchester
To illustrate spatial autocorrelation in lattice data, we are going to hop across the globe to Manchester, UK. To follow, read in the file
burglary_manchester_lsoa.geojson which contains burglaries per Lower Layer Super Output Area (LSOA), which is the data with which we will be working here.
burglary <- st_read("data/burglary_manchester_lsoa.geojson",
quiet=TRUE)
To have a look at the data, we can plot it quickly with
tmap:
tm_shape(burglary) +
tm_bubbles("burglary", border.lwd=NA, perceptual = TRUE) +
tm_borders(alpha=0.1) +
tm_layout(legend.position = c("right", "bottom"),
legend.title.size = 0.8,
legend.text.size = 0.5)
Do you see any patterns? Are burglaries randomly spread around the map? Or would you say that areas that are closer to each other tend to be more alike? Is there evidence of clustering? Do burglaries seem to appear in certain pockets of the map? To answer these questions, we will need to consider whether we see similar values in neighbouring areas.
Subsection 9.2.2 What is a neighbour?
Previously we asked whether areas are similar to their neighbours or to areas that are close. But what is a neighbour? Or what do we mean by close? How can one define a set of neighbours for each area? If we want to know if what we measure in a particular area is similar to what happens on its neighbouring areas, we first need to establish what we mean by a neighbour.
Subsubsection 9.2.2.1 Contiguity-based neighbours
There are various ways of defining a neighbour. Most rely on topological or geometrical relationships among the areas. First, we can say that two areas are neighbours if they share boundaries, if they are next to each other. In this case we talk of neighbours by contiguity. Contiguous can, at the same time, mean all areas that share common boundaries (what we call contiguity using the rook criteria, like in chess) or areas that share common boundaries and common corners, that is, that have any point in common (and we call this contiguity using the queen criteria). When we use this criteria we can refine our definition by defining the intensity of neighbourliness "as a function of the length of the common border as a proportion of the total border" [178] (p. 89). Given that in criminology we mostly work with irregular lattice data (such as police districts or census geographies), the queen criteria makes more sense. There is little theoretical justification for why, say, a police district would only exhibit dependence according to the rook criteria.
When defining neighbours by contiguity, we may also specify the order of contiguity. First order contiguity means that we are focusing on areas immediately contiguous (those in dark blue in the figure below). Second order means that we consider neighbours only those areas that are immediately contiguous to our first order neighbours (the light blue areas in the figure below) and you could go on and on. Look at the graphic below for clarification:
Subsubsection 9.2.2.2 Distance-based neighbours
We may also define neighbours by geographical distance. You can consider neighbours those areas that distance-wise are close to each other (regardless of whether boundaries are shared). The distance metric often will be Euclidean, but other commonly used distance metrics could also be employed (for example, Manhattan). Often one takes the centroid of the polygon as the location to take the distance from. More sophisticated approaches may involve taking a population weighted centroid, where the location of the centroid depends on the distribution of the population within an area.
Other approaches to define neighbours are graph-based, attribute-based, or interaction-based (e.g., flow of goods or people between areas). Contiguity and distance are the two most commonly used methods, though.
You will come across the term spatial weight matrix at some point or, using mathematical notation, \(W\text{.}\) Essentially the spatial weight matrix is a \(n\) by \(n\) matrix with, often, ones and zeroes (in the case of contiguity-based definitions) identifying if any two observations are neighbours. So you can think of the spatial weight matrix as the new data table that we are constructing with our definition of neighbours (whichever definition that is).
When working with contiguity measures, the problem of "islands" sometimes arises. The spatial weight matrix in this case would have zeroes for the row representing this "island". This is not permitted in the subsequent statistical use of this kind of matrix. A common solution is to "join" islands to the spatial units in the "mainland" that are closer.
So how do you build such a matrix with R? Well, let’s turn to that. But to make things a bit simpler, let’s focus not on the whole of Manchester, but just in the LSOAs within the city centre. Calculating a spatial weights matrix is a computationally intensive process, which means it takes a long time. The larger area you have (which will have more LSOAs), the longer this will take. We will use code we covered in Chapter 2 to clip the spatial object with the counts of burglaries to only those that intersect with the City Centre ward. If any of the following are unclear, visit Chapter 2: Basic Geospatial Operations in R in this book.
# Read a geojson file with Manchester wards and select only the city centre ward
city_centre <- st_read("data/manchester_wards.geojson", quiet = TRUE) %>%
filter(wd16nm == "City Centre")
# Reproject the burglary data
burglary <- st_transform(burglary, 27700)
# Intersect
cc_intersects <- st_intersects(city_centre, burglary)
cc_burglary <- burglary[unlist(cc_intersects),]
# Remove redundant objects
rm(city_centre)
rm(cc_intersects)
We can now have a look at our data. Let’s use the
view mode of tmap to enable us to interact with the map.
# Plot with tmap
tmap_mode("view")
tm_shape(cc_burglary) +
tm_fill("burglary", style = "quantile", palette = "Reds", id="code") +
tm_borders() +
tm_layout(main.title = "Burglary counts", main.title.size = 0.7 ,
legend.position = c("right", "top"), legend.title.size = 0.8)
So now we have a new spatial object "cc_burglary" with the 23 LSOA units that compose the City Centre of Manchester. By focusing in a smaller subset of areas, we can visualise perhaps a bit better what comes next.
Subsection 9.2.3 Creating a list of neighbours based on contiguity
Since it is an interactive map, if you are following along, you should be able to click on each LSOA, and see its unique ID. You could, manually go through and find the neighbours for each one this way. However, it would be very, very tedious having to identify the neighbours of all the areas in our study area by hand. That’s why we love computers. We can automate tedious work so that they do it, and we have more time to do fun stuff.
To calculate our neighbours, we are first going to turn our
sf object into spatial (sp) objects, so we can make use of the functions from the sp package that allow us to create a list of neighbours.
#We coerce the sf object into a new sp object
cc_burglary_sp <- as(cc_burglary, "Spatial")
In order to identify neighbours, we will use the
poly2nb() function from the spdep package that we loaded at the beginning of our session. The spdep package provides basic functions for building neighbour lists and spatial weights and tests for spatial autocorrelation for areal data like Moran’s I (more on this below). The poly2nb() function builds a neighbours list based on regions with contiguous boundaries. If you look at the documentation, you will see that you can pass a queen argument that takes TRUE or FALSE as options. If you do not specify this argument the default is set to true; that is, if you don’t specify queen = FALSE, this function will return a list of first order neighbours using the Queen criteria.
nb_queen <- poly2nb(cc_burglary_sp, row.names=cc_burglary_sp$code)
class(nb_queen)
## [1] "nb"
We see that we have created a
nb, neighbour list object called nb_queen. We can get some idea of what’s there if we ask for a summary.
summary(nb_queen)
## Neighbour list object: ## Number of regions: 23 ## Number of nonzero links: 100 ## Percentage nonzero weights: 18.9 ## Average number of links: 4.348 ## Link number distribution: ## ## 2 3 4 5 6 7 8 ## 3 4 6 5 3 1 1 ## 3 least connected regions: ## E01033673 E01033674 E01033684 with 2 links ## 1 most connected region: ## E01033658 with 8 links
This is basically telling us that using this criteria each LSOA polygon has an average of 4.3 neighbours (when we just focus on the city centre) and that all areas have some neighbours (there are no islands). The link number distribution gives you the number of links (neighbours) per area. So here we have 3 polygons with 2 neighbours, 3 with 3, 6 with 4, and so on. The summary function here also identifies the areas sitting at both extreme of the distribution. We can graphically represent the links using the following code:
# We first plot the boundaries
plot(cc_burglary_sp, col='white', border='blue', lwd=2)
# Then we use the coordinates function to obtain
# the coordinates of the polygon centroids
xy <- coordinates(cc_burglary_sp)
#Then we draw lines between the polygon centroids
# for neighbours that are listed as linked in nb_bur
plot(nb_queen, xy, col='red', lwd=2, add=TRUE)
Subsection 9.2.4 Generating the weight matrix
We can transform the object
nb_queen, the list of neighbours, into a spatial weights matrix. A spatial weights matrix reflects the intensity of the geographic relationship between observations. For this we use the spdep function nb2mat() (make a matrix out of a list of neighbours). The argument style sets what kind of matrix we want. B stands for the basic binary coding, you get a 1 for neighbours, and a 0 if the areas are not neighbours.
wm_queen <- nb2mat(nb_queen, style='B')
You can view this matrix by autoprinting the object
wm_queen which we just created:
wm_queen
This matrix has values of 0 or 1 indicating whether the elements listed in the rows are adjacent (using our definition, which in this case was the Queen criteria) with each other. The diagonal is full of zeroes. An area cannot be a neighbour of itself. So, if you look at the first two and the second column, you see a 1. That means that the LSOA with the code E01005065 is a neighbour of the second LSOA (as listed in the rows) which is E01005066. You will have zeroes for many of the other columns because this LSOA only has 4 neighbours.
In many computations we will see that the matrix is row standardised. We can obtain a row standardise matrix changing the code, specifically setting the
style to 'W' (the indicator for row standardised in this function):
wm_queen_rs <- nb2mat(nb_queen, style='W')
Row standardisation of a matrix ensures that the sum of the rows adds up to 1. So, for example, if you have 4 neighbours and that has to add up to 4, you need to divide 1 by 4, which gives you 0.25. So in the columns for a polygon with 4 neighbours you will see 0.25 in the column representing each of the neighbours.
Subsection 9.2.5 Spatial weight matrix based on distance
Subsubsection 9.2.5.1 Defining neighbours based on a critical threshold of distance
As noted earlier, we can also build the spatial weight matrix using some metric of distance. Remember that metrics such as Euclidean distance only makes sense if you have projected data. This is so because it does not take into account the curvature of the earth. If you have non-projected data, you can either transform the data into a projected system (what we have been doing in earlier examples) or you need to use distance metrics that work with the curvature of the Earth (e.g., arc-distance or great-circle distance). The simplest spatial distance-based weights matrix is obtained defining a given radius within each area and only considers as neighbours other areas that fall within such radius. So that we avoid "islands", the radius needs to be chosen in a way that guarantees that each location has at least one neighbour. So the first step would be to find an appropriate radius for establishing neighbourliness.
Previous to this, we need the centroids for the polygons. We already showed above how to obtain this using
coordinates() function in the sp package.
xy <- coordinates(cc_burglary_sp)
Then we use the
knearneigh() function from the spdep package to find the nearest neighbour to each centroid. This function will just do that, which is to look for the centroid closest to each of the other centroids. Then we turn the created object with this information into a nb object as defined earlier using knn2nb() also in spdep package:
# find nearest neighbours for each centroid coordinate pair
nearest_n <- knearneigh(xy)
# transform into nb class object
nb_nearest_n <- knn2nb(nearest_n)
Now we have a
nb list that for each area defines its nearest neighbour. What we need to do is to compute the distance from each centroid to its nearest neighbour (nbdists() will do this for us) and then find out what is the maximum distance (we can use max()), so that we use this as the critical threshold to define the radius we will use (to avoid the appearance of isolates in our matrix).
radius <- max(unlist(nbdists(nb_nearest_n, xy)))
radius
## [1] 690.9
Now that we know the maximum radius, we can use it to avoid isolates. The function
dnearneigh() will generate a list of neighbours within the radius we will specify. Notice the difference between knearneigh(), that we used earlier, and dnearneigh(). The first function generates a knn class object with the single nearest neighbour to each area \(i\text{,}\) whereas dnearneigh() generates a nb object with all the areas \(j\) within a given radius of area \(i\) being defined as neighbours. For dnearneigh() we need to define as parameters the object with the coordinates, the lower distance bound and the upper distance bound (the "radius" we just calculated). If you were to be working with non-projected data, you would need to specify longlat as TRUE. The function assumes we have projected data and, thus, we don’t need to specify this argument in our case, since our data are indeed using the British Grid projection. As we did with the nb objects generated through contiguity, we can use summary() to get a sense of the resulting list of neighbours.
nb_distance <- dnearneigh(xy, 0, radius)
summary(nb_distance)
## Neighbour list object: ## Number of regions: 23 ## Number of nonzero links: 88 ## Percentage nonzero weights: 16.64 ## Average number of links: 3.826 ## Link number distribution: ## ## 1 2 3 4 5 6 7 ## 3 4 4 4 2 3 3 ## 3 least connected regions: ## 1 17 23 with 1 link ## 3 most connected regions: ## 5 12 18 with 7 links
And as earlier, we can also plot the results:
# We first plot the boundaries
plot(cc_burglary_sp, col='white', border='blue', lwd=2)
# Then we plot the distance-based neighbours
plot(nb_distance, xy, lwd=2, col="red", add=TRUE)
We could now use this
nb object to generate a spatial weight matrix in the same way as earlier:
wm_distance <- nb2mat(nb_distance, style='B')
If we we autoprint the object, we will see this is a matrix with 1 and 0 just as the contiguity-based matrix we generated earlier.
Subsubsection 9.2.5.2 Using the inverse distance weights
Instead of just using a critical threshold of distance to define a neighbour, we could use a function to define the weight assigned to the distance between two areas. Typical functions used are the inverse distance, the inverse distance raised to a power, and the negative exponential. All these functions essentially do the same; they decrease the weight as the distance increases. They impose a distance decay effect. Sometimes these functions are combined with a distance threshold [178].
The first steps in the procedure simply reproduce what we did earlier up to the point where we generated the
bur_dist_1 object, so we won’t repeat them here and simply invoke this object. We adapt code and instructions provided by the tutorials of the Center for Spatial Data Science of the University of Chicago prepared by [182]. Once we get to this point, we need to generate the inverse distances. This requires us to take the following steps:
-
Calculate the distance between all neighbours
-
Compute the inverse distance for the calculated distances
-
Assign them as weight values
Let’s go one step at the time. First we use
nbdists() from the spdep package to generate all the distances:
distances <- nbdists(nb_distance, xy)
Now that we have an object with all the distances, we need to compute the inverse distances. Taking the inverse simply involves dividing 1 by each of the distances. We can easily do that with the following code:
inv_distances <- lapply(distances, function(x) (1/x))
If you look inside the generated object, you will see all the inverse distances. You will notice that the values are rather small. The projection we use measures distance in metres. Considering our polygons are census geographies "close" to the ideas of neighbourhoods you will have distances between their centroids that are large. If we take the inverse, as we did, we will unavoidably produce small values very close to zero. We could rescale this to obtain values that are not this small. We can change our units to something larger, like kilometres, or if we want, we can just apply some maths here:
inv_distances_r <- lapply(distances, function(x) (1/(x/100)))
Once we have generated the inverse distance, we can move to generate the weights based on these using
nb2mat().
wm_inverse_dist <- nb2mat(nb_distance,
glist = inv_distances_r,
style='B')
The
glist argument identifies the inverse distances, and the style set to "B" ensures we are not using row standardisation. If we view this object, we will see that rather than zeroes and ones we see the inverse distances used as weights. In this case, all areas neighbour each other (but not themselves!) and their "neighbouringness" is defined with this distance measure.
Subsection 9.2.6 Spatial weight matrix based on k-neighbours
[180] suggest that using the distance-based spatial weight matrix tend to result in too few neighbours for rural areas and many neighbours for urban areas. In this case they point out that the k-nearest neighbourhood structure may be more appropriate. K-neighbours also avoids the problem of isolates without having to specify a radius for a critical threshold.
Computing the list of neighbours is pretty straightforward and uses functions we are familiar with by now. In the first place we use
knearneigh() to get an object of class knn. In this instance, though, we pass a parameter k identifying the number of neighbours we want for each area. We will set this to 3. Then we use knn2nb() to generate the list of neighbours.
nb_k3 <- knn2nb(knearneigh(xy, k = 3))
As before, we can use this list of neighbours to create a matrix and we can visualise the relationships with a plot:
# Get the matrix
wm_k3 <- nb2mat(nb_k3, style='B')
# We first plot the boundaries
plot(cc_burglary_sp, col='white', border='blue', lwd=2)
# Then we plot the distance-based neighbours
plot(nb_k3, xy, lwd=2, col="red", add=TRUE)
