Required current contributed CRAN packages:

I am running R 4.1.0, with recent update.packages().

needed <- c("rgeoda", "digest", "ggplot2", "tmap", "spdep", "spData", "sp", "sf")

Script

Script and data at https://github.com/rsbivand/UAM21_II/raw/main/UAM21_II_210608.zip. Download to suitable location, unzip and use as basis.

Seminar introduction

Schedule

Time Topic
Monday 7/6
09.00-12.00 How may we integrate spatial data from different sources? How may we aggregate spatial data? Which data structures are helpful in handling spatial data?
13.00-16.00 Observations of spatial data are related to each other either by distance, or by a graph of edges linking observations seen as being neighbours. How may they be constructed? How may we address the issues raised by the probable presence of spatial autocorrelation in the spatial data that we are using?
Tuesday 8/6
09.00-12.00 How can we measure global spatial autocorrelation?
13.00-16.00 How can we measure local spatial autocorrelation?
Monday 14/6
09.00-12.00 How may we fit regression models to spatial data in the presence of spatial autocorrelation? Maximum likelihood and spatial filtering, case weights.
13.00-16.00 How should we interpret the coefficients or impacts of spatial regression models? How may we predict from spatial regession models?
Tuesday 15/6
09.00-12.00 Multi-level models: at which level in the data may we fit spatial processes?
13.00-16.00 Spatial filtering, hierarchical GLM, GAM, etc., spatial epidemiological approaches
Monday 21/6
09.00-12.00 Presentations/consultations/discussion
13.00-16.00 Presentations/consultations/discussion 
library(sf)
## Linking to GEOS 3.10.0dev, GDAL 3.3.0, PROJ 8.0.1
data(pol_pres15, package="spDataLarge")
pol_pres15 |> 
    subset(select=c(TERYT, name, types)) |> 
    head()
## Simple feature collection with 6 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 235157.1 ymin: 366913.3 xmax: 281431.7 ymax: 413016.6
## Projected CRS: ETRS89 / Poland CS92
##    TERYT                name       types                       geometry
## 1 020101         BOLESŁAWIEC       Urban MULTIPOLYGON (((261089.5 38...
## 2 020102         BOLESŁAWIEC       Rural MULTIPOLYGON (((254150 3837...
## 3 020103            GROMADKA       Rural MULTIPOLYGON (((275346 3846...
## 4 020104        NOWOGRODZIEC Urban/rural MULTIPOLYGON (((251769.8 37...
## 5 020105          OSIECZNICA       Rural MULTIPOLYGON (((263423.9 40...
## 6 020106 WARTA BOLESŁAWIECKA       Rural MULTIPOLYGON (((267030.7 38...
library(spdep)
## Loading required package: sp
## Loading required package: spData
pol_pres15 |> poly2nb(queen=TRUE) -> nb_q
nb_q |> nb2listw(style="B") -> lw_q_B
nb_q |> nb2listw(style="W") -> lw_q_W
pol_pres15 |> st_geometry() |> st_centroid(of_largest_polygon=TRUE) -> coords 
coords |> dnearneigh(0, 18300) -> nb_d183
nb_d183 |> nbdists(coords) |> lapply(\(x) 1/(x/1000)) -> gwts
nb_d183 |> nb2listw(glist=gwts, style="B") -> lw_d183_idw_B

Measures of spatial autocorrelation

When analysing areal data, it has long been recognised that, if present, spatial autocorrelation changes how we may infer, relative to the default position of independent observations. In the presence of spatial autocorrelation, we can predict the values of observation \(i\) from the values observed at \(j \in N_i\), the set of its proximate neighbours. Early results (Moran 1948; Geary 1954), entered into research practice gradually, for example the social sciences (Duncan, Cuzzort, and Duncan 1961). These results were then collated and extended to yield a set of basic tools of analysis (Cliff and Ord 1973, 1981).

Cliff and Ord (1973) generalised and extended the expression of the spatial weights matrix representation as part of the framework for establishing the distribution theory for join count, Moran’s \(I\) and Geary’s \(C\) statistics. This development of what have become known as global measures, returning a single value of autocorrelation for the total study area, has been supplemented by local measures returning values for each areal unit (Getis and Ord 1992; L. Anselin 1995).

Measures and process mis-specification

Measures of spatial autocorrelation unfortunately pick up other mis-specifications in the way that we model data (Schabenberger and Gotway 2005; Daniel P. McMillen 2003). For reference, Moran’s \(I\) is given as (Cliff and Ord 1981, 17):

\[ I = \frac{n \sum_{(2)} w_{ij} z_i z_j}{S_0 \sum_{i=1}^{n} z_i^2} \] where \(x_i, i=1, \ldots, n\) are \(n\) observations on the numeric variable of interest, \(z_i = x_i - \bar{x}\), \(\bar{x} = \sum_{i=1}^{n} x_i / n\), \(\sum_{(2)} = \stackrel{\sum_{i=1}^{n} \sum_{j=1}^{n}}{i \neq j}\), \(w_{ij}\) are the spatial weights, and \(S_0 = \sum_{(2)} w_{ij}\). First we test a random variable using the Moran test, here under the normality assumption (argument randomisation=FALSE, default TRUE). Inference is made on the statistic \(Z(I) = \frac{I - E(I)}{\sqrt{\mathrm{Var}(I)}}\), the z-value compared with the Normal distribution for \(E(I)\) and \(\mathrm{Var}(I)\) for the chosen assumptions; this x does not show spatial autocorrelation with these spatial weights:

set.seed(1)
glance_htest <- function(ht) c(ht$estimate, 
    "Std deviate"=unname(ht$statistic), 
    "p.value"=unname(ht$p.value))
(pol_pres15 |> 
        nrow() |> 
        rnorm() -> x) |> 
    moran.test(lw_q_B, randomisation=FALSE, alternative="two.sided") |> 
    glance_htest()
## Moran I statistic       Expectation          Variance       Std deviate 
##     -0.0047715202     -0.0004009623      0.0001400449     -0.3693203351 
##           p.value 
##      0.7118889701

The test however fails to detect a missing trend in the data as a missing variable problem, finding spatial autocorrelation instead:

beta <- 0.0015
coords |> 
    st_coordinates() |> 
    subset(select=1, drop=TRUE) |> 
    (\(x) x/1000)() -> t
(x + beta * t -> x_t) |> 
    moran.test(lw_q_B, randomisation=FALSE, alternative="two.sided") |> 
    glance_htest()
## Moran I statistic       Expectation          Variance       Std deviate 
##      0.0434026880     -0.0004009623      0.0001400449      3.7014905445 
##           p.value 
##      0.0002143366

If we test the residuals of a linear model including the trend, the apparent spatial autocorrelation disappears:

lm(x_t ~ t) |> 
    lm.morantest(lw_q_B, alternative="two.sided") |> 
    glance_htest()
## Observed Moran I      Expectation         Variance      Std deviate 
##    -0.0047772213    -0.0007891850     0.0001397877    -0.3373064280 
##          p.value 
##     0.7358859145

A comparison of implementations of measures of spatial autocorrelation shows that a wide range of measures is available in R in a number of packages, chiefly in the spdep package, and that differences from other implementations can be attributed to design decisions (Roger S. Bivand and Wong 2018). The spdep package also includes the only implementations of exact and Saddlepoint approximations to global and local Moran’s I for regression residuals (Tiefelsdorf 2002; R. S. Bivand, Müller, and Reder 2009).

Global measures

We will begin by examining join count statistics, where joincount.test() takes a factor vector of values fx= and a listw object, and returns a list of htest (hypothesis test) objects defined in the stats package, one htest object for each level of the fx= argument. The observed counts are of neighbours with the same factor levels, known as same-colour joins.

args(joincount.test)
## function (fx, listw, zero.policy = NULL, alternative = "greater", 
##     sampling = "nonfree", spChk = NULL, adjust.n = TRUE) 
## NULL

The function takes an alternative= argument for hypothesis testing, a sampling= argument showing the basis for the construction of the variance of the measure, where the default "nonfree" choice corresponds to analytical permutation; the spChk= argument is retained for backward compatibility. For reference, the counts of factor levels for the type of municipality or Warsaw borough are:

(pol_pres15 |> 
        st_drop_geometry() |> 
        subset(select=types, drop=TRUE) -> Types) |> 
    table()
## 
##          Rural          Urban    Urban/rural Warsaw Borough 
##           1563            303            611             18

Since there are four levels, we re-arrange the list of htest objects to give a matrix of estimated results. The observed same-colour join counts are tabulated with their expectations based on the counts of levels of the input factor, so that few joins would be expected between for example Warsaw boroughs, because there are very few of them. The variance calculation uses the underlying constants of the chosen listw object and the counts of levels of the input factor. The z-value is obtained in the usual way by dividing the difference between the observed and expected join counts by the square root of the variance.

