Mapping And Modelling Geographic Data In R

Session 4: Spatial weights and neighbours

Introduction

This session demonstrates some of R’s capability for spatial analysis by showing how to create spatial weightings based on contiguity, nearest neighbours and by distance. Methods of inverse distance weighting are also introduced and used to produce Moran plots and tests of spatial dependency in data, including the results of regression analysis.

To get started:

if(!"sp" %in% installed.packages()) install.packages("sp")
if(!"sf" %in% installed.packages()) install.packages("sf")
rm(list=ls()) # Deletes everything in the workspace - use with care!
require(sp)
load("spatial-weights.RData")
ls()

This will load an R workspace containing two objects, combined.map and districts. The first is if class,

class(combined.map)

and contains (synthetic) information about land parcels in Beijing. Each land parcel is represented by a point centroid. The second is of class

class(districts)

and contains some information about districts within Beijing. The sp objects can be converted to spatial features using

require(sf)
st_as_sf(combined.map) -> sf_combined.map
st_as_sf(districts) -> sf_districts

and we can take a glimpse of the data:

if(!"tidyverse" %in% installed.packages()) install.packages("tidyverse")
require(tidyverse)
sf_combined.map %>%
  glimpse()
sf_districts %>%
  glimpse()

The variables are as follow:

LNPRICE: Log of the land price: RMB per square metre LNAREA: Log of the land parcel size DCBD: distance to the CBD on a log scale DELE: distance to the nearest elementary school on a log scale DRIVER: distance to the nearest river on a log scale DPARK: distance to the nearest park on a log scale Y0405: A dummy variable (1 or 0) indicating whether the land was sold in the period 2004-5 Y0607: A dummy variable indicating whether the land was sold in the period 2006-7 Y0809: A dummy variable indicating whether the land was sold in the period 2007-8 (All of the remaining land parcels were sold before 2004) SP_ID: An ID for the polygon (district) in which the land parcel is located POPDEN: the population density in each district. data source: the fifth census data in 2000. JOBDEN the jobden density in each district. data source: the fifth census data in 2000.

Identifying Contiguous neighbours

A contiguity matrix is one that identifies polygons that share boundaries and (in what is called the Queen’s case) corners too. In other words, it identifies neighbouring areas. To do this we use the poly2nb(…) function – polygons to an object of class neighbours - in the spdep (spatial dependency) library. A current problem is that spdep expects sp objects not sf ones! Undoubtedly this will change but, for now, the conversion of the sp object to an sf one was premature so we shall use the originals.

if(!"spdep" %in% installed.packages()) install.packages("spdep")
require(spdep)
contig <- poly2nb(districts)
summary(contig)

The summary of the neighbours object (which is of class nb) shows that there are 134 regions (districts) – the same number as when we glimpsed the data, above – with each being linked to 5.46 others, on average. There are two regions with no links and two regions with 10 links. The Queen’s case is assumed by default, see ?poly2nb.

It is helpful to learn a little more about the structure of the contiguity object. It is an object of class nb which is itself a type of list. The function str() is sometimes useful to take a look at the structure of an object in R.

class(contig)
typeof(contig)
str(contig)

If we look at the first item in the list, we obtain the following:

contig[[1]]

What this is showing is that district 1 has two contiguous neighbours, which are the 52nd and 54th entries in the list. Reciprocally,

contig[[52]]

contig[[54]]

which both include 1 (as they should – the relationship is symmetric: if \(a\) shares a border with \(b\) then \(b\) shares a border with \(a\)).

Going back to the summary above, region 20 has no links. However,

contig[[20]]

suggests that it actually has four. The confusion here is that polygon 20 in the original sp object (and the shapefile before that) is actually the 21st entry in the list. We can get a sense of this if we look at,

head(districts$SP_ID)

and note that the numbering begins with zero, not one. Note also that the tidyverse is not implemented with sp classes but is with sf. Hence, the following generates an error,

districts %>%
  dplyr::select(SP_ID) %>%
  print(n = 3)

but the following – with a conversion to sf - works fine:

districts %>%
  st_as_sf(.) %>%
  dplyr::select(SP_ID) %>%
  print(n = 3)

