Local Moran’s I and Moran scatterplot

Section 10.4 Local Moran’s I and Moran scatterplot

Subsection 10.4.1 Computing the local Moran’s I

The Local Getis and Ord statistic is a local indicator of spatial association. In Chapter 9 of this book, we looked at using Moran’s I as a global measure of spatial autocorrelation. There is also a local version of Moran’s I, which allows for the decomposition of global indicators of spatial autocorrelation such as Moran’s I discussed in Chapter 9, into the contribution of each observation. This allows us to elaborate on the location and nature of the clustering within our study area, once we have found that autocorrelation occurs globally.

🔗

Let’s first look at Moran’s scatterplot for our data. Remember from last chapter that Moran’s scatterplot plots the values of the variable of interest, in this case the burglary count on each micro-grid, against the spatial lagged value. For creating Moran’s scatterplot we generally use row standardisation in the weight matrix so that the y and the x-axes are more comparable.

🔗

# create row standardised weights
lw_queen_rs <-  nb2listw(nb_queen_cc, style='W')

# build Moran's scatterplot
m_plot_burglary <- moran.plot(as.vector(scale(grid_cc$layer)),
                              lw_queen_rs, cex=1, pch=".",
                              xlab="burglary count",  ylab="lagged burglary count")

A scatterplot with a fitted line, with the vertical axis labeled as ’lagged burglary count’, ranging from negative point five to two. The horizontal axis is labeled burglary count, and ranges from negative point five to a little over eight. One dashed vertical and one horizontal lines appear at each zero value, which the fitted line intersects. The plotted points form vertical lines, with most left of the zero dashed line, with values lower than one. More points appear to the right of the dashed line, but end near two. In addition, numerous crossed diamonds with numeric labels appear outside of this area. Many are left of the vertical zero dashed line above the value one, while another concentration is close to the horizontal zero dashed line between two and four. More appear in more extreme locations. — Figure 10.4.1. Moran’s scatterplot for burglary in Manchester
🔗

Notice how the plot is split into four quadrants structured around the centered mean for the value of burglary and its spatial lag. The top right corner belongs to areas that have high level of burglary and are surrounded by other areas that have above the average level of burglary. These are high-high locations. The bottom left corner belongs to the low-low areas. These are areas with low level of burglary and surrounded by areas with below average levels of burglary. Both the high-high and low-low represent clusters. A high-high cluster is what you may refer to as a hot spot. And the low-low clusters represent cold spots.

🔗

In the opposite diagonal we have spatial outliers. They are not outliers in the standard statistical sense, as in extreme observations. Instead, they are outliers in that they are surrounded by areas that are very unlike them. So you could have high-low spatial outliers, areas with high levels of burglary and low levels of surrounding burglary, or low-high spatial outliers, areas that have themselves low levels of burglary (or whatever else we are mapping) and that are surrounded by areas with above average levels of burglary. This would suggest pockets of non-stationarity. [Anselin (1996)] warns these can sometimes surface as a consequence of problems with the specification of the weight matrix.

🔗

The slope of the line in the scatterplot gives you a measure of the global spatial autocorrelation. As we noted earlier, the local Moran’s I helps you to identify the locations that weigh more heavily in the computation of the global Moran’s I. The object generated by the moran.plot() function includes a measure of this leverage, a hat value, for each location in your study area and we can plot these values if we want to geographically visualise those locations:

🔗

# We extract the influence measure and place it in the sp object
grid_cc_sf$lm_hat <- m_plot_burglary$hat
# Then we can plot this measure of influence
tm_shape(grid_cc_sf) + tm_fill("lm_hat")

To compute the local Moran’s I, we can use the localmoran() function from the spdep package.

🔗

locm_burglary <- localmoran(grid_cc_sf$layer, lw_queen_rs, alternative = "two.sided")

As with the local \(G^*\) it is necessary to adjust for multiple comparisons. The unadjusted p-value for the test is stored in the fifth column of the object of class localmoran that we generated with the localmoran() function. As before we can extract this element and pass it as an argument to the p.adjust() function in order to obtain p-values with some method of adjustment for the problem of multiple comparisons.

🔗

p_values_lm <- locm_burglary[,5]
sum(p_values_lm < 0.05)

## [1] 54

NOTE: It is important to note that by default, the localmoran() function does not adjust for the problem of multiple comparisons. This needs to be kept in mind, as it might lead us to erroneously conclude we have many hot or cold spots.

🔗

In our example, if we do not adjust for multiple comparisons, then, there would be 33 areas with a significant p value. If we apply again the Bonferroni correction, these get reduced to only 11.

🔗

bonferroni_lm <- p.adjust(p_values_lm, "bonferroni")
sum(bonferroni_lm < 0.05)

## [1] 10

Another approach is to use conditional permutation for inference purposes, instead of the analytical methods used above. In this case we would need to invoke the localmoran_perm() function also from the spdep package. This approach is based on simulations, so we also have to specify an additional argument nsim setting the number of simulations. [Bivand and Wong (2018)] suggest that "increasing the number of draws beyond 999 has no effect" (p. 741). The iseed argument ensures we use the same random seed and get reproducible results. As before adjusting for multiple comparisons notably reduces the number of significant clusters.