The join count test was subsequently adapted for multi-colour join counts (Upton and Fingleton 1985). The implementation as joincount.mult() in spdep returns a table based on nonfree sampling, and does not report p-values.

Types |> joincount.multi(listw=lw_q_B)
##                                Joincount   Expected   Variance  z-value
## Rural:Rural                   3087.00000 2793.92018 1126.53420   8.7320
## Urban:Urban                    110.00000  104.71854   93.29937   0.5468
## Urban/rural:Urban/rural        656.00000  426.52553  331.75903  12.5986
## Warsaw Borough:Warsaw Borough   41.00000    0.35018    0.34743  68.9646
## Urban:Rural                    668.00000 1083.94086  708.20864 -15.6297
## Urban/rural:Rural             2359.00000 2185.76854 1267.13133   4.8665
## Urban/rural:Urban              171.00000  423.72864  352.18954 -13.4669
## Warsaw Borough:Rural            12.00000   64.39253   46.45991  -7.6865
## Warsaw Borough:Urban             9.00000   12.48300   11.75800  -1.0158
## Warsaw Borough:Urban/rural       8.00000   25.17200   22.35382  -3.6320
## Jtot                          3227.00000 3795.48557 1496.39842 -14.6959

So far, we have used binary weights, so the sum of join counts multiplied by the weight on that join remains integer. If we change to row standardised weights, where the weights are not unity in all cases, the counts, expectations and variances change, but there are few major changes in the z-values.

Using an inverse distance based listw object does, however, change the z-values markedly, because closer centroids are upweighted relatively strongly:

Types |> joincount.multi(listw=lw_d183_idw_B)
##                                Joincount   Expected   Variance  z-value
## Rural:Rural                   3.4648e+02 3.6123e+02 4.9314e+01  -2.1003
## Urban:Urban                   2.9045e+01 1.3539e+01 2.2281e+00  10.3877
## Urban/rural:Urban/rural       4.6498e+01 5.5145e+01 9.6134e+00  -2.7891
## Warsaw Borough:Warsaw Borough 1.6822e+01 4.5275e-02 6.6083e-03 206.3805
## Urban:Rural                   2.0206e+02 1.4014e+02 2.3645e+01  12.7338
## Urban/rural:Rural             2.2517e+02 2.8260e+02 3.5892e+01  -9.5852
## Urban/rural:Urban             3.6499e+01 5.4784e+01 8.8586e+00  -6.1434
## Warsaw Borough:Rural          5.6502e+00 8.3253e+00 1.7260e+00  -2.0362
## Warsaw Borough:Urban          9.1800e+00 1.6139e+00 2.5392e-01  15.0150
## Warsaw Borough:Urban/rural    3.2676e+00 3.2545e+00 5.5180e-01   0.0177
## Jtot                          4.8183e+02 4.9072e+02 4.1570e+01  -1.3781

The implementation of Moran’s \(I\) in spdep in the moran.test() function has similar arguments to those of joincount.test(), but sampling= is replaced by randomisation= to indicate the underlying analytical approach used for calculating the variance of the measure. It is also possible to use ranks rather than numerical values (Cliff and Ord 1981, 46). The drop.EI2= agrument may be used to reproduce results where the final component of the variance term is omitted as found in some legacy software implementations.

args(moran.test)
## function (x, listw, randomisation = TRUE, zero.policy = NULL, 
##     alternative = "greater", rank = FALSE, na.action = na.fail, 
##     spChk = NULL, adjust.n = TRUE, drop.EI2 = FALSE) 
## NULL

The default for the randomisation= argument is TRUE, but here we will simply show that the test under normality is the same as a test of least squares residuals with only the intercept used in the mean model. The spelling of randomisation is that of Cliff and Ord (1973).

(pol_pres15 |> 
        st_drop_geometry() |> 
        subset(select=I_turnout, drop=TRUE) -> z) |> 
    moran.test(listw=lw_q_B, randomisation=FALSE) |> 
    glance_htest()
## Moran I statistic       Expectation          Variance       Std deviate 
##      0.6914339743     -0.0004009623      0.0001400449     58.4613487704 
##           p.value 
##      0.0000000000

The lm.morantest() function also takes a resfun= argument to set the function used to extract the residuals used for testing, and clearly lets us model other salient features of the response variable (Cliff and Ord 1981, 203). To compare with the standard test, we are only using the intercept here, and as can be seen, the results are the same.

lm(I_turnout ~ 1, pol_pres15) |> 
    lm.morantest(listw=lw_q_B) |> 
    glance_htest()
## Observed Moran I      Expectation         Variance      Std deviate 
##     0.6914339743    -0.0004009623     0.0001400449    58.4613487704 
##          p.value 
##     0.0000000000

The only difference between tests under normality and randomisation is that an extra term is added if the kurtosis of the variable of interest indicates a flatter or more peaked distribution, where the measure used is the classical measure of kurtosis. Under the default randomisation assumption of analytical randomisation, the results are largely unchanged.

(z |> 
    moran.test(listw=lw_q_B) -> mtr) |> 
    glance_htest()
## Moran I statistic       Expectation          Variance       Std deviate 
##      0.6914339743     -0.0004009623      0.0001400522     58.4598351527 
##           p.value 
##      0.0000000000

Of course, from the very beginning, interest was shown in Monte Carlo testing, also known as a Hope-type test and as a permutation bootstrap. By default, moran.mc() returns a "htest" object, but may simply use boot::boot() internally and return a "boot" object when return_boot=TRUE. In addition the number of simulations of the variable of interest by permutation, that is shuffling the values across the observations at random, needs to be given as nsim=.

set.seed(1)
z |> 
    moran.mc(listw=lw_q_B, nsim=999, return_boot = TRUE) -> mmc

The bootstrap permutation retains the outcomes of each of the random permutations, reporting the observed value of the statistic, here Moran’s \(I\), the difference between this value and the mean of the simulations under randomisation (equivalent to \(E(I)\)), and the standard deviation of the simulations under randomisation.

If we compare the Monte Carlo and analytical variances of \(I\) under randomisation, we typically see few differences, arguably rendering Monte Carlo testing unnecessary.

c("Permutation bootstrap"=var(mmc$t), 
  "Analytical randomisation"=unname(mtr$estimate[3]))
##    Permutation bootstrap Analytical randomisation 
##             0.0001441596             0.0001400522

Geary’s global \(C\) is implemented in geary.test() largely following the same argument structure as moran.test(). The Getis-Ord \(G\) test includes extra arguments to accommodate differences between implementations, as Bivand and Wong (2018) found multiple divergences from the original definitions, often to omit no-neighbour observations generated when using distance band neighbours. It is given by (Getis and Ord 1992, 194). For \(G_*\), the \(\sum_{(2)}\) constraint is relaxed by including \(i\) as a neighbour of itself (thereby also removing the no-neighbour problem, because all observations have at least one neighbour).

Finally, the empirical Bayes Moran’s \(I\) takes account of the denominator in assessing spatial autocorrelation in rates data (Assunção and Reis 1999). Until now, we have considered the proportion of valid votes cast in relation to the numbers entitled to vote by spatial entity, but using EBImoran.mc() we can try to accommodate uncertainty in extreme rates in entities with small numbers entitled to vote. There is, however, little impact on the outcome in this case.

Global measures of spatial autocorrelation using spatial weights objects based on graphs of neighbours are, as we have seen, rather blunt tools, which for interpretation depend critically on a reasoned mean model of the variable in question. If the mean model is just the intercept, the global measures will respond to all kinds of mis-specification, not only spatial autocorrelation. A key source of mis-specification will typically also include the choice of entities for aggregation of data.

ACS CV data set

