pkgs <- c("sf", "mapview", "spdep", "spatialreg", "tmap", "viridisLite") # note: load spdep first, then spatialreg
lapply(pkgs, require, character.only = TRUE)3 Spatial Relationships W
Required packages
Session info
sessionInfo()R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)
Matrix products: default
locale:
[1] LC_COLLATE=English_United Kingdom.utf8
[2] LC_CTYPE=English_United Kingdom.utf8
[3] LC_MONETARY=English_United Kingdom.utf8
[4] LC_NUMERIC=C
[5] LC_TIME=English_United Kingdom.utf8
time zone: Europe/London
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
[1] viridisLite_0.4.2 tmap_3.3-3 spatialreg_1.2-9
[4] Matrix_1.5-4.1 spdep_1.2-8 spData_2.2.2
[7] mapview_2.11.0 sf_1.0-13
loaded via a namespace (and not attached):
[1] xfun_0.39 raster_3.6-20 htmlwidgets_1.6.2
[4] lattice_0.21-8 vctrs_0.6.3 tools_4.3.1
[7] crosstalk_1.2.0 LearnBayes_2.15.1 generics_0.1.3
[10] parallel_4.3.1 stats4_4.3.1 sandwich_3.0-2
[13] tibble_3.2.1 proxy_0.4-27 fansi_1.0.4
[16] pkgconfig_2.0.3 KernSmooth_2.23-21 satellite_1.0.4
[19] RColorBrewer_1.1-3 leaflet_2.1.2 webshot_0.5.4
[22] lifecycle_1.0.3 compiler_4.3.1 deldir_1.0-9
[25] munsell_0.5.0 terra_1.7-29 leafsync_0.1.0
[28] codetools_0.2-19 stars_0.6-1 htmltools_0.5.5
[31] class_7.3-22 pillar_1.9.0 MASS_7.3-60
[34] classInt_0.4-9 lwgeom_0.2-13 wk_0.7.3
[37] abind_1.4-5 boot_1.3-28.1 multcomp_1.4-24
[40] nlme_3.1-162 tidyselect_1.2.0 digest_0.6.31
[43] mvtnorm_1.2-2 dplyr_1.1.2 splines_4.3.1
[46] fastmap_1.1.1 grid_4.3.1 colorspace_2.1-0
[49] expm_0.999-7 cli_3.6.1 magrittr_2.0.3
[52] base64enc_0.1-3 dichromat_2.0-0.1 XML_3.99-0.14
[55] survival_3.5-5 utf8_1.2.3 TH.data_1.1-2
[58] leafem_0.2.0 e1071_1.7-13 scales_1.2.1
[61] sp_1.6-1 rmarkdown_2.22 zoo_1.8-12
[64] png_0.1-8 coda_0.19-4 evaluate_0.21
[67] knitr_1.43 tmaptools_3.1-1 s2_1.1.4
[70] rlang_1.1.1 Rcpp_1.0.10 glue_1.6.2
[73] DBI_1.1.3 jsonlite_1.8.5 R6_2.5.1
[76] units_0.8-2 Reload data from pervious session
load("_data/msoa2_spatial.RData")3.1 Spatial interdependence
We can not only use coordinates and geo-spatial information to connect different data sources, we can also explicitly model spatial (inter)dependence in the analysis of our data. In many instance, accounting for spatial dependence might even be necessary to avoid biased point estimates and standard errors, as observations are often not independent and identically distributed.
Tobler’s first law of geography has been used extensively (11,584 citation in 2023-06) to describe spatial dependence: ‘Everything is related to everything else, but near things are more related than distant things’ (Tobler 1970).
Note
Tobler’s first law is a bit of story
And it has been labeled as an excuse to not think too much about the reasons for spatial dependence or auto-correlation. For instance, measurement error, omitted variables, or inappropriate levels of aggregation are among reasons for auto-correlation (Pebesma and Bivand 2023).
We will come back to the reasons of spatial dependence. However, for now, we are interested in some tools to detect and analyse spatial relations.
To analyse spatial relations, we first need to define some sort of connectivity between units (e.g. similar to network analysis). There are some obvious candidates that be used to define these relations here: adjacency and proximity.
3.2 \(\boldsymbol{\mathbf{W}}\): Connectivity between units
The connectivity between units is usually represented in a matrix \(\boldsymbol{\mathbf{W}}\). There is an ongoing debate about the importance of spatial weights for spatial econometrics and about the right way to specify weights matrices (LeSage and Pace 2014; Neumayer and Plümper 2016). The following graph shows some possible options in how to define connectivity between units.
In spatial econometrics, the spatial connectivity (as shown above) is usually represented by a spatial weights matrix \({\boldsymbol{\mathbf{W}}}\): \[ \begin{equation} \boldsymbol{\mathbf{W}} = \begin{bmatrix} w_{11} & w_{12} & \dots & w_{1n} \\ w_{21} & w_{22} & \dots & w_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ w_{n1} & w_{n2} & \dots & w_{nn} \end{bmatrix} \end{equation} \] The spatial weights matrix \(\boldsymbol{\mathbf{W}}\) is an \(N \times N\) dimensional matrix with elements \(w_{ij}\) specifying the relation or connectivity between each pair of units \(i\) and \(j\).
Note: The diagonal elements \(w_{i,i}= w_{1,1}, w_{2,2}, \dots, w_{n,n}\) of \(\boldsymbol{\mathbf{W}}\) are always zero. No unit is a neighbour of itself. This is not true for spatial multiplier matrices (as we will see later).
3.2.1 Contiguity weights
A very common type of spatial weights. Binary specification, taking the value 1 for neighbouring units (queens: sharing a common edge; rook: sharing a common border), and 0 otherwise.
Contiguity weights \(w_{i,j}\), where
\[\begin{equation} w_{i,j} = \begin{cases} 1 & \text{if $i$ and $j$ neighbours}\\ 0 & \text{otherwise} \end{cases} \end{equation}\]A contiguity weights matrix with three units, where unit 1 and unit 3 are neighbours, while unit 2 has no neighbours would look like this:
\[ \begin{equation} \boldsymbol{\mathbf{W}} = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{bmatrix} \nonumber \end{equation} \]
Sparse matrices
Problem of `island’: units without neighbours (if I calculate an average of their neigbours, would that be zero, or NA, or a mean?)
Lets create a contiguity weights matrix (Queens neighbours) for the London MSOAs: we create a neighbours list (nb) using poly2nb(), which is an efficient way of storing \({\boldsymbol{\mathbf{W}}}\). A snap of 1 meter accounts for potential lacks of accuracy between lines and points.
# Contiguity (Queens) neighbours weights
queens.nb <- poly2nb(msoa.spdf,
queen = TRUE, # a single shared boundary point meets the contiguity condition
snap = 1) # we consider points in 1m distance as 'touching'
summary(queens.nb)Neighbour list object:
Number of regions: 983
Number of nonzero links: 5648
Percentage nonzero weights: 0.5845042
Average number of links: 5.745677
Link number distribution:
2 3 4 5 6 7 8 9 10 11 12 13
9 39 130 264 273 169 66 19 5 6 2 1
9 least connected regions:
160 270 475 490 597 729 755 778 861 with 2 links
1 most connected region:
946 with 13 links# Lets plot that
plot(st_geometry(msoa.spdf), border = "grey60")
plot(queens.nb, st_centroid(st_geometry(msoa.spdf)),
add = TRUE, pch = 19, cex = 0.6)
# We can also transform this into a matrix W
W <- nb2mat(queens.nb, style = "B")
print(W[1:10, 1:10]) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
1 0 0 0 0 0 0 0 0 0 0
2 0 0 1 0 0 0 0 0 0 0
3 0 1 0 0 1 0 0 0 0 0
4 0 0 0 0 0 1 0 0 0 1
5 0 0 1 0 0 1 1 0 0 0
6 0 0 0 1 1 0 1 0 1 1
7 0 0 0 0 1 1 0 1 1 0
8 0 0 0 0 0 0 1 0 0 0
9 0 0 0 0 0 1 1 0 0 1
10 0 0 0 1 0 1 0 0 1 0
Question
Among those first 10 units that you see above, which are the neighbours of unit number 6?
Why is the diagonal of this matrix all zero?
Overall, the matrix W has dimensions \(N \times N\), a row and a column for each observation. The value in a cell shows how units \(i\) (row number) and \(j\) (column number) are related to each other.
dim(W)[1] 983 983The row and column sums indicate the number of neighbours of each observation.
rowSums(W)[1:10] 1 2 3 4 5 6 7 8 9 10
11 6 7 5 5 6 6 6 6 5 colSums(W)[1:10] [1] 11 6 7 5 5 6 6 6 6 5Adjacency or graph-based neighbour’s weights matrices are usually symmetric. If unit 1 is a neighbour of unit 55, then unit 55 is also a neighbour of unit 1.
Higher Order Neighbours
Your neighbours have neighbours too, and they are called higher (second) order neighbours. The neighbours of your neighbour’s neighbours are third order neighbours.
You can use nblag() to calculate higher order neighbour relations.
3.2.2 Distance based weights
Another common type uses the distance \(d_{ij}\) between each unit \(i\) and \(j\).
- Inverse distance weights \(w_{i,j} = \frac{1}{d_{ij}^\alpha}\), where \(\alpha\) define the strength of the spatial decay.
\[ \begin{equation} \boldsymbol{\mathbf{W}} = \begin{bmatrix} 0 & \frac{1}{d_{ij}^\alpha} & \frac{1}{d_{ij}^\alpha} \\ \frac{1}{d_{ij}^\alpha} & 0 & \frac{1}{d_{ij}^\alpha} \\ \frac{1}{d_{ij}^\alpha} & \frac{1}{d_{ij}^\alpha} & 0 \end{bmatrix} \nonumber \end{equation} \]
Dense matrices
Specifying thresholds may be useful (to get rid of very small non-zero weights)
For now, we will just specify a neighbours list with a distance threshold of 3km using dnearneigh(). An alternative would be k nearest neighbours using knearneigh(). We will do the inverse weighting later.
# Crease centroids
coords <- st_geometry(st_centroid(msoa.spdf))Warning: st_centroid assumes attributes are constant over
geometries# Neighbours within 3km distance
dist_3.nb <- dnearneigh(coords, d1 = 0, d2 = 3000)
summary(dist_3.nb)Neighbour list object:
Number of regions: 983
Number of nonzero links: 22086
Percentage nonzero weights: 2.285652
Average number of links: 22.46796
Link number distribution:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
4 3 7 13 11 14 14 17 26 22 26 30 33 34 46 34 59 43 38 30 25 19
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
22 15 21 14 23 17 17 23 28 19 26 24 29 24 27 25 22 18 8 10 12 5
45 46 47
3 2 1
4 least connected regions:
158 160 463 959 with 1 link
1 most connected region:
545 with 47 links# Lets plot that
plot(st_geometry(msoa.spdf), border = "grey60")
plot(dist_3.nb, coords,
add = TRUE, pch = 19, cex = 0.6)
And you can see that the matrix is not so sparse anymore:
W2 <- nb2mat(dist_3.nb, style = "B")
W2[1:10, 1:10] [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
1 0 0 0 0 0 0 0 0 0 0
2 0 0 1 0 1 0 0 0 0 0
3 0 1 0 0 1 1 1 0 0 0
4 0 0 0 0 1 1 1 0 1 1
5 0 1 1 1 0 1 1 1 1 1
6 0 0 1 1 1 0 1 1 1 1
7 0 0 1 1 1 1 0 1 1 1
8 0 0 0 0 1 1 1 0 1 0
9 0 0 0 1 1 1 1 1 0 1
10 0 0 0 1 1 1 1 0 1 03.3 Normalization of \({\boldsymbol{\mathbf{W}}}\)
Normalizing ensures that the parameter space of the spatial multiplier is restricted to \(-1 < \rho > 1\), and the multiplier matrix is non-singular (don’t worry, more on this later).
The main message: Normalizing your weights matrix is always a good idea. Otherwise, the spatial parameters might blow up – if you can estimate the model at all. It also ensure easy interpretation of spillover effects.
Again, how to normalize a weights matrix is subject of debate (LeSage and Pace 2014; Neumayer and Plümper 2016).
3.3.1 Row-normalization
Row-normalization divides each non-zero weight by the sum of all weights of unit \(i\), which is the sum of the row.
\[ \frac{w_{ij}}{\sum_j^n w_{ij}} \]
With contiguity weights, spatially lagged variables contain mean of this variable among the neighbours of \(i\)
Proportions between units such as distances get lost (can be bad!)
Can induce asymmetries: \(w_{ij} \neq w_{ji}\)
For instance, we can use row-normalization for the Queens neighbours created above, and create a neighbours list with spatial weights.
queens.lw <- nb2listw(queens.nb,
style = "W") # W ist row-normalization
summary(queens.lw)Characteristics of weights list object:
Neighbour list object:
Number of regions: 983
Number of nonzero links: 5648
Percentage nonzero weights: 0.5845042
Average number of links: 5.745677
Link number distribution:
2 3 4 5 6 7 8 9 10 11 12 13
9 39 130 264 273 169 66 19 5 6 2 1
9 least connected regions:
160 270 475 490 597 729 755 778 861 with 2 links
1 most connected region:
946 with 13 links
Weights style: W
Weights constants summary:
n nn S0 S1 S2
W 983 966289 983 355.1333 4017.47To see what happened, let’s look at our example in matrix format again.
# transform into matrix with row-normalization
W_norm <- nb2mat(queens.nb, style = "W")
print(W_norm[1:10, 1:10]) [,1] [,2] [,3] [,4] [,5] [,6] [,7]
1 0 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
2 0 0.0000000 0.1666667 0.0000000 0.0000000 0.0000000 0.0000000
3 0 0.1428571 0.0000000 0.0000000 0.1428571 0.0000000 0.0000000
4 0 0.0000000 0.0000000 0.0000000 0.0000000 0.2000000 0.0000000
5 0 0.0000000 0.2000000 0.0000000 0.0000000 0.2000000 0.2000000
6 0 0.0000000 0.0000000 0.1666667 0.1666667 0.0000000 0.1666667
7 0 0.0000000 0.0000000 0.0000000 0.1666667 0.1666667 0.0000000
8 0 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1666667
9 0 0.0000000 0.0000000 0.0000000 0.0000000 0.1666667 0.1666667
10 0 0.0000000 0.0000000 0.2000000 0.0000000 0.2000000 0.0000000
[,8] [,9] [,10]
1 0.0000000 0.0000000 0.0000000
2 0.0000000 0.0000000 0.0000000
3 0.0000000 0.0000000 0.0000000
4 0.0000000 0.0000000 0.2000000
5 0.0000000 0.0000000 0.0000000
6 0.0000000 0.1666667 0.1666667
7 0.1666667 0.1666667 0.0000000
8 0.0000000 0.0000000 0.0000000
9 0.0000000 0.0000000 0.1666667
10 0.0000000 0.2000000 0.0000000
Question
Overall, how many neighbours does unit 9 have (including all columns)? How do you know?
rowSums(W)[9]We can also use the nb object to see which ones the neighbours are. Here, for instance, neighbours of unit 6:
queens.nb[6][[1]]
[1] 4 5 7 9 10 462This fits to what we see in the matrix above.
Warning
Note that row-normalization has some undesirable properties when we use some non-contigutiy based neighbour relations, such as distance based neighbours.
The problem: It obscures the proportion due to dividing by a row-specific value.
Let’s construct a hypothetical example
# Subset of 5 units
sub.spdf <- msoa.spdf[c(4, 5, 6, 102, 150), ]
mapview(sub.spdf)We construct the inverse-distance weighted 2 nearest neighbors.
# 2 closest neighbours
sub.coords <- st_geometry(st_centroid(sub.spdf))Warning: st_centroid assumes attributes are constant over
geometriesknn.nb <- knearneigh(sub.coords,
k = 2) # number of nearest neighboursWarning in knearneigh(sub.coords, k = 2): k greater than one-third
of the number of data pointsknn.nb <- knn2nb(knn.nb)
summary(knn.nb)Neighbour list object:
Number of regions: 5
Number of nonzero links: 10
Percentage nonzero weights: 40
Average number of links: 2
Non-symmetric neighbours list
Link number distribution:
2
5
5 least connected regions:
1 2 3 4 5 with 2 links
5 most connected regions:
1 2 3 4 5 with 2 links# listw with inverse-distance based weights
sub.lw <- nb2listwdist(knn.nb,
x = sub.coords, # needed for idw
type = "idw", # inverse distance weighting
alpha = 1, # the decay parameter for distance weighting
style = "raw") # without normalization
W_sub <- listw2mat(sub.lw)
formatC(W_sub, format = "f", digits = 6) [,1] [,2] [,3] [,4] [,5]
1 "0.000000" "0.000414" "0.000723" "0.000000" "0.000000"
2 "0.000414" "0.000000" "0.000962" "0.000000" "0.000000"
3 "0.000723" "0.000962" "0.000000" "0.000000" "0.000000"
4 "0.000000" "0.000033" "0.000032" "0.000000" "0.000000"
5 "0.000049" "0.000000" "0.000049" "0.000000" "0.000000"As you can see, units 1, 2, 3 have relatively proximate neighbours (.e.g inverse distance 0.000962: 3 zeros). Units 4 and 5, in contrast, have only very distant neighbours (e.g. inverse distance 0.000049: 4 zeros).
Now, see what happens when we use row-normalization.
sub.lw <- nb2listwdist(knn.nb,
x = sub.coords, # needed for idw
type = "idw", # inverse distance weighting
alpha = 1, # the decay parameter for distance weighting
style = "W") # for row normalization
W_sub <- listw2mat(sub.lw)
formatC(W_sub, format = "f", digits = 6) [,1] [,2] [,3] [,4] [,5]
1 "0.000000" "0.364083" "0.635917" "0.000000" "0.000000"
2 "0.300879" "0.000000" "0.699121" "0.000000" "0.000000"
3 "0.429123" "0.570877" "0.000000" "0.000000" "0.000000"
4 "0.000000" "0.507955" "0.492045" "0.000000" "0.000000"
5 "0.499360" "0.000000" "0.500640" "0.000000" "0.000000"All rows sum up to 1, but the strength of the relation is now similar for the distant units 4 and 5, and the proximate units 1, 2, 3.
3.3.2 Maximum eigenvalues normalization
Maximum eigenvalues normalization divides each non-zero weight by the overall maximum eigenvalue \(\lambda_{max}\). Each element of \(\boldsymbol{\mathbf{W}}\) is divided by the same scalar parameter, which preserves the relations.
\[ \frac{\boldsymbol{\mathbf{W}}}{\lambda_{max}} \]
Interpretation may become more complicated
Keeps proportions of connectivity strengths across units (relevant esp. for distance based \(\boldsymbol{\mathbf{W}}\))
We use eigenvalue normalization for the inverse distance neighbours. We use nb2listwdist() to create weight inverse distance based weights and normalize in one step.
coords <- st_geometry(st_centroid(msoa.spdf))Warning: st_centroid assumes attributes are constant over
geometriesidw.lw <- nb2listwdist(dist_3.nb,
x = coords, # needed for idw
type = "idw", # inverse distance weighting
alpha = 1, # the decay parameter for distance weighting
style = "minmax") # for eigenvalue normalization
summary(idw.lw)Characteristics of weights list object:
Neighbour list object:
Number of regions: 983
Number of nonzero links: 22086
Percentage nonzero weights: 2.285652
Average number of links: 22.46796
Link number distribution:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
4 3 7 13 11 14 14 17 26 22 26 30 33 34 46 34 59 43 38 30 25 19
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
22 15 21 14 23 17 17 23 28 19 26 24 29 24 27 25 22 18 8 10 12 5
45 46 47
3 2 1
4 least connected regions:
158 160 463 959 with 1 link
1 most connected region:
545 with 47 links
Weights style: minmax
Weights constants summary:
n nn S0 S1 S2
minmax 983 966289 463.6269 23.92505 1117.636Examples from above: See how this keeps the proportions in our example. Instead of transforming values to sum up to 1 in each row, we now have much smaller values for 4 and 5 then we have for the proximate units 1, 2, 3.
sub.lw <- nb2listwdist(knn.nb,
x = sub.coords, # needed for idw
type = "idw", # inverse distance weighting
alpha = 1, # the decay parameter for distance weighting
style = "minmax") # for eigenvalue normalization
W_sub <- listw2mat(sub.lw)
formatC(W_sub, format = "f", digits = 6) [,1] [,2] [,3] [,4] [,5]
1 "0.000000" "0.245687" "0.429123" "0.000000" "0.000000"
2 "0.245687" "0.000000" "0.570877" "0.000000" "0.000000"
3 "0.429123" "0.570877" "0.000000" "0.000000" "0.000000"
4 "0.000000" "0.019663" "0.019047" "0.000000" "0.000000"
5 "0.029099" "0.000000" "0.029174" "0.000000" "0.000000"3.4 Islands / missings
In practice, we often have a problem with islands. If we use contiguity based or distance based neighbour definitions, some units may end up with empty neighbours sets: they just do not touch any other unit and do not have a neighbour within a specific distance. This however creates a problem: what is the value in the neighbouring units?
The zero.policy option in spdep allows to proceed with empty neighbours sets. However, many further functions may run into problems and return errors. It often makes sense to either drop islands, to choose weights which always have neighbours (e.g. k nearest), or impute empty neighbours sets by using the nearest neighbours.