🔗

locm_burglary_perm <- localmoran_perm(grid_cc_sf$layer,
                                      lw_queen_rs,
                                      nsim=499,
                                      alternative="two.sided",
                                      iseed=1)
p_values_lmp <- locm_burglary_perm[,5]
sum(p_values_lmp < 0.05)
bonferroni_lmp <- p.adjust(p_values_lmp, "bonferroni")
sum(bonferroni_lmp < 0.05)

## [1] 54
## [1] 9

We go from 54 to 9 significant regions! This really illustrates the perils of not correcting for multiple comparisons!

🔗

Subsection 10.4.2 Creating a LISA map

In order to map these Local Indicators of Spatial Association (LISA), we need to do some prep work.

🔗

First we must scale the variable of interest. As noted earlier, when we scale burglary, what we are doing is re-scaling the values so that the mean is zero. We use the scale() function, which is a generic function whose default method centers and/or scales the variable.

🔗

# Scale the count of burglary
grid_cc_sf$s_burglary <- scale(grid_cc_sf$layer)
# Creates p values for Bonferroni adjustment
grid_cc_sf$localmp_bonf <- p.adjust(p_values_lmp, "bonferroni")

To produce the LISA maps, we also need to generate a spatial lag: the average value of the burglary count in the areas that are considered neighbours of each LSOA. For this we need our listw object, which is the object created earlier, when we generated the list with weights using row standardisation. We then pass this listw object into the lag.listw() function, which computes the spatial lag of a numeric vector using a listw sparse representation of a spatial weights matrix.

🔗

#create a spatial lag variable and save it to a new column
grid_cc_sf$lag_s_burglary <- lag.listw(lw_queen_rs, grid_cc_sf$s_burglary)

Make sure to check the summaries to ensure the numbers are in line with our expectations, and nothing weird is going on.

🔗

summary(grid_cc_sf$s_burglary)
summary(grid_cc_sf$lag_s_burglary)

We are now going to create a new variable to identify the quadrant in which each observation falls within Moran’s scatterplot, so that we can tell apart the high-high, low-low, high-low, and low-high areas. We will only identify those that are significant according to the p-value that was provided by the local Moran function adjusted for multiple comparisons with our Bonferroni adjustment. The data we need for each observation, in order to identify whether it belongs to the high-high, low-low, high-low, or low-high quadrants are the standardised burglary score (positive or negative), the spatial lag score, and the p-value.

🔗

Essentially all we’ll be doing is assigning a variable values based on where in the plot it is. So for example, if it’s in the upper right, it is high-high and has values larger than 0 for both the burglary and the spatial lag values. If it’s in the upper left, it’s low-high and has a value larger than 0 for the spatial lag value, but lower than 0 on the burglary value. And so on, and so on.

🔗

A scatterplot divided into four quadrants, the vertical axis being ’Lagged X’ from negative two to two, and horizontal axis ’X’ from negative two to two. The four quadrants are shaded in blue and green, and each has a white box with a label indicating whether each ’Lagged X’ and ’Standardized X’ are greater than or less than zero. A handful of crossed diamonds are plotted, spreading from the top left to bottom right, with most of them lying in those two quadrants. — Figure 10.4.2. Quadrants of a Moran scatterplot
🔗

So let’s first initialise this variable. In this instance we are creating a new column in the sp object and calling it "quad_sig".

🔗

grid_cc_sf <- grid_cc_sf %>%
  mutate(quad_sig = case_when(localmp_bonf >= 0.05 ~ "Not significant",
                              s_burglary[,1] > 0 &
                                lag_s_burglary[,1] > 0 &
                                localmp_bonf < 0.05 ~ "High-high",
                              s_burglary[,1] < 0 &
                                lag_s_burglary[,1] < 0 &
                                localmp_bonf < 0.05 ~ "Low-low",
                              s_burglary[,1] < 0 &
                                lag_s_burglary[,1] > 0 &
                                localmp_bonf < 0.05 ~ "Low-high",
                              s_burglary[,1] > 0 &
                                lag_s_burglary[,1] < 0 &
                                localmp_bonf < 0.05 ~ "High-low" ))

Now we can have a look at what this returns:

🔗

table(grid_cc_sf$quad_sig)

## 
##       High-high        Low-high Not significant 
##               5               4             513

So the 9 significant clusters we found split into 5 high-high and 4 low-high.

🔗

Let’s put them on a map, using the standard colours used in this kind of maps:

🔗

tm_shape(grid_cc_sf) +
  tm_fill("quad_sig",
          palette= c("red","blue","white"),
          labels = c("High-High","Low-High", "non-significant"),
          alpha=0.5) +
  tm_borders(alpha=.5) +
  tm_layout(frame = FALSE,
            legend.position = c("right", "bottom"),
            legend.title.size = 0.8,
            legend.text.size = 0.5,
            main.title = "LISA with the p-values",
            main.title.position = "centre",
            main.title.size = 1.2)

There we have it, our LISA map of a grid overlay over Manchester’s City Centre ward, identifying local spatial autocorrelation in burglaries.

🔗

Prev Top Next