df_tracts <- st_read("df_tracts.gpkg")
## Reading layer `df_tracts' from data source 
##   `/home/rsb/und/uam21/units_2/df_tracts.gpkg' using driver `GPKG'
## Simple feature collection with 71353 features and 28 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.7332 ymin: 24.5447 xmax: -66.9499 ymax: 49.38436
## Geodetic CRS:  NAD83
(nb_subset <- readRDS("nb_subset.rds"))
## Neighbour list object:
## Number of regions: 71353 
## Number of nonzero links: 438146 
## Percentage nonzero weights: 0.008605862 
## Average number of links: 6.140541 
## 17 regions with no links:
## 12745 26628 29077 29120 32398 32775 43104 46344 46390 52303 58113 58261 70000 70055 70058 70457 70458
nc_nb_subset <- n.comp.nb(nb_subset)
nc_nb_subset$nc
## [1] 31
table(table(nc_nb_subset$comp.id))
## 
##     1     2     4    13    15    17   107 71156 
##    17     4     5     1     1     1     1     1
set.ZeroPolicyOption(TRUE)
## [1] FALSE
lw <- nb2listw(nb_subset, style="W")
system.time(mt_med_inc_cv <- moran.test(df_tracts$med_inc_cv, lw, alternative="two.sided"))
##    user  system elapsed 
##   2.590   0.000   2.596
mt_med_inc_cv
## 
##  Moran I test under randomisation
## 
## data:  df_tracts$med_inc_cv  
## weights: lw  n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = 85.539, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      1.887741e-01     -1.401836e-05      4.871047e-06
form <- log(med_inc_cv) ~ log1p(vacancy_rate) + log1p(old_rate) + log1p(black_rate) + log1p(hisp_rate) + log1p(group_pop) + log1p(dens)
lm_mod <- lm(form, data=df_tracts)
summary(lm_mod)
## 
## Call:
## lm(formula = form, data = df_tracts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.14822 -0.29445  0.01038  0.30068  2.99263 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -6.2284926  0.0355078 -175.41   <2e-16 ***
## log1p(vacancy_rate)  0.4039100  0.0248345   16.26   <2e-16 ***
## log1p(old_rate)      0.4069179  0.0040408  100.70   <2e-16 ***
## log1p(black_rate)    0.5478924  0.0110139   49.75   <2e-16 ***
## log1p(hisp_rate)     0.5854893  0.0121825   48.06   <2e-16 ***
## log1p(group_pop)     0.0216907  0.0007901   27.45   <2e-16 ***
## log1p(dens)          0.0354277  0.0018286   19.37   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4526 on 71346 degrees of freedom
## Multiple R-squared:  0.2296, Adjusted R-squared:  0.2295 
## F-statistic:  3543 on 6 and 71346 DF,  p-value: < 2.2e-16
lm.morantest(lm_mod, lw)
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = form, data = df_tracts)
## weights: lw
## 
## Moran I statistic standard deviate = 50.637, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     1.117328e-01    -6.347763e-05     4.874291e-06

Local measures

Building on insights from the weaknesses of global measures, local indicators of spatial association began to appear in the first half of the 1990s (L. Anselin 1995; Getis and Ord 1992, 1996). In addition, the Moran plot was introduced, plotting the values of the variable of interest against their spatially lagged values, typically using row-standardised weights to make the axes more directly comparable (L. Anselin 1996). The moran.plot() function also returns an influence measures object used to label observations exerting more than propotional influence on the slope of the line representing global Moran’s \(I\). We can see that there are many spatial entities exerting such influence. These pairs of observed and lagged observed values make up in aggregate the global measure, but can also be explored in detail. The quadrants of the Moran plot also show low-low pairs in the lower left quadrant, high-high in the upper right quadrant, and fewer low-high and high-low pairs in the upper left and lower right quadrants.

z |> moran.plot(listw=lw_q_W, labels=pol_pres15$TERYT, cex=1, pch=".", xlab="I round turnout", ylab="lagged turnout") -> infl_W

If we extract the hat value influence measure from the returned object, the figure suggests that some edge entities exert more than proportional influence (perhaps because of row standardisation), as do entities in or near larger urban areas.

pol_pres15$hat_value <- infl_W$hat
library(tmap)
tm_shape(pol_pres15) + tm_fill("hat_value")

Bivand and Wong (2018) discuss issues impacting the use of local indicators, such as local Moran’s \(I\) and local Getis-Ord \(G\). Some issues affect the calculation of the local indicators, others inference from their values. Because \(n\) statistics may be being calculated from the same number of observations, there are multiple comparison problems that need to be addressed. Although the apparent detection of hotspots from values of local indicators has been quite widely adopted, it remains fraught with difficulty because adjustment of the inferential basis to accommodate multiple comparisons is not often chosen, and as in the global case, mis-specification also remains a source of confusion. Further, interpreting local spatial autocorrelation in the presence of global spatial autocorrelation is challenging (J. K. Ord and Getis 2001; Tiefelsdorf 2002; R. S. Bivand, Müller, and Reder 2009). The mlvar= and adjust.x= arguments to localmoran() are discussed in Bivand and Wong (2018), and permit matching with other implementations. The p.adjust.method= argument uses an untested speculation that adjustment should only take into account the cardinality of the neighbour set of each observation when adjusting for multiple comparisons; using stats::p.adjust() is preferable.

Taking "two.sided" p-values because these local indicators when summed and divided by the sum of the spatial weights, and thus positive and negative local spatial autocorrelation may be present, we obtain:

z |> 
    localmoran(listw=lw_q_W, alternative="two.sided") -> locm
all.equal(sum(locm[,1])/Szero(lw_q_W), unname(moran.test(z, lw_q_W)$estimate[1]))
## [1] TRUE

Using stats::p.adjust() to adjust for multiple comparisons, we see that almost 29% of the local measures have p-values < 0.05 if no adjustment is applied, but only 12% using Bonferroni adjustment, with two other choices also shown:

pva <- \(pv) cbind("none"=pv, "bonferroni"=p.adjust(pv, "bonferroni"), "fdr"=p.adjust(pv, "fdr"), "BY"=p.adjust(pv, "BY"))
locm |> 
    subset(select="Pr(z != 0)", drop=TRUE) |> 
    pva() -> pvsp
f <- \(x) sum(x < 0.05)
apply(pvsp, 2, f)
##       none bonferroni        fdr         BY 
##        715        297        576        424

In the global measure case, bootstrap permutations could be used as an alternative to analytical methods for possible inference. In the local case, conditional permutation may be used, retaining the value at observation \(i\) and randomly sampling from the remaining \(n-1\) values to find randomised values at neighbours, and is provided as localmoran_perm(), which will use multiple nodes to sample in parallel if provided, and permits the setting of a seed for the random number generator across the compute nodes:

library(parallel)
set.coresOption(ifelse(detectCores() == 1, 1, detectCores()-1L))
## NULL
system.time(z |> 
        localmoran_perm(listw=lw_q_W, nsim=499, alternative="two.sided", iseed=1) -> locm_p)
##    user  system elapsed 
##   0.675   0.300   0.251

The outcome is that almost 32% of observations have two sided p-values < 0.05 without multiple comparison adjustment, and under 3% with Bonferroni adjustment.

locm_p |> 
    subset(select="Pr(z != 0)", drop=TRUE) |> 
    pva() -> pvsp
apply(pvsp, 2, f)
##       none bonferroni        fdr         BY 
##        783         68        466        161

We can see what is happening by tabulating counts of the standard deviate of local Moran’s \(I\), where the two-sided \(\alpha=0.05\) bounds would be \(0.025\) and \(0.975\), but Bonferroni adjustment is close to \(0.00001\) and \(0.99999\). Without adjustment, almost 800 observations are significant, with Bonferroni adjustment, only 68 in the conditional permutation case:

brks <- qnorm(c(0, 0.00001, 0.0001, 0.001, 0.01, 0.025, 0.5, 0.975, 0.99, 0.999, 0.9999, 0.99999, 1))
(locm_p |> 
        subset(select=Z.Ii, drop=TRUE) |> 
        cut(brks) |> 
        table()-> tab)
## 
##  (-Inf,-4.26] (-4.26,-3.72] (-3.72,-3.09] (-3.09,-2.33] (-2.33,-1.96] 
##             0             0             1             4             3 
##     (-1.96,0]      (0,1.96]   (1.96,2.33]   (2.33,3.09]   (3.09,3.72] 
##           471          1241           182           320           142 
##   (3.72,4.26]   (4.26, Inf] 
##            63            68
sum(tab[c(1:5, 8:12)])
## [1] 783
sum(tab[c(1, 12)])
## [1] 68

Recent development: (Sauer et al. 2021)

z |> 
    localmoran(listw=lw_q_W, conditional=TRUE, alternative="two.sided") -> locm_c
locm_c |> 
    subset(select="Pr(z != 0)", drop=TRUE) |> 
    pva() -> pvsp
apply(pvsp, 2, f)
##       none bonferroni        fdr         BY 
##        789         69        468        156
pol_pres15$locm_Z <- locm[, "Z.Ii"]
pol_pres15$locm_c_Z <- locm_c[, "Z.Ii"]
pol_pres15$locm_p_Z <- locm_p[, "Z.Ii"]
tm_shape(pol_pres15) + tm_fill(c("locm_Z", "locm_c_Z", "locm_p_Z"), breaks=brks, midpoint=0, title="Standard deviates of\nLocal Moran's I") + tm_facets(free.scales=FALSE, ncol=2) + tm_layout(panel.labels=c("Analytical total", "Analytical conditional", "Conditional permutation"))

