[2] Arbia , Giuseppe, (2014). A Primer for Spatial Econometrics : With Applications in R. Palgrave Macmillan UK. (available at UCT library online resource)
[3] Elhorst , Paul. J. (2014). Spatial Econometrics : From Cross-Sectional Data to Spatial Panels. Springer Berlin
[4] Gelfand, A.E., Diggle, P., Guttorp , P., Fuentes, M. and Fitzmaurice, G. (2010). Handbook of Spatial Statistics. Chapman & Hall/CRC. (available at UCT library online resource) Brilliant resource to have. Part 4 is dedicated to SPPs. *
[6] Haining , Robert, (2004). Spatial Data Analysis, Theory and Practice. Cambridge University Press. *
Throughout these lectures we use the Boston house-price data, Columbus crime data, and EU regional data collected at NUTS 2 level.
Boston data: classic Harrison and Rubinfeld (1978) housing data, with 506 census tracts and 23 variables.
Columbus data: spatial crime data for Columbus.
EU2 data: regional data for the second level of the Nomenclature of Territorial Units for Statistics (NUTS 2), widely used by Eurostat and other EU bodies.
A spatial data set includes the geometry of the spatial units, spatial coordinates or other spatial indexing information, and the attributes associated with each unit. These attributes may include economic, demographic, social, environmental, or other variables.
The three example data sets illustrate the same lattice-data idea: each observation belongs to a spatial unit, and each unit has both geometry and attributes. The geometry lets us construct neighbourhood relationships; the attributes are the variables we analyse.
Tip — click for more
A lattice data set is not just a data frame with coordinates. Its geometry carries the information needed to decide which areas are adjacent or close. In modern R workflows this geometry is usually stored in an sf object.
library(sf)library(spData)example(columbus)
colmbs> if (requireNamespace("sf", quietly = TRUE)) {
colmbs+ columbus <- sf::st_read(system.file("shapes/columbus.gpkg", package="spData")[1])
colmbs+ plot(sf::st_geometry(columbus))
colmbs+ }
Reading layer `columbus' from data source
`C:\Users\01438475\AppData\Local\R\win-library\4.4\spData\shapes\columbus.gpkg'
using driver `GPKG'
Simple feature collection with 49 features and 20 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 5.874907 ymin: 10.78863 xmax: 11.28742 ymax: 14.74245
Projected CRS: Undefined Cartesian SRS with unknown unit
A major pleasure in working with spatial data is their visualisation. Maps are amongst the most compelling graphics.
The Columbus data provide a spatially referenced example used throughout the lecture.
plot(sf::st_geometry(columbus))
library(tmap)tm_shape(columbus) +tm_polygons()
Plot of Crime Variable
tm_shape(columbus) +tm_polygons(style="quantile", col ="CRIME") +tm_legend(outside =TRUE, text.size = .8)
EU Nuts2 Dataset
EU NUTS 2 level GDP per capita (Euros), 2009.
Analysing Spatial Lattice Data
Tobler’s first law of geography:
“Everything is related to everything else, but close things are more related than things that are far apart.”
From this definition, we need to identify:
what we mean by spatial dependence;
how we measure spatial dependence;
who the neighbours of each spatial unit are; and
what the strength of the relationship is.
For each unit we therefore need a list of neighbours and a corresponding set of spatial weights.
Tip — click for more
Tobler’s law is a motivation, not a statistical model. Spatial statistics makes it operational by defining which locations are close, how strongly they are connected, and whether neighbouring values are more similar than would be expected under a suitable null model.
Definition of Spatial Dependence
Spatial dependence in a collection of sample observations refers to the fact that an observation associated with location \(i\) depends on observations at locations \(j\), where \(i\neq j\).
Why does spatial dependence occur? Two main reasons are highlighted by LeSage (1998, pp. 4–5):
measurement error related to the data-collection process for geographical units;
the spatial dimension of the research problem may itself be important in explaining the process under study.
Tip — click for more
Spatial dependence may be substantive—because processes spill across borders—or partly induced by how data are measured and aggregated. This is one reason the choice of spatial unit and spatial weights needs substantive justification.
Positive and Negative Spatial Dependence
Positive spatial dependence
Objects—points, areas, or lines—that are close in space tend to assume similar values. This produces spatial clustering of similar values.
Negative spatial dependence
Objects that are close in space tend to assume dissimilar values. This can produce an alternating or checkerboard-like arrangement.
Spatial dependence can therefore take the form of spatial correlation.
Ref: Giuseppe Arbia lecture notes, SEAI 2011.
Tip — click for more
Positive and negative dependence describe the sign of local association. They should not be confused with high and low values themselves: a low value surrounded by low values is still positive spatial association.
Exploring Spatial Dependence
The level of spatial dependence can be explored using graphical techniques and statistical measures:
graphical representations of spatial dependence;
Moran’s \(I\);
Geary’s \(C\).
Cliff and Ord (1972, 1975, 1981) found Moran’s \(I\) to be a particularly useful measure in many settings.
Tip — click for more
Moran’s I and Geary’s C both combine attribute information with spatial relationships, but they emphasise different aspects. Moran’s I works with cross-products of deviations from the mean; Geary’s C works with squared pairwise differences.
Global Moran’s I Statistic
Global Moran’s \(I\) combines attribute similarity and locational similarity.
For a variable \(y\), let
\[
z_i = y_i-\bar y,
\]
and let \(W^s=(w^s_{ij})\) denote the row standardizesd spatial weights matrix. The general form of Moran’s \(I\) is
\[
I =
\frac{Z'W^sZ}{Z'Z}
\]
where
\[
\sum_{i,j=1}^{n}w^s_{ij}=n.
\]
For a row-standardised weights matrix, the spatial lag \(W^sZ\) is the weighted mean of neighbouring standardised observations.
Tip — click for more
The spatial weights matrix is not a minor technical choice. It defines the spatial comparison being made. Changing \(W\) can change the value, interpretation, and significance of Moran’s \(I\).
How to Construct Weights’ Matrix
For geographic lattice data, there are two broad choices for \(W\) (Arbia, 2006, pp.37-38, Fischer et al., 2010, pp.66-68, LeSage, 1999, pp.11-13):
Boundary-based weights
Rook
Bishop
Queen
Distance-based weights
\(k\)-nearest-neighbour connectivity
critical cut-off connectivity
inverse-distance connectivity
distance-decay functions
combinations of these approaches
Tip — click for more
Boundary-based weights are often natural for administrative polygons. Distance-based weights can be preferable when interaction does not stop at administrative borders, or when the units are points/centroids rather than contiguous polygons.
Defining Spatial Neighbours: Boundary Based
Neighbourhood relationships can be defined in several ways (Arbia, 2005, p. 37).
For polygon data, a simple definition is contiguity: polygons \(i\) and \(j\) are neighbours when they share part of their boundary.
The three forms considered here are:
Rook: common side;
Bishop: common vertex only;
Queen: common side or common vertex.
Tip — click for more
Robustness matters. If conclusions change dramatically when moving from rook to queen contiguity, or from one distance threshold to another, that sensitivity is itself an important empirical finding.
Spatial Weights Based on Boundaries
A hypothetical five-region configuration is used to construct and interpret a binary contiguity matrix \(W\).
Spatial Weights Based on Boundaries: Rook Contiguity
Rook contiguity: define
\[
w_{ij} =
\begin{cases}
1, & \text{if regions } i \text{ and } j \text{ share a common side},\\
0, & \text{otherwise}.
\end{cases}
\]
The diagonal is normally set to zero, \(w_{ii}=0\), because a region is not treated as its own neighbour.
Tip — click for more
Rook contiguity is stricter than queen contiguity because touching at a corner is not enough. For grid cells this corresponds to the four orthogonal directions.
library(sf)library(spdep)# For a polygon sf object called x:rook_nb <-poly2nb(columbus, queen =FALSE)rook_w <-nb2listw(rook_nb, style ="W")head(rook_nb)
where \(d_{ij}\) is the adopted distance measure and \(d^*\) is a chosen critical cut-off.
Tip — click for more
For \(k\)-nearest neighbours every observation can be given exactly \(k\) outgoing neighbour links, which avoids isolates. A fixed distance threshold instead preserves a common physical scale but may create very different neighbour counts across dense and sparse areas.
Spatial Weights Based on Distances: 1st Nearest Neighbour – Obs (1)
For observation (1), calculate pairwise centroid distances and identify the closest observation. The first-nearest-neighbour link is then entered in the neighbour structure.
Spatial Weights Based on Distances: 1st Nearest Neighbour – Obs (2)
For observation (2), calculate pairwise centroid distances and identify the closest observation. The first-nearest-neighbour link is then entered in the neighbour structure.
Spatial Weights Based on Distances: 1st Nearest Neighbour – Obs (3)
For observation (3), calculate pairwise centroid distances and identify the closest observation. The first-nearest-neighbour link is then entered in the neighbour structure.
Spatial Weights Based on Distances: 1st Nearest Neighbour – Obs (4)
For observation (4), calculate pairwise centroid distances and identify the closest observation. The first-nearest-neighbour link is then entered in the neighbour structure.
Spatial Weights Based on Distances: 1st Nearest Neighbour – Obs (5)
For observation (5), calculate pairwise centroid distances and identify the closest observation. The first-nearest-neighbour link is then entered in the neighbour structure.
Spatial Weights Based on Distances: 1st Nearest Neighbour
Combining the first-nearest-neighbour relationships for all five observations gives the complete first-nearest-neighbour graph.
Neighbour list object:
Number of regions: 49
Number of nonzero links: 252
Percentage nonzero weights: 10.49563
Average number of links: 5.142857
2 disjoint connected subgraphs
Link number distribution:
1 2 3 4 5 6 7 8 9 10 11
4 8 6 2 5 8 6 2 6 1 1
4 least connected regions:
6 10 21 47 with 1 link
1 most connected region:
28 with 11 links
If you examine this neighbourhood list object a bit more closely you will see that, Polygon1 is neighbours with Poly2 and Poly3; Polygon2 is neighbours with Poly1 and Poly4; Polygon3 is neighbours with Poly1, Poly4 and Poly5 and so on:
dnbTresh1[1:3]
[[1]]
[1] 2 3
[[2]]
[1] 1 4
[[3]]
[1] 1 4 5
In calculation of the weights that will be used to create the spatially lagged variables, only the neighbours values are going to be used. For example for the 1st Polygon, the values of neighbouring Poly2 and Poly3 will be used. Before doing this the nb values need to be standardized, so that the row
Neighbours to Weights Matrix
Here we will work with the Critical Cut-off nearest neighbours
Choosing a threshold of \(68\) therefore guarantees that every centroid has at least one neighbour. As you can see row sum of the \(W^s\) is always 1, we row-standardized the W matrix.
Tip — click for more
Choosing the maximum first-nearest-neighbour distance as the threshold is a practical way to prevent isolates. It does not mean that this threshold is substantively optimal; it is a connectivity-based choice that should still be checked against the research problem.
Tip — click for more
Row standardisation makes every row sum to one, so each area’s neighbours collectively receive the same total weight. A consequence is that a region with two neighbours gives each neighbour more weight than a region with ten neighbours, all else equal.
The level of spatial dependence (spatial autocorrelation) can be explored with:
Moran’s \(I\);
Geary’s \(C\);
graphical representations such as the Moran scatter plot, correlogram, and variogram.
Tip — click for more
A spatial autocorrelation test asks whether the observed arrangement of values is more spatially structured than expected under a null model. Analytical tests and permutation tests use different routes to construct that reference distribution.
Some notation
To simplify notation:
\(n\): number of spatial areas;
\(i,j\): two areal units;
\(y_i\): value of the variable of interest at area \(i\);
\(z_i=y_i-\bar y\): mean-centred value at area \(i\);
\(w_{ij}\): measure of spatial connection between \(i\) and \(j\);
Global Measures of Spatial Autocorrelation: Moran’s I Statistic
Moran’s \(I\) is a global measure: it provides one summary of spatial autocorrelation for the entire spatial pattern.
\[
I =
\frac{Z'W^sZ}{Z'Z}
\]
With row-standardised \(W^s\), this is closely related to the cross-product between the centred variable \(Z'=Y-\bar{Y}\) and its spatial lag \(W^sZ\).
Tip — click for more
The numerator of Moran’s \(I\) is large and positive when connected areas tend to have deviations from the mean with the same sign. It becomes negative when neighbours tend to fall on opposite sides of the mean.
Moran’s I Expected Value and Variance
Under the randomisation null hypothesis, the expected value of Moran’s \(I\) is
\[
E(I)=-\frac{1}{n-1}.
\]
Thus, under the null, Moran’s \(I\) is not centred exactly at zero for a finite sample. As \(n\) increases,
\[
E(I)\longrightarrow 0.
\]
A standardised test statistic can be constructed as
Because \(E(I)=-1/(n-1)\), comparing an observed \(I\) directly with zero is only an approximation. Formal inference compares the observed statistic with its expected value and variability under the selected null model.
Comparing Moran’s I
Two broad approaches are used to assess whether an observed Moran’s \(I\) departs significantly from the null hypothesis of no spatial autocorrelation:
A Moran’s \(I\) value above its expected value indicates positive spatial autocorrelation.
A Moran’s \(I\) value below its expected value indicates negative spatial autocorrelation.
Moran’s \(I\) is a global measure: it does not by itself identify where local clusters or outliers occur.
Tip — click for more
The magnitude of a raw Moran’s \(I\) should not be compared casually across different data sets or different weights matrices. For inference, use the reference distribution appropriate to the selected \(W\) and null hypothesis.
Calculation of Moran’s I for the Columbus Dataset with Cut-off neighbourhood
For the Columbus CRIME variable, first centre the observations:
\[
z_i = y_i-\bar y.
\]
The first part of the centred vector in the original analysis is:
with a standard deviate of \(5.6206\) and a very small \(p\)-value.
moran.test(columbus$CRIME,dnbTresh1.listw)
Moran I test under randomisation
data: columbus$CRIME
weights: dnbTresh1.listw
Moran I statistic standard deviate = 5.6206, p-value = 9.514e-09
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.55182569 -0.02083333 0.01038068
Tip — click for more
The compact matrix expression is useful because it shows Moran’s \(I\) as a normalised spatial cross-product. It is also the bridge to spatial regression models, where \(Wy\) becomes the spatially lagged response.
Moran Scatter Plot: Y_s vs lagY_s
The Moran scatter plot places the centred or standardised variable on the horizontal axis and its spatial lag on the vertical axis:
\[
Z_i \quad \text{versus} \quad (WZ)_i.
\]
When \(W\) and \(z\) are standardised consistently, the slope of the fitted line is Moran’s \(I\).
The plot therefore links the global statistic to the local behaviour of individual observations.
lm_eqn <-function(df, y, x){ m <-lm(y ~ x, df); eq <-substitute(italic(y) == a + b %.%italic(x)*","~~italic(r)^2~"="~r2, list(a =format(unname(coef(m)[1]), digits =2),b =format(unname(coef(m)[2]), digits =2),r2 =format(summary(m)$r.squared, digits =3)))as.character(as.expression(eq));}dataMoransI =data.frame(Y_s, lagY_s)library(ggplot2)ggplot(dataMoransI, aes(x = Y_s, y = lagY_s)) +geom_point() +geom_smooth(method ="lm", se=FALSE)+geom_text(x =-10, y =20, label =lm_eqn(dataMoransI,lagY_s,Y_s), parse =TRUE)+geom_vline(xintercept=0, linetype="dashed", color ="red")+geom_hline(yintercept=0.78, linetype="dashed", color ="red")
Tip — click for more
The Moran scatter plot is descriptive. Influential observations can strongly affect the fitted slope, so ordinary regression \(t\)-tests on that plotted line are not a substitute for a spatial randomisation test.
The Moran scatter plot is partitioned into four quadrants:
Quadrant
Observation
Neighbours
Interpretation
High–High
High
High
positive association / high-value cluster
Low–Low
Low
Low
positive association / low-value cluster
High–Low
High
Low
negative association / spatial outlier
Low–High
Low
High
negative association / spatial outlier
The High–High and Low–Low quadrants represent positive spatial autocorrelation.
The High–Low and Low–High quadrants represent negative spatial autocorrelation.
Tip — click for more
Quadrants describe the sign pattern of a location and its spatial lag. Statistical significance is a separate question. A High–High point is not automatically a significant hot spot until a local inferential procedure is applied.
Geary’s C
Geary’s \(C\) focuses on squared differences between neighbouring observations:
Interpretation is opposite in direction to Moran’s \(I\):
\(C<1\): positive spatial autocorrelation;
\(C\approx1\): spatial randomness;
\(C>1\): negative spatial autocorrelation.
Tip — click for more
Geary’s C reacts strongly to local neighbour-to-neighbour differences. Its null reference value is approximately 1, which is why its interpretation runs in the opposite numerical direction to Moran’s \(I\).
Geary C test under randomisation
data: columbus$CRIME
weights: dnbTresh1.listw
Geary C statistic standard deviate = 5.2484, p-value = 1
alternative hypothesis: Expectation less than statistic
sample estimates:
Geary C statistic Expectation Variance
0.44927847 1.00000000 0.01101077
Spatial Correlogram
A non-parametric spatial correlogram is an alternative way of studying global spatial dependence.
Rather than beginning with one predefined weights matrix, it examines the association between pairs of observations as a function of their separation distance.
\(d_{ij}\) is the pairwise distance between locations \(i\) and \(j\);
\(f(\cdot)\) is a non-parametric function estimated from the data;
a LOWESS or kernel smoother may be used.
This allows us to ask how spatial association changes as distance increases.
Tip — click for more
A distance-based correlogram helps reveal the scale of spatial dependence. Positive association at short distances may decay, disappear, or even become negative as distance increases.
Local Measures and Tests for Spatial Autocorrelation
Local Indicators of Spatial Association (LISA) were introduced by Anselin (1995) to decompose global spatial autocorrelation into contributions from individual observations.
A common form of the local Moran statistic for area \(i\) is
\[
I_i
=
\frac{z_i}{m_2}
\sum_j w_{ij}z_j,
\]
where
\[
z_i=y_i-\bar y,
\qquad
m_2=\frac{1}{n}\sum_i z_i^2.
\]
The summation runs over the neighbourhood of area \(i\).
The collection of local Moran statistics is related to the global Moran statistic: their sum is proportional to global Moran’s \(I\).
Tip — click for more
Global Moran’s \(I\) can be significant while hiding substantial local heterogeneity. LISA statistics help reveal where the global pattern is being generated and where spatial outliers occur.
Identification of Local Spatial Clusters
Local spatial clusters—often called hot spots or cold spots—may be identified as locations, or sets of contiguous locations, for which a LISA statistic is significant.
The null hypothesis is no local spatial association.
Because the exact distribution of a generic LISA statistic can be difficult to obtain, a common alternative is a conditional randomisation / permutation procedure that produces empirical pseudo-\(p\)-values.
Tip — click for more
In conditional permutation for local Moran’s \(I_i\), the focal value at location \(i\) is held fixed while values at other locations are permuted. The observed \(I_i\) is then compared with this empirical reference distribution.
# Conditional permutation inference for local Moran statisticslocal_perm <-localmoran_perm( columbus$CRIME, dnbTresh1.listw,nsim =999)head(local_perm)
The output provides a local statistic for every spatial unit, together with quantities used for inference.
Tip — click for more
Modern spdep::localmoran() returns several columns rather than only \(I_i\). When teaching, it is useful to separate the local statistic itself, its expectation/variance under the chosen framework, and the resulting inferential quantities.