In any case, if region 20 is actually the 21st in the list then the following should have zero neighbours,

contig[[21]]

… and it does! We can confirm this by plotting the source data,

plot(districts[-21,], col="red")
plot(districts[21,], col="blue", add=T)

Or, using the simple feature version of the object that was created earlier,

plot(districts[-21,0], col="red")
plot(districts[21,0], col="blue", add=T)

• If you are eagle-eyed, you may notice that the code block above includes a zero (e.g. [-21,0]) where as the code above does not (e.g. [21,]). Try deleting the 0 in the code block above and see what happens!

K nearest neighbours (knn)

It does not make sense sense to evaluate contiguity of point data because points are zero dimensional (and have no borders). Instead, we could find, for example, the six nearest neighbours to each point. There are various ways to do this, including using the RANN package:

if(!"RANN" %in% installed.packages()) install.packages("RANN")
require(RANN)
pts <- coordinates(combined.map) 
knn6 <- nn2(pts, pts, k = 6)

Looking at the attributes of the object knn6, it shows that nn2 creates two outputs, one with name nn.idx, the other with name nn.dists.

attributes(knn6)
str(knn6)

If we are interested in identifying the six nearest neighbours to point 1 then they are,

knn6$nn.idx[1,]

and their distances from the point are

knn6$nn.dists[1,]

Note that the six nearest neighbours to point 1 include point 1 (and the six nearest neighbours to point 2 include point 2, etc.). This isn’t necessarily what is wanted and a problem that spdep’s knearneigh() function avoids. It also can work directly on the sp object:

knear6 <- knearneigh(combined.map, k=6)
attributes(knear6)
knear6$nn[1,]

However, unlike RANN (which, actually, it is using), it doesn’t retain the distances to the nearest neighbours.

Spatial dependency exists in a set of data if the value at one location is not independent of values at surrounding locations – if there is a geographical pattern in the data. Imagine we are interested in calculating the correlation between some variable (call it \(x\)) at each of the points and at each of the points’ closest neighbour. Using the simulated data about the price of land parcels in Beijing, we can calculate this correlation as follows:

x <- combined.map$LNPRICE
cor(x, x[knear6$nn[,1]])

Similarly, the correlation with the second nearest neighbour is

x <- combined.map$LNPRICE
cor(x, x[knear6$nn[,2]])

and the correlations from the first to sixth neighbour are

r <- sapply(1:6, function(i) {
  cor(x, x[knear6$nn[,i]])
})
df <- data.frame(k = 1:6, r = r)
print(df)