The figure shows that conditional permutation scales back the proportion of standard deviate values taking extreme values, especially positive values. As we will see below, the analytical standard deviates of local Moran’s \(I\) should probably not be used if alternatives are available.

In presenting local Moran’s \(I\), use is often made of “hotspot” maps. Because \(I_i\) takes high values both for strong positive autocorrelation of low and high values of the input variable, it is hard to show where “clusters” of similar neighbours with low or high values of the input variable occur. The quadrants of the Moran plot are used, by creating a categorical quadrant variable interacting the input variable and its spatial lag split at their means. The quadrant categories are then set to NA if, for the chosen standard deviate and adjustment, \(I_i\) would be considered insignificant. Here, for the conditional permutation standard deviates, Bonferroni adjusted, 14 observations belong to “Low-Low clusters,” and 54 to “High-High clusters”:

quadr <- interaction(cut(infl_W$x, c(-Inf, mean(infl_W$x), Inf), labels=c("Low X", "High X")), cut(infl_W$wx, c(-Inf, mean(infl_W$wx), Inf), labels=c("Low WX", "High WX")), sep=" : ")
a <- table(quadr)
pol_pres15$hs_an_q <- quadr
is.na(pol_pres15$hs_an_q) <- !(pol_pres15$locm_Z < brks[6] | pol_pres15$locm_Z > brks[8])
b <- table(pol_pres15$hs_an_q)
pol_pres15$hs_cp_q <- quadr
is.na(pol_pres15$hs_cp_q) <- !(pol_pres15$locm_p_Z < brks[2] | pol_pres15$locm_p_Z > brks[12])
c <- table(pol_pres15$hs_cp_q)
pol_pres15$hs_ac_q <- quadr
is.na(pol_pres15$hs_ac_q) <- !(pol_pres15$locm_c_Z < brks[2] | pol_pres15$locm_c_Z > brks[12])
d <- table(pol_pres15$hs_ac_q)
t(rbind("Moran plot quadrants"=a, "Unadjusted analytical total"=b, "Bonferroni analytical cond."=d, "Bonferroni cond. perm."=c))
##                  Moran plot quadrants Unadjusted analytical total
## Low X : Low WX                   1040                         370
## High X : Low WX                   264                           3
## Low X : High WX                   213                           0
## High X : High WX                  978                         342
##                  Bonferroni analytical cond. Bonferroni cond. perm.
## Low X : Low WX                            14                     14
## High X : Low WX                            0                      0
## Low X : High WX                            0                      0
## High X : High WX                          55                     54
tm_shape(pol_pres15) + tm_fill(c("hs_an_q", "hs_ac_q", "hs_cp_q"), colorNA="grey95", textNA="Not significant", title="Turnout hotspot status\nLocal Moran's I") + tm_facets(free.scales=FALSE, ncol=2) + tm_layout(panel.labels=c("Unadjusted analytical total", "Bonferroni analytical cond.", "Cond. perm. with Bonferroni"))

The figure shows the impact of using analytical or conditional permutation standard deviates, and no or Bonferroni adjustment, reducing the counts of observations in “Low-Low clusters” from 370 to 14, and “High-High clusters” from 342 to 54; the “High-High clusters” are metropolitan areas.

The local Getis-Ord \(G\) measure is reported as a standard deviate, and may also take the \(G^*\) form where self-neighbours are inserted into the neighbour object using include.self(). The observed and expected values of local \(G\) with their analytical variances may also be returned if return_internals=TRUE.

system.time(z |> 
        localG(lw_q_W) -> locG)
##    user  system elapsed 
##   0.004   0.000   0.004
system.time(z |> 
        localG_perm(lw_q_W, nsim=499, iseed=1) -> locG_p)
##    user  system elapsed 
##   0.246   0.149   0.200

Once again we face the problem of multiple comparisons, with the count of areal unit p-values < 0.05 being reduced by an order of magnitude when employing Bonferroni correction:

locG |> 
    c() |> 
    abs() |> 
    pnorm(lower.tail = FALSE) |> 
    (\(x) x*2)() |> 
    pva() -> pvsp
apply(pvsp, 2, f)
##       none bonferroni        fdr         BY 
##        789         69        468        156
library(ggplot2)
p1 <- ggplot(data.frame(Zi=locm_c[,4], Zi_perm=locm_p[,4])) + geom_point(aes(x=Zi, y=Zi_perm), alpha=0.2) + xlab("Analytical conditional") + ylab("Permutation conditional") + coord_fixed() + ggtitle("Local Moran's I")
p2 <- ggplot(data.frame(Zi=c(locG), Zi_perm=c(locG_p))) + geom_point(aes(x=Zi, y=Zi_perm), alpha=0.2) + xlab("Analytical conditional") + ylab("Permutation conditional") + coord_fixed() + ggtitle("Local G")
gridExtra::grid.arrange(p1, p2, nrow=1)

The figure shows that, keeping fixed aspect in both panels, conditional permutation changes the range and distribution of the standard deviate values for \(I_i\), but that for \(G_i\), the two sets of standard deviates are equivalent.

pol_pres15$locG_Z <- c(locG)
pol_pres15$hs_G <- cut(c(locG), c(-Inf, brks[2], brks[12], Inf), labels=c("Low", "Not significant", "High"))
table(pol_pres15$hs_G)
## 
##             Low Not significant            High 
##              14            2426              55
m1 <- tm_shape(pol_pres15) + tm_fill(c("locG_Z"), midpoint=0, title="Standard\ndeviate")
m2 <- tm_shape(pol_pres15) + tm_fill(c("hs_G"), title="Bonferroni\nhotspot status")
tmap_arrange(m1, m2, nrow=1)

As can be seen, we do not need to contrast the two estimation methods, and showing the mapped standard deviate is as informative as the “hotspot” status for the chosen adjustment. In the case of \(G_i\), the values taken by the measure reflect the values of the input variable, so a “High cluster” is found for observations with high values of the input variable, here high turnout in metropolitan areas.

ACS data set

system.time(localI_med_inc_cv <- localmoran(df_tracts$med_inc_cv, lw, alternative="greater"))
##    user  system elapsed 
##   1.442   0.002   1.455
system.time(localI_med_inc_cv_cond <- localmoran(df_tracts$med_inc_cv, lw, conditional=TRUE, alternative="greater"))
##    user  system elapsed 
##   1.542   0.001   1.568
library(parallel)
set.coresOption(detectCores()-1L)
# [1] 5
get.mcOption()
# [1] TRUE
system.time(localI_med_inc_cv_perm <- localmoran_perm(df_tracts$med_inc_cv, lw, nsim=499, alternative="greater", iseed=1))
#    user  system elapsed
# 113.378   4.688  26.053
#saveRDS(localI_med_inc_cv_perm, file="localI_med_inc_cv_perm.rds")
ACS_infl <- moran.plot(df_tracts$med_inc_cv, lw, labels=df_tracts$GEOID, cex=1, pch=".", xlab="CV", ylab="lagged CV")

quadr <- interaction(cut(ACS_infl$x, c(-Inf, mean(ACS_infl$x), Inf), labels=c("Low X", "High X")), cut(ACS_infl$wx, c(-Inf, mean(ACS_infl$wx), Inf), labels=c("Low WX", "High WX")), sep=" : ")
(a <- table(quadr))
## quadr
##   Low X : Low WX  High X : Low WX  Low X : High WX High X : High WX 
##            30022            11794            14117            15420
df_tracts$hs_an_q <- quadr
is.na(df_tracts$hs_an_q) <- !(localI_med_inc_cv[, 5] < 0.05)
b <- table(df_tracts$hs_an_q)
df_tracts$hs_ac_q <- quadr
is.na(df_tracts$hs_ac_q) <- !(localI_med_inc_cv_cond[, 5] < 0.05)
c <- table(df_tracts$hs_ac_q)
df_tracts$hs_acb_q <- quadr
is.na(df_tracts$hs_acb_q) <- !(localI_med_inc_cv_cond[, 5] < 0.0000227795)
d <- table(df_tracts$hs_acb_q)
t(rbind("Moran plot quadrants"=a, "Unadjusted analytical total"=b, "Unadjusted analytical cond."=c, "Bonferroni analytical cond."=d))
##                  Moran plot quadrants Unadjusted analytical total
## Low X : Low WX                  30022                        1782
## High X : Low WX                 11794                           5
## Low X : High WX                 14117                           0
## High X : High WX                15420                        3970
##                  Unadjusted analytical cond. Bonferroni analytical cond.
## Low X : Low WX                          3759                          12
## High X : Low WX                            5                           5
## Low X : High WX                            0                           0
## High X : High WX                        4589                         698
chicago_MA <- read.table("Chicago_MA.txt", colClasses=c("character", "character"))
chicago_MA_tracts <- !is.na(match(substring(df_tracts$GEOID, 1, 5), chicago_MA$V2))
tm_shape(df_tracts[chicago_MA_tracts,]) + tm_fill(c("hs_an_q", "hs_ac_q", "hs_acb_q"), colorNA="grey95", textNA="Not significant", title="CV hotspot status\nLocal Moran's I") + tm_facets(free.scales=FALSE, ncol=3) + tm_layout(panel.labels=c("Unadjusted total", "Unadjusted conditional", "Bonferroni conditional"))

