Section 9.3 Choosing the correct matrix
So what is the best approach to creating your spatial weights matrix? As we will see later, the spatial weight matrix not only is used for exploring spatial autocorrelation, it also plays a key role in the context of spatial regression. How we specify the matrix is important.
As noted above, contiguity and distance are the two most commonly used methods. In earlier applications researchers would select a particular matrix without a strong rationale for it. You should be guided by theory in how you define the neighbours in your study; unfortunately "substantive theory regarding spatial effects remain an underdeveloped research area in most of the social sciences" [179] (p. 23) or as [180] put it "there is limited theory available regarding how to select which spatial matrix to use."
In this context, it may make sense to explore various definitions and see if our findings are robust to these specifications. [181] developed an algorithm (AMOEBA) which is a design for the construction of a spatial weights matrix using empirical data that can also simultaneously identify the geometric form of spatial clusters. [180] propose to use this as a possible starting point in the comparisons, though the procedure is not yet fully implemented in R (there is an
AMOEBA package, but it only generates a vector with the spatial clusters, not the matrix). They also suggest using the spatial weight matrix corresponding to high spatial autocorrelation and high statistical significance as more appropriate for exploratory data analysis of a given variable, and for the purpose of regression to use the matrix with the highest correlation and significance when exploring the residuals of the model.
In any case, the choice of spatial weights should be made carefully, as it will affect the results of subsequent analyses based on it. But now, that we know how to build one, we can move to implementing it in order to measure whether or not we can observe spatial autocorrelation in our data.
Subsection 9.3.1 Measuring global autocorrelation
Subsubsection 9.3.1.1 Moran’s I
The most well-known measure of spatial autocorrelation is Moran’s I. It was developed by Patrick Alfred Pierce Moran, an Australian statistician. Like the more familiar Pearson’s \(r\) correlation coefficient, which is used to measure the dependency between a pair of variables, Moran’s I aims to measure dependency. This statistic measures how one area is similar to others surrounding it in a given attribute or variable. When the values of an attribute in areas defined as neighbours in our spatial weight matrix are similar, Moran’s I will be large and positive (closer to 1). When the values of an attribute in neighbouring areas are dissimilar Moran’s I will be negative, tending to -1. In the first case, when there is positive spatial dependence, we speak of clustering. In the latter, what we observe is spatial dispersion (areas alike, in an attribute, "repel" each other). In social science applications, positive spatial autocorrelation is more likely to be expected. There are fewer social processes that would lead us to expect negative spatial autocorrelation. When Moran’s I is close to zero, we fail to find evidence for either dispersion or clustering. The null hypothesis in a test of spatial autocorrelation is that the values of the attribute in question are randomly distributed across space.
Moran’s I is a weighted product-moment correlation coefficient, where the weights reflect geographic proximity. Unlike autocorrelation in timeseries (which is unidimensional: the future cannot affect the past), spatial autocorrelation is multi-dimensional, the dependence between an area and its neighbours is either simultaneous or reciprocal. Moran’s I is defined as:
\begin{equation*}
I = \frac{N}{W} \cdot \frac{\sum_{i} \sum_{j} w_{ij}(x_{i}-\bar{x})(x_{j}-\bar{x})}{\sum_{i}(x_{i}-\bar{x})^2}
\end{equation*}
where \(N\) is the number of spatial areas indexed by \(i\) and \(j\text{;}\) \(x\) is our variable of interest; \(\bar{x}\) is the mean of \(x\text{;}\) \(w_{ij}\) is the spatial weight matrix with zeroes in the diagonal; and \(W\) is the sum of all of \(w_{ij}\text{.}\)
We can easily compute this statistic with the
moran.test() function in spdep package. This function does not take as an input the spatial weight matrix that we create using nb2mat(). Instead it prefers to read this information from a listw class object. This kind of object includes the same information than the matrix, but it does it in a more efficient way that helps with the computation (see [185], p. 230 for details). So the first step to use moran.test() is to turn our matrix into a listw class object.
Before we can use the functions from
spdep to compute the global Moran’s I, we need to create a listw type spatial weights object (instead of the matrix we used above). Let’s use the objects we created using the various criteria. To get the same value as above, we use style='B' to use binary (TRUE/FALSE) weights.
lw_queen <- nb2listw(nb_queen, style='B')
lw_distance <- nb2listw(nb_distance, style='B')
lw_k3 <- nb2listw(nb_k3, style='B')
lw_inverse_d <- nb2listw(nb_distance,
glist = inv_distances_r,
style='B')
moran.test(cc_burglary_sp$burglary, lw_queen)
## ## Moran I test under randomisation ## ## data: cc_burglary_sp$burglary ## weights: lw_queen ## ## Moran I statistic standard deviate = 1.6, ## p-value = 0.06 ## alternative hypothesis: greater ## sample estimates: ## Moran I statistic Expectation Variance ## 0.12033 -0.04545 0.01111
So the Moran’s I here is 0.12, which does not look very large. With this statistic, we have a benchmarking problem. Whereas with standard correlation the value is always restricted to lie within the range of -1 to 1, the range of Moran’s I varies with the W matrix. How can we establish the range for our particular matrix? It has been shown that the maximum and minimum values of \(I\) are the maximum and minimum values of the eigenvalues of \(\frac{(W + W^T)}{2}\) (see [185], p. 232 for details). Don’t worry if you don’t know what an eigenvalue is at this point. We can use this understanding to calculate Moran’s range using the following code:
# Adapted from Brunsdon and Comber (2015)
moran_range <- function(lw) {
wmat <- listw2mat(lw)
return(range(eigen((wmat + t(wmat))/2)$values))
}
moran_range(lw_queen)
## [1] -2.541 5.023
On absolute scale, thus, this does not look like a very large value of spatial autocorrelation.
The second issue is whether this statistic is significant. If we assume a null hypothesis of independence, what is the probability of obtaining a value as extreme as 0.12. The spatial autocorrelation (Global Moran’s I) test is an inferential statistic, which means that the results of the analysis are always interpreted within the context of its null hypothesis. For Global Moran’s I statistic, the null hypothesis states that the attribute being analysed is randomly distributed among the features in your study area; said another way, the spatial process promoting the observed pattern of values is random chance. Imagine that you could pick up the values for the attribute you are analysing and throw them down onto your features, letting each value fall where it may. This process (picking up and throwing down the values) is an example of a random chance spatial process. When the p-value returned by this tool is statistically significant, you can reject the null hypothesis of complete spatial randomness.
The
moran.test() function can compute the variance of Moran’s I under the assumption of normality or under the randomisation assumption (which is the default at least we set the argument randomisation to FALSE). We can also use a Monte Carlo procedure. The way Monte Carlo works is that the values of burglary are randomly assigned to the polygons, and Moran’s I is computed. This is repeated several times to establish a distribution of expected values under the null hypothesis. The observed value of Moran’s I is then compared with the simulated distribution to see how likely it is that the observed values could be considered a random draw.
We use the function
moran.mc() to run a permutation test for Moran’s I statistic calculated by using some number of random permutations of our numeric vector, for the given spatial weighting scheme, to establish the rank of the observed statistic in relation to the simulated values. [125] argue that if we compare the Monte Carlo and analytical variances of \(I\) under randomisation, we typically see few differences, arguably rendering Monte Carlo testing unnecessary.
In
moran.mc() we need to specify our variable of interest (burglary), the listw object we created earlier (ww), and the number of permutations we want to run (here we choose 99).
set.seed(1234)
burg_moranmc_results <- moran.mc(cc_burglary_sp$burglary, lw_queen, nsim=99)
burg_moranmc_results
## ## Monte-Carlo simulation of Moran I ## ## data: cc_burglary_sp$burglary ## weights: lw_queen ## number of simulations + 1: 100 ## ## statistic = 0.12, observed rank = 92, p-value = ## 0.08 ## alternative hypothesis: greater
So, the probability of observing this Moran’s I if the null hypothesis was true is 0.08. This is higher than our alpha level of 0.05. In this case, we cannot conclude that there is significant global spatial autocorrelation. We can use this test to check for our various specifications of the spatial weight matrix:
moran.test(cc_burglary_sp$burglary, lw_distance)
moran.test(cc_burglary_sp$burglary, lw_k3)
moran.test(cc_burglary_sp$burglary, lw_inverse_d)
## ## Moran I test under randomisation ## ## data: cc_burglary_sp$burglary ## weights: lw_distance ## ## Moran I statistic standard deviate = 3, p-value ## = 0.002 ## alternative hypothesis: greater ## sample estimates: ## Moran I statistic Expectation Variance ## 0.28763 -0.04545 0.01266 ## ## Moran I test under randomisation ## ## data: cc_burglary_sp$burglary ## weights: lw_k3 ## ## Moran I statistic standard deviate = 1.9, ## p-value = 0.03 ## alternative hypothesis: greater ## sample estimates: ## Moran I statistic Expectation Variance ## 0.18567 -0.04545 0.01423 ## ## Moran I test under randomisation ## ## data: cc_burglary_sp$burglary ## weights: lw_inverse_d ## ## Moran I statistic standard deviate = 2.4, ## p-value = 0.008 ## alternative hypothesis: greater ## sample estimates: ## Moran I statistic Expectation Variance ## 0.23700 -0.04545 0.01380
As we can see, using the different specifications of neighbours produces results that are more suggestive of dependence. Such results further highlight the importance of picking appropriate and justified spatial weights for your analysis.
We can also make a "Moran scatterplot" to visualize spatial autocorrelation. This type of scatterplot was introduced by [183] as a tool to visualise the relationship between the value of interest and its spatial lag, but also to explore local clusters and areas of non-stationarity (we will discuss this in greater detail in Chapter 10). Note the row standardisation of the weights matrix.
wm_queen_rs <- mat2listw(wm_queen, style='W')
# Checking if rows add up to 1 (they should)
mat <- listw2mat(wm_queen_rs)
#This code simply sums each row to see if they sum to 1
# we are only displaying the first 15 rows in the matrix
apply(mat, 1, sum)[1:15]
## E01005065 E01005066 E01005128 E01005212 E01033653 ## 1 1 1 1 1 ## E01033654 E01033655 E01033656 E01033658 E01033659 ## 1 1 1 1 1 ## E01033661 E01033662 E01033664 E01033667 E01033672 ## 1 1 1 1 1
Now we can plot. Anselin (1996) propose to center the variables around their mean, which is not the default taken by
moran.plot(), so we need to specify this ourselves when plotting by virtue of using scale(). This simply substract the mean from each value, so that the mean of the resulting variable will be 0.
moran.plot(as.vector(scale(cc_burglary_sp$burglary)), wm_queen_rs)
The x-axis represents the values of our burglary variable in each unit (each LSOA) and the y-axis represents a spatial lag of this variable. A spatial lag in this context is simply the average value of the burglary count in the areas that are considered neighbours of each LSOA. So we are plotting the value of burglary against the average value of burglary in the neighbours. The slope of the fitted line represents Moran’s I. And you can see the correlation is almost flat here. We will further explain this scatterplot in Chapter 10. [183] suggests that the effective interpretation of this scatterplot should focus on the extent to which the fitted regression line that we see plotted reflects well the overall pattern of association between the variable of interest and its spatial lag: "the indication of observations that do not follow the overall trend represents useful information on local instability or non-stationarity" (p. 116-117).
As with any correlation measure, you could get positive spatial autocorrelation, which would mean that as you move further to the right in the x-axis, you have higher levels of burglary in the surrounding area. This is what we see here. But the correlation is fairly low and as we saw is not statistically significant. You can also obtain negative spatial autocorrelation. That would mean that areas with high level of crime tend (it’s all about the global average effect!) to be surrounded by areas with low levels of crime. This is clearly not what we see here.
Subsubsection 9.3.1.2 Geary’s C
Another measure of global autocorrelation is Geary’s C. Unlike Moran’s I, which measures similarity, Geary’s C measures dissimilarity. The value is achieved from considering the squared difference for attribute similarity. In the equation
\begin{equation*}
C = \frac{(N-1) \cdot \sum_i\sum_j w_{ij}(x_i - x_j)^2}{2 \cdot \sum_i\sum_j w_{ij} \cdot \sum_i(x_i - \bar{x})^2}
\end{equation*}
the \((x_i - x_j)^2\) measures the dissimilarity between observations \(x_i\) and neighbours \(x_j\text{.}\)
This means that the interpretation for Geary’s C is essentially the inverse of that of Moran’s I. Positive autocorrelation is seen when \(C \lt 1\) (or \(z \lt 0\)), while negative autocorrelation is seen when \(C \gt 1\) (or \(z \gt 0\)). The nice thing about Geary’s C is that it doesn’t need any assumptions about linearity to be made (unlike Moran’s I).
So how do we go about implementing this in R? For obtaining this statistic, we can use
geary.test() function from the spdep package. As with Moran’s I, we specify the variable of interest (cc_burglary_sp$burglary) and the spatial weights matrix (lw_queen):
geary.test(cc_burglary_sp$burglary, lw_queen)
## ## Geary C test under randomisation ## ## data: cc_burglary_sp$burglary ## weights: lw_queen ## ## Geary C statistic standard deviate = -1.7, ## p-value = 1 ## alternative hypothesis: Expectation greater than statistic ## sample estimates: ## Geary C statistic Expectation Variance ## 1.39340 1.00000 0.05395
Note that for Geary’s C, the test assumes that the weights matrix is symmetric. For inherently non-symmetric matrices, such as k-nearest neighbour matrices discussed above, you can make use of the
listw2U() helper function which constructs a weights list object to make the matrix symmetric.
It is very important to understand that global statistics like the spatial autocorrelation (Global Moran’s I, Geary’s C) tool assess the overall pattern and trend of your data. They are most effective when the spatial pattern is consistent across the study area. In other words, you may have clusters (local pockets of autocorrelation), without having clustering (global autocorrelation). This may happen if the signs of the clusters negate each other. In general, Moran’s I is a better measure of global correlation, as Geary’s C is more sensitive to local correlations.
Subsection 9.3.2 The rgeoda and sfweight packages
Traditionally, the generation of a weight matrix was done with functions included in the
spdep program, and this is what we had introduced above. We do want to note, however, that more recently the functionality of GeoDa, a very well-known point and click interface for spatial statistics originally developed by Luc Anselin (a key figure in the field of spatial statistics), has been brought to R via the rgeoda package. This package is fast, though "less flexible when modifications or enhancements are desired" [125]. The nomenclature of the functions is slightly more intuitive, and there are subtle variations in some of the formulas. It works with sf objects, and so may be useful for integrating our analysis to these objects. It is enough now that you are aware of the content of this chapter, and one path to implementing what we learned with the sp package. This should equip you to be able to explore these other packages as well, which will have different ways of achieving these outcomes.
For example, to create queen contiguity weights, you can use the
queen_weights() function from the rgeoda package, and as an argument you only need to pass the name of your polygon object (cc_burglary).
queen_w_geoda <- queen_weights(cc_burglary)
summary(queen_w_geoda)
## name value ## 1 number of observations: 23 ## 2 is symmetric: TRUE ## 3 sparsity: 0.189035916824197 ## 4 # min neighbors: 2 ## 5 # max neighbors: 8 ## 6 # mean neighbors: 4.34782608695652 ## 7 # median neighbors: 4 ## 8 has isolates: FALSE
There are similar functions for other types of weights:
# For example:
# For our k 3 nearest neighbours
knn_weights(cc_burglary, 3)
# For distance
distance_weights(cc_burglary, min_distthreshold(cc_burglary))
More details are available in the documentation for
rgeoda.
It is also worth mentioning
sfweight. The sfweight package, still only available from GitHub at the time of this writing, and therefore we could consider it as "an ongoing development" that aims "to provide a more streamlined interface to the spdep" package. The goal of this package is to create a simpler workflow for creating neighbours, spatial weights, and spatially lagged variables. You can read more about it in its GitHub repository.