(All this is a function where the variable \(i\) takes the values from 1 to 6, in turn, so, in the first instance, it is calculating, cor(x, x[knear6$nn[,1]]) then cor(x, x[knear6$nn[,2]]) and so forth up to cor(x, x[knear6$nn[,6]])

Plotting these,

if(!"ggplot2" %in% installed.packages()) install.packages("ggplot2")
require(ggplot2)
df %>%
  ggplot(aes(x = k, y = r)) +
  geom_line() +
  xlab("kth nearest neighbour") +
  ylab("Correlation, r")

What these values suggest is that even at the sixth nearest neighbour, the value of a land parcel at any given point tends to be similar to the value of the land parcels around it –- an example of positive spatial autocorrelation. However, the level of similarity is highest for the nearest neighbour and decreases the further the neighbour is away.

An issue is that the threshold of six nearest neighbours is purely arbitrary. An interesting question is how far – how many neighbours away – we can typically go from a point and still find a similarity in the land price values. One way to determine this would be to carry on with the calculations above, repeating the procedure until we get to, say, the 250th nearest neighbour.

knear250 <- knearneigh(combined.map, k=250)
r <- sapply(1:250, function(i) {
  cor(x, x[knear250$nn[,i]])
})
data.frame(k = 1:250, r = r) %>%
  ggplot(aes(x = k, y = r)) +
  geom_line() +
  geom_smooth(se = FALSE) +
  xlab("kth nearest neighbour") +
  ylab("Correlation, r")

Looking at the plot we find that the land prices become more dissimilar (less correlated) the further away we go from each point, dropping to zero correlation from after about the 200th nearest neighbour. The rate of decrease in the correlation is greatest to about the 35th neighbour, after which it begins to flatten.

We can also determine the p-values associated with each of these correlations and identify which are not significant at a 99% confidence:

pvals <- sapply(1:250, function(i) {
  cor.test(x, x[knear250$nn[,i]])$p.value
})
which(pvals > 0.01)

It is from about the 100th neighbour that the correlations begin to become regularly insignificant. Whether this is usefully information or not is a moot point: a measure of statistical significance is really only an indirect measure of the sample size. It may be better to make a decision about the threshold at which the neighbours are not substantively correlated based on the actual correlations (the effect sizes) rather than their p-values. Whilst it remains a subjective choice, here we will here use the 35th neighbour as the limit, before which the correlations are typically equal to r = 0.20 or greater.

knear35 <- knearneigh(combined.map, k=35)
class(knear35)

The object is of class knn. To convert it to class nb (the same as the contiguity object from earlier),

knear35nb <- knn2nb(knear35)
knear35nb

Note that this is a non-symmetric list – just because \(a\) is one of the 35 nearest neighbours to \(b\), it does not follow that \(b\) must be one of the 35 nearest neighbours to \(a\).

Identifying neighbours by (physical) distance apart

It is also possible to identify the neighbours of points by their Euclidean or Great Circle distance apart using the function dnearneigh(). For example, if we wanted to identify all points between 100 and 1000 metres of each other:

d100to1000 <- dnearneigh(combined.map, 100, 1000)
d100to1000

• On average, how many regions are between 100 and 1000 metres of another?

Creating a Spatial Weights matrix

What we have done so far is use various methods to create a list of neighbours with flexibility to decide about what counts as a neighbour. The next stage is to convert that list into a spatial weights matrix so we can use it for various methods of spatial analysis. This extra stage of conversion may seem like an unnecessary additional chore. However, the creation of the spatial weights matrix allows us to define the strength of relationship between neighbours. For example, we may want to give more weight to neighbours that are located closer together and less weight to those that are further apart (decreasing to zero beyond a certain threshold). The nb object defines which are neighbours. Converting it to a spatial weights object says how strongly they are neighbours.

Creating a binary list of weights

We could create a simple binary ‘matrix’ from any of our existing lists of neighbours. In principle:

spcontig <- nb2listw(contig, style="B")

Note, however, the error message, which arises because two of the Chinese districts do not share a boundary with others. Correcting for this,

spcontig <- nb2listw(contig, style="B", zero.policy=T)

The same problem does not arise for the one hundred nearest neighbours (by definition, it cannot – each point has neighbours) but it does for the distance based list if empty sets are not permitted (they are permitted by setting zero.policy = T):

spknear35 <- nb2listw(knear35nb, style="B")
spd100to1000 <- nb2listw(d100to1000, style="B", zero.policy=T)

What we create are objects of class listw (a list of spatial weights). Looking at the first of these objects we can see how it has been constructed. It contains binary weights (style B); district 1 has two neighbours, districts 52 and 54; and both of those have been given a weight of one (all other districts therefore have a weight of zero with district 1). Similarly, district 2 has neighbours 3, 4, 6, 999 and 100, each with a weight of one.

class(spcontig)
typeof(spcontig)
names(spcontig)
spcontig$style
spcontig$neighbours[[1]]
spcontig$weights[[1]]
spcontig$neighbours[[2]]
spcontig$weights[[2]]

Creating a row-standardised list of weights

Using binary weights where (1) indicates two places are neighbours, and (0) indicates they are not, may create a problem when different places have different numbers of neighbours (as is the case for both the contiguity and distance-based approaches). Imagine a calculation where the result is in some way dependent upon the sum of the weights involved. All things being equal, we expect places with more neighbours to generate larger values of y simply because they have more non-zero values contributing to the sum. A way around this problem is to scale the weights so that for any one place they sum to one – a process known as row-standardisation and actually the default option:

spcontig <- nb2listw(contig, zero.policy=T)
names(spcontig)
spcontig$style
spcontig$neighbours[[1]]
spcontig$weights[[1]]
spcontig$neighbours[[2]]
spcontig$weights[[2]]

In the case of the contiguity matrix, district 1 still has neighbours 52 and 54, and district 2 still has neighbours 3, 4, 6, 99 and 100 but the weights are now row-standardised (style W) and in each case they sum to one.

Creating an inverse distance weighting (IDW)

A more ambitious undertaking is to decrease the weighting given to two points according to their distance apart, reducing to zero, for example, beyond the 35th nearest neighbour. To achieve this, we can use the RANN package from earlier to find the distances (another possibility is spDists – see ?spDists) then assign the weights based on the inverse of those distances.

pts <- coordinates(combined.map) # Get the centroid of each district
knn <- nn2(pts, pts, k = 36)
# Recall that using RANN includes each district as its own neighbour
d <- knn$nn.dists   # The distances
d <- d[,-1]         # Delete the first column of zero distances
idw <- 1/d          # the inverse of distance 
glist <- lapply(1:nrow(idw), function(i) idw[i,])
# This is just a way of converting the inverse distances into weights
knear <- knearneigh(as(combined.map, "SpatialPoints"), k=35)
# For compariability, this is also using the centroid of each district
knearnb <- knn2nb(knear)  # Create a neighbourhood list
spknear35IDW <- nb2listw(knearnb, glist=glist)
# Combine the neighbourhood list with the inverse distance weighting

That isn’t the most efficient way of doing it but it appears to more-or-less work because if we look at the weights assigned to the 35 neighbours of the first district they are decreasing with distance:

spknear35IDW$weights[[1]]

Why more-or-less? Because the weights are row-standardised! This means that the inverse distance weighting is a function of the local distribution of points around each point not just how far away they are. For example, consider a point where all its neighbours are quite far from it. Using a strictly distance-based weighting each of those neighbours should receive a low weighting. However, once row standardisation is applied those low weights will be scaled upwards to sum to one. Reciprocally, imagine a point where all its neighbours are very close. Using a distance-based weighting those neighbours should receive a relatively high weighting; in effect, though, they will be scaled downwards by the row standardisation. This may sound undesirable and counter to the objectives of inverse distance weighting, and can be prevented by changing the weights style,

spknear35IDW <- nb2listw(knearnb, glist=glist, style = "C")
spknear35IDW$weights[[1]]

However, imagine the points sample across both urban and rural areas. The distances between points will most likely be smaller in the urban regions (where the density of points is greater, reflecting the greater population density), with greater distances between points in the rural regions. If row standardisation is not applied then the net result will be to give more weight to the urban parts of the region such that any subsequent calculation dependent upon the sum of the weights will be more strongly influenced by the urban areas than by the rural ones. Therefore careful though needs to be given to the style of weights to use.

Variants of IDW

Common forms of inverse distance weighting include the bisquare and Gaussian functions. Adapting the code above, a bisquare version is

pts <- coordinates(combined.map)
knn <- nn2(pts, pts, k = 36)
d <- knn$nn.dists
d <- d[,-1]         
glist <- lapply(1:nrow(d), function(i) {
  max.d <- d[i,35]
  1 - (d[i,] / max.d)
  # The above calculates the bisquare weighting
})
knear <- knearneigh(as(combined.map, "SpatialPoints"), k=35)
knearnb <- knn2nb(knear)
spknear35bisq <- nb2listw(knearnb, glist=glist, style="C")

… and a Gaussian version is

pts <- coordinates(combined.map)
knn <- nn2(pts, pts, k = 36)
d <- knn$nn.dists
d <- d[,-1]         
glist <- lapply(1:nrow(d), function(i) {
  max.d <- d[i,35]
  exp(-0.5*d[i,]^2 / max.d^2)   
  # The above calculates the bisquare weighting
})
knear <- knearneigh(as(combined.map, "SpatialPoints"), k=35)
knearnb <- knn2nb(knear)
spknear35gaus <- nb2listw(knearnb, glist=glist, style="C")

A useful illustration of the types of weighting and how they differ can be found as Figure 1 of GWmodel: An R Package for Exploring Spatial Heterogeneity Using Geographically Weighted Models. Generally the shape (the rate of decay) matters less than the number of neighbours that are included.

Using Spatial Weights

Creating a spatially lagged variable

Once we have the spatial weights we can use them to create a spatially lagged variable. For example, if \(x_i\) is the value of the land parcel at point \(i\), then its spatial lag is the mean value of the land parcels that are the neighbours of \(i\), where those neighbours are defined by the spatial weights. More precisely, it is the weighted mean value if, for example, inverse distance weighting has been employed. It is straightforward to calculate the spatially lagged variable. For example,

x <- combined.map$LNPRICE
lagx <- lag.listw(spknear35gaus, x)

Having done so, the correlation between points and their neighbours can be calculated,

cor.test(x, lagx)

Here there is evidence of significant positive spatial autocorrelation -– that the land price at one point tends to be similar to the land prices of its neighbours. This can be seen if we plot the two variables on a scatter plot, although the relationship is also somewhat noisy and may not be linear.

data.frame(x = x, lagx = lagx) %>%
  ggplot(aes(x=x, y=lagx)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", se = FALSE) +
  geom_vline(xintercept = mean(x), linetype = "dotted") +
  geom_hline(yintercept = mean(lagx), linetype = "dotted")

A Moran plot and test

The plot created above is known as a Moran plot. A more direct way of producing it is to use the moran.plot() function,

moran.plot(x, spknear35gaus)

which flags potential outliers / influential observations. To suppress their labelling, include the argument labels=F – try it!

The Moran coefficient and related test provide a measure of the spatial autocorrelation in the data, given the spatial weightings.

moran.test(x, spknear35gaus)

Essentially the Moran statistic is a correlation value, although it need not be exactly zero in the presence of no correlation (here the expected value is not zero but slightly negative) and need not be in the range -1 to +1. The interpretation though is that the price of the land parcels and their neighbours are positively correlated: there is a tendency for like-near-like values.

More strictly, we should acknowledge the note found under ?moran.test() that the derivation of the tests assumes the matrix is symmetric, which spknear35gaus is not. To correct for this we can use a helper function, listw2U(),

moran.test(x, listw2U(spknear35gaus))

What is the range of the Moran’s value?

XXXX observe that the potential range for a Moran’s value is dependent upon the spatial weights matrix. Using their code:

moran.range <- function(lw) {
  wmat <- listw2mat(lw)
  return(range(eigen((wmat + t(wmat))/2)$values))
}
moran.range(spknear35gaus)
moran.range(spknear35bisq)
moran.range(spcontig)
moran.range(nb2listw(contig, zero.policy = TRUE))

None of these inspire confidence that Moran values typically lie in the range -1 to +1 (either that or the method to determine the extremes is wrong but I have no evidence that is the case). Increasingly I wonder if it would be easier to use, e.g.

x <- combined.map$LNPRICE
r <- cor(x, lag.listw(spknear35gaus, x))
r

… which is just the Pearson correlation between the variable and its lagged equivalent and will lie in the range from -1 to +1. Calculating the significance of the value can be done by simulation. For example,

listw <- spknear35gaus
x <- combined.map$LNPRICE
nsims <- 1000

vals <- sapply(1:nsims, function(i) {
  xx <- sample(x, length(x), replace = FALSE)  
  cor(xx, lag.listw(listw, xx, zero.policy = TRUE))
})

… calculates the correlation 1000 times with random re-arrangements of the data (where the same set of measurements is used, the neighbours remain the same as does their weighting, but the data are shuffled around geographically). If the actual correlation, \(r\), is greater than 99% of the randomly generated values then we have reason to believe that the positive spatial autocorrelation between the values and their neighbours (as defined using the particular spatial weights matrix) is ‘significant’ at above a 99% level.

r > quantile(vals, probs = 0.99)

Checking a regression model for spatially correlated errors

Assume we fit the regression model,

model1 <- lm(LNPRICE ~ DCBD + DELE + DRIVER + DPARK + POPDEN + JOBDEN + Y0405 + Y0607 + Y0809, data=combined.map)

An assumption of regression is that the residuals are supposed to be left over and random ‘noise’. If that is true, they should be independent of one another and not display any spatial pattern. Is that the case here? One way to check is by mapping them:

resids <- rstudent(model1)
resids <- cut(resids, breaks = c(-Inf, -2.58, -1, 1, 2.58, Inf))
combined.map$resids <- resids
require(tmap)
tmap_style("col_blind")
tm_shape(districts) +
  tm_polygons() +
tm_shape(combined.map) +
  tm_bubbles("resids", size = 0.05)

I would say that it is hard to judge. An alternative is to use a Moran’s test, ideally the version of the test intended for use on regression residuals.

lm.morantest(model1, listw2U(spknear35gaus))
# listw2U is used because the weights are not symmetrical

The estimated correlation is about 0.097. Not huge, perhaps but significant enough to question the assumption of independence. Note that this result is dependent on the spatial weightings. If we change them, then the results of the Moran test will change also. For example, using the bi-square weightings,

lm.morantest(model1, listw2U(spknear35bisq))

Here the change is slight, largely because the rate of decay of the inverse distance weighting matters less than the number of neighbours it decays to.

Localised Moran values

The Moran’s value is an example of a global statistic, which is one calculated for the whole study region. It is a single (an average) measure of spatial autocorrelation between neighbours. It can, however, by ‘unpacked’ to give a series of local Moran values – one for each location in the data:

localm <- localmoran(combined.map$LNPRICE, listw2U(spknear35gaus))
head(localm)

This produce a matrix, where the first column contains the local Moran values. A statistical test of significance also is included.

A positive local Moran value indicates a local cluster of alike values; a negative Moran value indicates a cluster of dissimilar values. What they do not distinguish between is whether, for example, high land price values are surrounded by other high land prices or whether low land prices are surrounded by other low land prices – both are examples of positive spatial autocorrelation (similar values clustering together). Equally, they do not distinguish between whether low prices are surrounded by high prices or high prices by low (both are examples of negative spatial autocorrelation).

To distinguish them, we need to refer to the original data (let’s call it \(x\)). If \(x\) is above the mean and the local Moran value is positive and significant then we might argue that it is a ‘hot spot’ of a high value surrounded by other high values: high-high. If \(x\) is below the mean and the local Moran value is positive and significant then we could say it is a ‘cold spot’ of a low value surrounded by other low values: low-low. If \(x\) is above the mean but the Moran value is negative and significant then high-low; if \(x\) is below the mean but the Moran value is negative and significant then low-high:

x <- combined.map$LNPRICE
localm <- localmoran(combined.map$LNPRICE, listw2U(spknear35gaus))
type <- rep(NA, times = nrow(combined.map))
type <- ifelse(x > mean(x) & localm[,1] > 0 & localm[,1] < 0.05,
       "High - High", type)
type <- ifelse(x > mean(x) & localm[,1] < 0 & localm[,1] < 0.05,
       "High - Low", type)
type <- ifelse(x < mean(x) & localm[,1] < 0 & localm[,1] < 0.05,
       "Low - Low", type)
type <- ifelse(x < mean(x) & localm[,1] > 0 & localm[,1] < 0.05,
       "Low - High", type)
combined.map$type <- type
sub.map <- combined.map[!is.na(type),]
tmap_style("cobalt")
tm_shape(districts) +
  tm_polygons() +
tm_shape(sub.map) +
  tm_bubbles("type", size = 0.05)

Moran values can also be calculated for regression residuals. For example,

localmr <- localmoran.exact(model1, nb = knn2nb(knearneigh(combined.map, k = 35)))

See ?localmoran.exact for further details.

Discussion and Conclusion

A spatial weights matrix is a way of encoding the relationship between neighbours and is a marker of spatial analysis. The difficulty is that the results obtained are a function of the spatial weights matrix used – change the matrix and you will changes the results too. How, then, to choose the right matrix? In practice, this is often given too little thought but, ideally, it should be based on some theoretical basis or through an analysis of the scale of the geographical patterning of the data. Most likely there will not be one correct choice so an element of trial-and-error and sensitivity analysis (how sensitive are the results to a change in the matrix?) may be required.

If you have time left over, read the Bivand article from the further reading below.

Further Reading

Bivand R, 2019. Creating Neighbours