New rgeoda

Very recently, Geoda has been wrapped for R as rgeoda (Li and Anselin 2021), and will provide very similar functionalities for the exploration of spatial autocorrelation in areal data as spdep. The active objects are kept as pointers to a compiled code workspace; using compiled code for all operations (as in Geoda itself) makes rgeoda perform fast, but leaves less flexible when modifications or enhancements are desired.

The contiguity neighbours it constructs are the same as those found by poly2nb(), as almost are the \(I_i\) measures. The difference is as established by Roger S. Bivand and Wong (2018), that localmoran() calculates the input variable variance divinding by \(n\), but Geoda uses \((n-1)\), as can be reproduced by setting mlvar=FALSE:

library(rgeoda)
## Loading required package: digest
## 
## Attaching package: 'rgeoda'
## The following object is masked from 'package:spdep':
## 
##     skater
system.time(Geoda_w <- queen_weights(pol_pres15))
## here
##    user  system elapsed 
##   0.066   0.000   0.067
summary(Geoda_w)
##                      name               value
## 1 number of observations:                2495
## 2          is symmetric:                 TRUE
## 3               sparsity: 0.00228786229774178
## 4        # min neighbors:                   1
## 5        # max neighbors:                  13
## 6       # mean neighbors:    5.70821643286573
## 7     # median neighbors:                   6
## 8           has isolates:               FALSE
system.time(lisa <- local_moran(Geoda_w, pol_pres15["I_turnout"], 
    cpu_threads=ifelse(parallel::detectCores() == 1, 1, parallel::detectCores()-1L), permutations=499, seed=1))
##    user  system elapsed 
##   0.150   0.005   0.064
all.equal(card(nb_q), lisa_num_nbrs(lisa), check.attributes=FALSE)
## [1] TRUE
all.equal(lisa_values(lisa), localmoran(pol_pres15$I_turnout, listw=lw_q_W, mlvar=FALSE)[,1], check.attributes=FALSE)
## [1] TRUE

Background to spatial econometrics

Antecedents

In the same way that (Fujita, Krugman, and Venables 1999) begin their study of the spatial economy by looking at the antecedents of their subject, it is helpful to place spatial econometrics in its temporal and academic context. This context is sufficiently different from the contemporary setting that it may be hard to grasp the background for many of the features of spatial econometrics that came into being during its earlier years. Indeed, the ranges of topics that were studied in economics in the 1960’s and 1970’s differ markedly from those in focus today. If we can sketch the context within which spatial econometrics was created, and its methods developed, we should be able to illuminate choices made then which influence our understanding and application of spatial econometric methods.

Critics of the practice of spatial econometrics, such as (Gibbons and Overman 2012), appear to overlook these antecedents, and consequently judge the potential of the field on a partial, perhaps anachronistic, understanding, viewing phenomena with a history in ahistorical way. Since we are attempting to provide an introduction to applied spatial econometrics, we need to throw light on the original motivations and concerns of the first scholars engaged in the field. (Luc Anselin 2010) indicates clearly and repeatedly [L. Anselin (1988); anselin:06; anselin:10a] that we should acknowledge Spatial Econometrics by (Paelinck and Klaassen 1979) of the Netherlands Economic Institute as our starting point, and so celebrates thirty years of spatial econometrics in 2009. This firm confirmation of the importance of Jean Paelinck’s contributions as scholar and community-builder is fully justified. We should then turn to the motivations given in (Paelinck and Klaassen 1979) to indicate which contextual factors were of importance at that time, and the breadth of the academic communities with which they were in contact.

In a recent short commentary, (Paelinck 2013) recalls his conviction, expressed in 1967, that “early econometric exercises \(\ldots\) relating only variables posessing the same regional index \(\ldots\) were inadequate to represent the correct spatial workings of the economy, which would then be reflected in the policy outcomes.” A year before, (Paelinck 2012) points to salient isomorphisms linking spatial regression models, simultaneous equation models and input-output models; these were known of and discussed in the early formative period of spatial econometrics. We will return in subsequent chapters to the ways in which spatial regression models may be specified, but for now, a simple presentation of these isomorphisms as perceived in the early period is sufficient:

\[ {\mathbf y} = {\mathbf A}{\mathbf y} + {\mathbf X}{\mathbf b} + {\mathbf \varepsilon} \]

is a spatial regression model where \({\mathbf A}\) is a matrix expressing the mutual first order spatial dependencies between regions — the similarity of this form and the Koyck distributed lag model is striking [Koyck (1954); Klein (1958); Griliches (1967)};

\[ {\mathbf A}{\mathbf y} + {\mathbf X}{\mathbf b} = {\mathbf \varepsilon} \]

is a simultaneous equation model where \({\mathbf A}\) is a matrix expressing the dependencies between the equations; and:

\[ {\mathbf y} = {\mathbf A}{\mathbf y} + {\mathbf f} \]

is an input-output model where \({\mathbf A}\) is a matrix of sectoral input-output coefficients, and \({\mathbf f}\) is final demand.

Input-output models, simultaneous equation models, and the importance of policy outcomes were all known intimately at the Netherlands Economic Institute at this time, and elsewhere among applied economists. The isomporhisms flowed from the known to the unknown, from the stuff of contemporary research and policy advice to doubts about the calibration of aspatial models, and on to what became termed spatial econometrics. If we compare these topics with those described for Regional Science by (Boyce 2004), we can see the outlines of research priorities at the time: including urban and regional models for planning, regional and interregional input-output models, transport and location models. During the 1960s and 1970s, many of these models were enhanced — matching needs for policy advice — to cover environmental questions, adding natural resources as inputs and pollution to outputs. Paelinck’s co-author in a key paper in spatial econometrics (Hordijk and Paelinck 1976), went on to work in environmental management and research.

Reading (Paelinck and Klaassen 1979), we see that the programme of research into the space economy undertaken at the Netherlands Economic Institute led first to the publication of (Paelinck and Nijkamp 1975), and then to (Klaassen, Paelinck, and Wagenaar 1979), published in the same year as Spatial Econometrics. All three books were published in the same series and appear to reflect the core concerns of economists at the Institute doing reasearch on regionalised national macro-economic models. The direct link to Jan Tinbergen is evident in the account of the context of economic research in the Netherlands given by (Theil 1964). If we take Paelinck at his word, he and his colleagues were aware that an aspatial regionalisation of national accounts, of input-output models, or transport models, might prejudice policy advice and outcomes through inadequate and inappropriate calibration.

(Klaassen, Paelinck, and Wagenaar 1979) is mainly concerned with model construction, while about a third of (Paelinck and Nijkamp 1975) is devoted to input-output analysis. Both books show sustained concern for economic measurement, especially of national accounts data, intersectoral transactions, and many other topics. Considerable attention is also paid to the data collection units, be they sectors or regions. The need to attempt to define regions that match the underlying economic realities was recognised clearly, and a key part of (Paelinck and Nijkamp 1975) is devoted to regionalisation, and the distinction between functional regions and homogeneous regional classifications is made. In the motivation for spatial econometric models given in (Paelinck and Klaassen 1979), consumption and investment in a region are modelled as depending on income both in the region itself and in its contiguous neighbours, termed a “spatial income-generating model.” It became important to be able to calibrate planning models of this kind to provide indications of the possible outcomes of alternative policy choices, hence the need for spatial econometrics.

Economic planning was widespread in Europe at the time, and was also central in the development of Regional Science, in particular input-output models; as (Boyce 2004) recounts, Walter Isard worked closely with Wassily Leontief. Operational and planning motivations for applied economics were unquestioned, as economists in the post-war period saw their role, beyond educating young economists, as providing rational foundations for economic policy. It is worth noting that Jean Paelinck participated actively in the Association de Science Régionale De Langue Française, becoming president in 1973–1976. The first president of the association was François Perroux, who had founded it with Walter Isard in 1961 (Bailly, Derycke, and Torre 2012).

Until the 1980s, it was not at all unusual to publish original results in other languages than English. French spatial economic research, for example (Ponsard 1983), while making little impact in Anglophone countries, was widely used in teaching and research elsewhere (Billot and Thisse 1992). They contrast, though, the “word wizardry of François Perroux with the rigour of Claude Ponsard” (Billot and Thisse 1992), echoing the views expressed by (Drèze 1964) with regard to the work of Perroux. Even if we accept that “word wizardry” deserves more rigour and recasting in normative and empirically testable forms, it is also part of the context within which spatial econometrics came into being. A reading of (Perroux 1950) is worthwhile, because it not only gives the reader a vignette of the context in the post-war period, but also provides a discussion of economic space, as opposed to banal, unreflected space — mere position — that has largely disappeared from our considerations.

The title of the journal: Regional and Urban Economics, Operational Methods , founded by Jean Paelinck in 1971, and which was renamed as Regional Science and Urban Economics in 1975 (Boyce 2004), points to the perceived importance of “operational methods,” a version of the term “operational theory and method” used in the title of (Paelinck and Nijkamp 1975). Spatial econometrics does not seem to have come into being as a set of estimation techniques as such, as perhaps we might think today, but rather as an approach addressing open research questions both in space economy and in the enhancement of interregional models to be used in offering policy advice.

Were motivations of this kind common during the 1960s and early 1970s? Not only was the spread of Regional Science extensive and firmly established (Boyce 2004), but public bodies were concerned to regionalise economic measurement and policy advice (Graham and Romans 1971). In Britain, Environment and Planning was started in 1969 with Alan Wilson as founding editor and published by Pion; he was assistant director at the Centre for Environmental Studies at this time before moving to the University of Leeds. In a recently published lecture series, (Wilson 2012) cites (Paelinck and Nijkamp 1975) as giving principles for contributions from economics to urban and regional analysis (Wilson 2000). The papers presented at annual Regional Science meetings were published in a series by Pion; the first number in the series included contributions by (Granger 1969) and (Cliff and Ord 1969).

In a contribution to a panel session at the 2006 annual meeting of the American Association of Geographers (co-panelists Luc Anselin and Daniel Griffith), Keith Ord pointed to the continued relevance of Granger’s remarks at the meeting almost fourty years earlier (J. Keith Ord 2010); we will return to these concerns below. As noted by (R. S. Bivand 2008), communities of researchers working in and near mathematical and theoretical geography was more integrated in the pre-internet and pre-photocopier age than one might expect, with duplicated working papers prepared using stencils circulating rapidly between collaborating academic centres. Knowledge of the preliminary results of other researchers then fed through into rapid innovation in an exciting climate for those with access to these meetings and working papers.

There was considerable overlap between quantitative geography and regional science, so that work like (Cliff and Ord 1969) is cited by (Hordijk 1974), and was certainly known at the Netherlands Economic Institute several years earlier. Although it has not been possible to find out who participated in the August 1968 conference of the British and Irish Section of the Regional Science Association at which (Cliff and Ord 1969) was read, it was not unusual for members of other sections to be present, and to return home with bundles of duplicated papers. Up to the 1990s, presenters at conferences handed out copies of their papers, and conference participants posted home parcels of these hand-outs, indexed using the conference programme.

Leslie Hepple was among the more thorough scholars working on the underpinnings of spatial econometrics prior to the publication of (Paelinck and Klaassen 1979). His wide-ranging review (Hepple 1974) is cited by (Bartels and Hordijk 1977), again demonstrating the close links between those working in this field. We will be returning to the review paper, and to (Hepple 1976), which studies methods of estimation for spatial econometrics models in some depth, building on and extending (J. K. Ord 1975).

(Hepple 1974), like (Cliff and Ord 1973), saw no distinction between spatial statistics and the antecedents to spatial econometrics. Obviously, spatial econometrics was strongly influenced by the research tasks undertaken by regional and urban economists and regional scientists. As (Griffith and Paelinck 2007) point out, spatial statistics and spatial econometrics continue to share most topics of interests, with each also possessing shorter lists of topics that have been of less concern to the other. They advocate a “non-standard” spatial econometrics, which is inclusive to wider concerns. It seems appropriate in this context to mention the somewhat heterodox position taken by (D. P. McMillen 2010), who draws attention to the crucial issue of functional form, which he argues may well lie behind observed spatial autocorrelation; we will return to this in later chapters.

Implementing methods in econometrics

Having described some of the contextual issues suggesting how specific research concerns influenced how spatial econometrics came into being, it may now be helpful to turn to broader econometrics. It will become clearer that the research concerns of broader econometrics at that time, apart from spatial interdependence and interaction, were generally similar to those of spatial econometrics. Indeed, econometrics and the provision of national economic data for analysis, modelling, and for the provision of policy advice were intimately linked, as were mathematical economics and econometrics. The place of econometrics within economics, as (Sandmo 2011) shows, was a matter of some contention from the very beginning.

Both (Morgan 1990) and (Qin 1993) conclude their historical accounts of the beginnings of econometrics in pessimistic ways. Econometrics had begun by addressing a range of research topics, including the expression of economic theories in mathematical form, the building of operational econometric models, the exploration of testing and estimation techniques, and statistical data preparation. “By the 1950s the founding ideal of econometrics, the union of mathematical and statistical economics into a truly synthetic economics, had collapsed” (Morgan 1990). The flavour of (Paelinck and Klaassen 1979) seems somewhat anachronistic compared to late 1950s and early 1960s “textbook” econometrics. One might conclude that while spatial econometrics was aligned with Haavelmo’s probabilistic revolution (Qin 1993), it retained, like Haavelmo, adherence to the founding ideals [Bjerkholt (2007); hendry+johansen:12].

It appears from the research programme culminating in (Paelinck and Klaassen 1979) that spatial econometrics could be seen as “creative juggling in which theory and data came together to find out about the real world” (Morgan 1990). (Morgan 1990) and (Qin 1993) indicate that, from the 1950s, the depth of econometrics was flattened , with “[D]ata taken less seriously as a source of ideas and information for econometric models, and the theory-development role of applied econometrics was downgraded to the theory-testing role” (Morgan 1990). The history of econometric ideas is of assistance in illuminating key components of what has come to be the practice of spatial econometrics:

With Haavelmo, the profession seemed to have arrived at a consensus: a statistical model of plural economic causes and errors in the relations. But in the meantime, the old scientific ideal had also changed, from a deterministic to a probabilistic view of the way the world worked. This left the explanatory level of Haavelmo’s model open to new doubts. Was it really based on underlying random behaviour of the economic variables (as in contemporary evolutionary biology and quantum mechanics), or was it, after all, only a convenient way of formally dealing with inference in the non-experimental framework? (Morgan 1990).

The study of the history of econometrics has been actively promoted by David Hendry. It is desirable to motivate our choice of the term “applied spatial econometrics”; to do this we turn to (Hendry 2009), the introductory chapter in the “Applied Econometrics” volume of the Palgrave Handbook of Econometrics. If we can shed light on how “applied econometrics” is understood, then we may be able to provide adequate underpinnings for what is to follow. This, however, turns out to be hard to do, as Hendry challenges simple definitions:

At the superficial level, “Applied Econometrics” is “any application of econometrics, as distinct from theoretical econometrics. \(\ldots\) Some applied econometricians would include any applications involving analyses of”real economic data" by econometric methods, making “Applied econometrics” synonymous with empirical econometrics. However, such a view leads to demarcation difficulties from applied economics on the one hand, and applied statistics on the other. \(\ldots\) Outsiders might have thought that “Applied Econometrics” was just the application of econometrics to data, but that is definitely not so \(\ldots\) Rather, the notion of mutual penetration dominates, but as a one-way street. Economic theory comes first, almost mandatorially. Perhaps this just arises from a false view of science, namely that theory precedes evidence, \(\ldots\) (Hendry 2009).

He butresses his views using the history of econometrics to caricature empirical econometric research in the light of his view that applied economics has a “false” view of science:

\(\ldots\) cumulative critiques \(\ldots\) led to an almost monolithic approach to empirical econometric research: first postulate an individualistic, intertemporal optimization theory; next derive a model therefrom; third, find data with the same names as the theory variables; then select a recipe from the econometrics cookbook that appropriately blends the model and the data, or if necessary, develop another estimator; finally report the newly forged economic law. \(\ldots\) Instead of progress, we find fashions, cycles and “schools” in research. \(\ldots\) At about the same time that a priori theory-based econometrics became dominant, data measurement and quality issues were also relegated as a central component of empirical publications (Hendry 2009).

Hendry’s scepticism with regard to the practice of applied econometrics finds ample support in the proposal by (Anderson et al. 2008) that econometric results be acknowledged only to the extent that they are open for replication. They point out that theoretical results in economics (and econometrics) are much easier to check, and that referees routinely question formal proofs, because the relevant equations are included in journal submissions. Further contributions to the discussion on reproducible econometric research results have been made by (Koenker and Zeileis 2009) and (Yalta and Yalta 2010). As we will see later on, spatial econometrics is in the fortunate position of having relatively many open source software implementations.

In order to complete this brief review of the antecedents to spatial econometrics, it makes sense to point to descriptions of the development of econometric software. Both in a handbook chapter (Renfro 2009a), and in book form (Renfro 2009b), we can benefit from Charles Renfro’s experience and insights. He does admit to expressing distinct views, but they are backed both by experience and evidence. One concern is that econometric software development has been seen less and less as academic achievement, but rather as a practical, technical concern only:

To take an interest in programming and econometric software development would seem therefore to be the graveyard of any academic economist’s professional ambitions, justifiable only as a minimally diverting hobby, spoken of only to trusted colleagues. (Renfro 2009b)

He is not alone in drawing attention to the view that economics suffers from a skewed distribution of attention to the actual components of knowledge creation, placing theory firmly in first place, with commensurately much less weight given to measurement, data preparation and programming. As he says, applied economics researchers would benefit from a more even distribution of academic acknowledgement to those who, in terms of the division of labour within the discipline, develop software:

As a software developer, the econometrician who incorporates new theoretical econometric results may therefore be faced with the often-difficult task of not only evaluating the relevance, hence the operational validity of the theoretical solution, but also implementing these results in a way that is contextually meaningful. This operationally focused econometrician consequently not only needs to understand the theoretical advances that are made but also to exercise independent judgment in the numerical implementation of new techniques, for in fact neither are blueprints provided nor are the building blocks prefabricated. (Renfro 2009b)

Since this division of labour is arguably more uneven in economics than in other subjects, it leads to the restriction of research questions that applied economists can address to those that match estimation functions available in software implementations for which they have licences and which they know how to use:

One of the consequences of this specialization has been to introduce an element of user dependence, now grown to a sufficiently great degree that for many economists whichever set of econometric operations can be performed by their choice of program or programs has for them in effect become the universal set. (Renfro 2009b)

This deplorable situation has arisen over time, and certainly did not characterise the research context when spatial econometrics came into being. At that time, graduate students simply regarded learning to program, often in Fortran, as an essential part of their preparation as researchers. This meant that researchers were moch “closer” to their tools, and could adapt them to suit their research needs. Nowadays, few economists feel confident as programmers despite the fact that modern high-level languages such as Matlab, Python, or R are easy to learn and very flexible, and many econometric and statistical software applications offer scripting languages (such as SAS, Stata, SPSS and specifically econometric programs).

The links between methodological advance and the evolution of spatial economic theory are only touched upon in (Luc Anselin 2010) — in that sense, his review is concerned with theoretical spatial econometrics (statistical methods) rather than applied spatial econometrics (economic models). Over time, applied spatial econometrics has tended to become synonymous with regression modelling applied to spatial data where spatial autocorrelation and spatial heterogeneity in particular are present and need to be accommodated. Its treatment of spatial effects reflects the growing “legitimization of space and geography” (Luc Anselin 2010) in the quantitative social sciences more generally. But the subfield perhaps needs to be more than that if it is to justify its separate identity from spatial statistics and fully justify its “econometric” label. A close link with mainstream economic theory would seem essential in order to provide economic legitimacy to models (systems of equations) within which geography and spatial relationships have been, in economic terms, rigorously embedded (Fingleton 2000). (Haining 2014)

R’s sessionInfo()

sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Fedora 34 (Workstation Edition)
## 
## Matrix products: default
## BLAS:   /home/rsb/topics/R/R410-share/lib64/R/lib/libRblas.so
## LAPACK: /home/rsb/topics/R/R410-share/lib64/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] rgeoda_0.0.8-1 digest_0.6.27  ggplot2_3.3.3  tmap_3.3-1     spdep_1.1-8   
## [6] spData_0.3.8   sp_1.4-5       sf_0.9-9      
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-152       gmodels_2.18.1     RColorBrewer_1.1-2 tools_4.1.0       
##  [5] bslib_0.2.5.1      utf8_1.2.1         R6_2.5.0           KernSmooth_2.23-20
##  [9] DBI_1.1.1          colorspace_2.0-1   raster_3.4-10      withr_2.4.2       
## [13] tidyselect_1.1.1   gridExtra_2.3      leaflet_2.0.4.1    compiler_4.1.0    
## [17] leafem_0.1.3       expm_0.999-6       labeling_0.4.2     sass_0.4.0        
## [21] scales_1.1.1       classInt_0.4-3     proxy_0.4-25       stringr_1.4.0     
## [25] dbscan_1.1-8       rmarkdown_2.8      base64enc_0.1-3    dichromat_2.0-0   
## [29] pkgconfig_2.0.3    htmltools_0.5.1.1  highr_0.9          htmlwidgets_1.5.3 
## [33] rlang_0.4.11       jquerylib_0.1.4    generics_0.1.0     farver_2.1.0      
## [37] jsonlite_1.7.2     crosstalk_1.1.1    gtools_3.8.2       dplyr_1.0.6       
## [41] magrittr_2.0.1     Matrix_1.3-3       Rcpp_1.0.6         munsell_0.5.0     
## [45] fansi_0.4.2        abind_1.4-5        lifecycle_1.0.0    stringi_1.6.2     
## [49] leafsync_0.1.0     yaml_2.2.1         MASS_7.3-54        tmaptools_3.1-1   
## [53] grid_4.1.0         gdata_2.18.0       crayon_1.4.1       deldir_0.2-10     
## [57] lattice_0.20-44    stars_0.5-2        splines_4.1.0      knitr_1.33        
## [61] pillar_1.6.1       boot_1.3-28        codetools_0.2-18   LearnBayes_2.15.1 
## [65] XML_3.99-0.6       glue_1.4.2         evaluate_0.14      vctrs_0.3.8       
## [69] png_0.1-7          gtable_0.3.0       purrr_0.3.4        assertthat_0.2.1  
## [73] xfun_0.23          lwgeom_0.2-6       e1071_1.7-7        coda_0.19-4       
## [77] class_7.3-19       viridisLite_0.4.0  tibble_3.1.2       units_0.7-1       
## [81] ellipsis_0.3.2     spDataLarge_0.5.1
Anderson, Richard G., William H. Greene, B. D. McCullough, and H. D. Vinod. 2008. “The Role of Data/Code Archives in the Future of Economic Research.” Journal of Economic Methodology 15: 99–119.
Anselin, L. 1988. Spatial Econometrics: Methods and Models. Dordrecht: Kluwer Academic Publishers.
Anselin, L. 1995. “Local indicators of spatial association - LISA.” Geographical Analysis 27 (2): 93–115.
Anselin, L. 1996. “The Moran Scatterplot as an ESDA Tool to Assess Local Instability in Spatial Association.” In Spatial Analytical Perspectives on GIS, edited by M. M. Fischer, H. J. Scholten, and D. Unwin, 111–25. London: Taylor & Francis.
Anselin, Luc. 2010. “Thirty Years of Spatial Econometrics.” Papers in Regional Science 89: 3–25.
Assunção, R. M., and E. A. Reis. 1999. “A New Proposal to Adjust Moran’s \(I\) for Population Density.” Statistics in Medicine 18: 2147–62.
Bailly, Antoine, Pierre-Henri Derycke, and André Torre. 2012. “50 Ans de Science régionale Francophone.” l’Association de Science Régionale de Langue Française. http://www.asrdlf.org/telechargements/50ansASRDLF.pdf.
Bartels, Cornelis P. A., and Leen Hordijk. 1977. “On the Power of the Generalized Moran Contiguity Coefficient in Testing for Spatial Autocorrelation Among Regression Disturbances.” Regional Science and Urban Economics 7: 83–101.
Billot, Antoine, and Jacques-Francois Thisse. 1992. “Claude Ponsard (1927-1990): A Biographical Essay.” Annals of Regional Science 26: 191–98.
Bivand, R. S. 2008. “Implementing Representations of Space in Economic Geography.” Journal of Regional Science 48: 1–27.
Bivand, R. S., W. Müller, and M. Reder. 2009. “Power Calculations for Global and Local Moran’s \(I\).” Computational Statistics and Data Analysis 53: 2859–72.
Bivand, Roger S., and David W. S. Wong. 2018. “Comparing Implementations of Global and Local Indicators of Spatial Association.” TEST 27 (3): 716–48. https://doi.org/10.1007/s11749-018-0599-x.
Bjerkholt, Olav. 2007. “Writing ‘The Probability Approach’ with Nowhere to Go: Haavelmo in the United States, 1939–1944.” Econometric Theory 23: 775–837.
Boyce, David. 2004. “A Short History of the Field of Regional Science.” Papers in Regional Science 83: 31–57.
Cliff, A. D., and J. K. Ord. 1969. “The Problem of Spatial Autocorrelation.” In Studies in Regional Science, edited by A. J. Scott, 25–55. London Papers in Regional Science. London: Pion.
———. 1973. Spatial Autocorrelation. London: Pion.
———. 1981. Spatial Processes. London: Pion.
Drèze, Jacques H. 1964. “Some Postwar Contributions of French Economists to Theory and Public Policy: With Special Emphasis on Problems of Resource Allocation.” The American Economic Review 54 (4): 2–64.
Duncan, O. D., R. P. Cuzzort, and B. Duncan. 1961. Statistical Geography: Problems in Analyzing Areal Data. Glencoe, IL: Free Press.
Fingleton, B. 2000. “Spatial Econometrics, Economic Geography, Dynamics and Equilibrium: A ‘Third Way?’.” Environment and Planning A 32: 1481–98.
Fujita, Masahisa, Paul Krugman, and Anthony J. Venables. 1999. The Spatial Economy. Cambridge, MA: MIT Press.
Geary, R. C. 1954. “The Contiguity Ratio and Statistical Mapping.” The Incorporated Statistician 5: 115–45.
Getis, A., and J. K. Ord. 1992. “The Analysis of Spatial Association by the Use of Distance Statistics.” Geographical Analysis 24 (2): 189–206.
———. 1996. “Local Spatial Statistics: An Overview.” In Spatial Analysis: Modelling in a GIS Environment, edited by P. Longley and M Batty, 261–77. Cambridge: GeoInformation International.
Gibbons, Stephen, and Henry G. Overman. 2012. “Mostly Pointless Spatial Econometrics?” Journal of Regional Science 52: 172–91.
Graham, Robert E., and J. Thomas Romans. 1971. “Regional Economic Accounting in the u.s. Office of Business Economics — an Appraisal.” Review of Income and Wealth 17: 1–34.
Granger, Clive W. J. 1969. “Spatial Data and Time Series Analysis.” In Studies in Regional Science, edited by A. J. Scott, 1–25. London Papers in Regional Science 1. London: Pion.
Griffith, Daniel A., and Jean H. P. Paelinck. 2007. “An Equation by Any Other Name Is Still the Same: On Spatial Econometrics and Spatial Statistics.” Annals of Regional Science 41: 209–27.
Griliches, Zvi. 1967. “Distributed Lags: A Survey.” Econometrica 35: 16–49.
Haining, Robert. 2014. “Spatial Data and Statistical Methods: A Chronological Overview.” In Handbook of Regional Science, edited by M. M. Fischer and P. Nijkamp, 1277–94. Berlin: Springer-Verlag.
Hendry, David. 2009. “The Methodology of Empirical Econometric Modelling: Applied Econometrics Through the Looking-Glass.” In Palgrave Handbook on Econometrics: Applied Econometrics, edited by T. C. Mills and K. Patterson, 3–67. Basingstoke: Palgrave Macmillan.
Hepple, Leslie W. 1974. “The Impact of Stochastic Process Theory Upon Spatial Analysis in Human Geography.” Progress in Geography 6: 89–142.
———. 1976. “A Maximum Likelihood Model for Econometric Estimation with Spatial Series.” In Theory and Practice in Regional Science, edited by I. Masser, 90–104. London Papers in Regional Science. London: Pion.
Hordijk, Leen. 1974. “Spatial Correlation in the Disturbances of a Linear Interregional Model.” Regional and Urban Economics 4: 117–40.
Hordijk, Leen, and Jean H. P. Paelinck. 1976. “Some Principles and Results in Spatial Econometrics.” Recherches Economiques de Louvain 42: 175–97.
Klaassen, Leo H., Jean H. P. Paelinck, and Sjoerd Wagenaar. 1979. Spatial Systems. Farnborough: Saxon House.
Klein, L. R. 1958. “The Estimation of Distributed Lags.” Econometrica 26: 553–65.
Koenker, Roger, and Achim Zeileis. 2009. “On Reproducible Econometric Research.” Journal of Applied Econometrics 24: 833–47.
Koyck, L. M. 1954. Distributed Lags and Investment Analysis. Amsterdam: North-Holland.
Li, Xun, and Luc Anselin. 2021. Rgeoda: R Library for Spatial Data Analysis. https://CRAN.R-project.org/package=rgeoda.
McMillen, D. P. 2010. “Issues in Spatial Data Analysis.” Journal of Regional Science 50: 119–41.
McMillen, Daniel P. 2003. “Spatial Autocorrelation or Model Misspecification?” International Regional Science Review 26: 208–17.
Moran, P. A. P. 1948. “The Interpretation of Statistical Maps.” Journal of the Royal Statistical Society, Series B (Methodological) 10 (2): 243–51.
Morgan, Mary S. 1990. The History of Econometric Ideas. Cambridge: Cambridge University Press.
Ord, J. K. 1975. “Estimation Methods for Models of Spatial Interaction.” Journal of the American Statistical Association 70: 120–26.
Ord, J. K., and A. Getis. 2001. “Testing for Local Spatial Autocorrelation in the Presence of Global Autocorrelation.” Journal of Regional Science 41 (3): 411–32.
Ord, J. Keith. 2010. “Spatial Autocorrelation: A Statistician’s Reflections.” In Perspectives on Spatial Data Analysis, edited by L. Anselin and S. J. Rey, 165–80. Berlin: Springer-Verlag.
Paelinck, Jean H. P. 2012. “Specification and Identification in Spatial Econometric Models.” Estadística Española 54: 35–52.
———. 2013. “Some Challenges for Spatial Econometricians.” In Spatial Econometrics and Regional Economic Analysis, edited by Bogdan Suchecki, 292:11–20. Acta Universitas Lodziensis, Folia Oeconomica. Wydawnictwa Uniwersytetu Łódzkiego.
Paelinck, Jean H. P., and Leo H. Klaassen. 1979. Spatial Econometrics. Farnborough: Saxon House.
Paelinck, Jean H. P., and Peter Nijkamp. 1975. Operational Theory and Method in Regional Economics. Farnborough: Saxon House.
Perroux, François. 1950. “Economic Space: Theory and Applications.” Quarterly Journal of Economics 64: 89–104.
Ponsard, Claude. 1983. History of Spatial Economic Theory. Berlin: Springer-Verlag.
Qin, Duo. 1993. The Formation of Econometrics: A Historical Perspective. Oxford: Oxford University Press.
Renfro, Charles G. 2009a. “Econometric Software.” In Handbook of Computational Econometrics, edited by David A. Belsley and Erricos Kontoghiorghes, 1–53. Chichester, U.K: John Wiley & Sons.
———. 2009b. The Practice of Econometric Theory: An Examination of the Characteristics of Econometric Computation. Berlin: Springer-Verlag.
Sandmo, Agnar. 2011. Evolving Economics: A History of Economic Thought. Princeton: Princeton University Press.
Sauer, Jeffery, Taylor M Oshan, Sergio Rey, and Levi J Wolf. 2021. “On Null Hypotheses and Heteroskedasticity.” OSF Preprints. https://doi.org/10.31219/osf.io/ugkhp.
Schabenberger, O., and C. A. Gotway. 2005. Statistical Methods for Spatial Data Analysis. Boca Raton/London: Chapman & Hall/CRC.
Theil, H. 1964. “Some Developments of Economic Thought in the Netherlands.” The American Economic Review 54: pp. 34–55.
Tiefelsdorf, M. 2002. “The Saddlepoint Approximation of Moran’s I and Local Moran’s \({I}_i\) Reference Distributions and Their Numerical Evaluation.” Geographical Analysis 34: 187–206.
Upton, G., and B. Fingleton. 1985. Spatial Data Analysis by Example: Point Pattern and Qualitative Data. New York: Wiley.
Wilson, Alan G. 2000. Complex Spatial Systems: The Modelling Foundations of Urban and Regional Analysis. Harlow: Prentice Hall.
———. 2012. The Science of Cities and Regions: Lectures on Mathematical Model Design. Dordrecht: Springer.
Yalta, A.Talha, and A.Yasemin Yalta. 2010. “Should Economists Use Open Source Software for Doing Research?” Computational Economics 35: 371–94.