10  Spatial Impacts

Required packages

pkgs <- c("sf", "mapview", "spdep", "spatialreg", "tmap", "viridisLite") # note: load spdep first, then spatialreg
lapply(pkgs, require, character.only = TRUE)

Session info

R version 4.5.1 (2025-06-13 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22631)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_United Kingdom.utf8 
[2] LC_CTYPE=English_United Kingdom.utf8   
[3] LC_MONETARY=English_United Kingdom.utf8
[4] LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.utf8    

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
[1] viridisLite_0.4.2 tmap_4.1          spatialreg_1.3-6 
[4] Matrix_1.7-3      spdep_1.3-13      spData_2.3.4     
[7] mapview_2.11.2    sf_1.0-21        

loaded via a namespace (and not attached):
 [1] xfun_0.52               raster_3.6-32          
 [3] htmlwidgets_1.6.4       lattice_0.22-7         
 [5] leaflet.providers_2.0.0 tools_4.5.1            
 [7] crosstalk_1.2.1         LearnBayes_2.15.1      
 [9] parallel_4.5.1          stats4_4.5.1           
[11] sandwich_3.1-1          proxy_0.4-27           
[13] KernSmooth_2.23-26      data.table_1.17.6      
[15] satellite_1.0.5         RColorBrewer_1.1-3     
[17] leaflet_2.2.2           lifecycle_1.0.4        
[19] compiler_4.5.1          farver_2.1.2           
[21] deldir_2.0-4            microbenchmark_1.5.0   
[23] terra_1.8-54            leafsync_0.1.0         
[25] codetools_0.2-20        leaflegend_1.2.1       
[27] stars_0.6-8             htmltools_0.5.8.1      
[29] class_7.3-23            MASS_7.3-65            
[31] classInt_0.4-11         lwgeom_0.2-14          
[33] wk_0.9.4                boot_1.3-31            
[35] abind_1.4-8             multcomp_1.4-28        
[37] nlme_3.1-168            digest_0.6.37          
[39] mvtnorm_1.3-3           splines_4.5.1          
[41] fastmap_1.2.0           grid_4.5.1             
[43] colorspace_2.1-1        cli_3.6.5              
[45] logger_0.4.0            magrittr_2.0.3         
[47] maptiles_0.10.0         base64enc_0.1-3        
[49] XML_3.99-0.18           dichromat_2.0-0.1      
[51] cols4all_0.8            survival_3.8-3         
[53] leafem_0.2.4            TH.data_1.1-3          
[55] e1071_1.7-16            scales_1.4.0           
[57] sp_2.2-0                rmarkdown_2.29         
[59] zoo_1.8-14              png_0.1-8              
[61] coda_0.19-4.1           evaluate_1.0.4         
[63] knitr_1.50              tmaptools_3.2          
[65] s2_1.1.9                rlang_1.1.6            
[67] Rcpp_1.0.14             glue_1.8.0             
[69] DBI_1.2.3               rstudioapi_0.17.1      
[71] jsonlite_2.0.0          R6_2.6.1               
[73] spacesXYZ_1.6-0         units_0.8-7            

Reload data from pervious session

load("_data/msoa2_spatial.RData")

10.1 Coefficient estimates `marginal’ effects

Warning

Do not interpret coefficients as marginal effects in SAR, SAC, and SDM!!

At first glance, the specifications presented above seem relatively similar in the way of modelling spatial effects. Yet, they differ in very important aspects.

First, models with an endogenous spatial term (SAR, SAC, and SDM) assume a very different spatial dependence structure than models with only exogenous spatial terms as SLX and SDEM specifications. While the first three assume global spatial dependence, the second two assume local spatial dependence (; ; ).

Second, the interpretation of the coefficients differs greatly between models with and without endogenous effects. This becomes apparent when considering the reduced form of the equations above. Exemplary using the SAR model, the reduced form is given by:

yρWy=Xβ+ε,(INρW)y=Xβ+ε,y=(INρW)1(Xβ+ε),

where IN is an N×N diagonal matrix (diagonal elements equal 1, 0 otherwise). This contains no spatially lagged dependent variable on the right-hand side.

If we want to interpret coefficient, we are usually in marginal or partial effects (the association between a unit change in X and Y). We obtain these effects by looking at the first derivative.

When taking the first derivative of the explanatory variable xk from the reduced form in (???) to interpret the partial effect of a unit change in variable xk on y, we receive

yxk=(INρW)1N×Nβk,

for each covariate k={1,2,...,K}. As can be seen, the partial derivative with respect to xk produces an N×N matrix, thereby representing the partial effect of each unit i onto the focal unit i itself and all other units .

Note that the diagonal elements of (INρW)1 are not zero anymore (as they are in W). Look at the following minimal example:

W~=(0101010101010101010101010),and normalized W=(00.500.500.3300.3300.3300.500.500.3300.3300.3300.500.50)

and

ρ=0.6,

then

ρW=(00.300.300.200.200.200.300.300.200.200.200.300.30).

If we want to get the total effect of X on Y we need to add the direct association within i and j and so on…

INρW=(1000001000001000001001001)(00.300.300.200.200.200.300.300.200.200.200.300.30)=(10.300.300.210.200.200.310.300.200.210.200.300.31).

And finally we take the inverse of that

(INρW)1=(10.300.300.210.200.200.310.300.200.210.200.300.31)1=(1.18750.468750.18750.468750.18750.31251.281250.31250.281250.31250.18750.468751.18750.468750.18750.31250.281250.31251.281250.31250.18750.468750.18750.468751.1875).

As you can see, (INρW)1 has >1: these are feedback loops. My X influences my Y directly, but my Y then influences my neigbour’s Y, which then influences my Y again (also also other neighbour’s Ys). Thus the influence of my X on my Y includes a spatial multiplier.

Check yourself:

I = diag(5)
rho = 0.6
W = matrix(c(0 , 0.5 , 0 , 0.5 , 0,
            1/3 , 0 , 1/3 , 0 , 1/3,
            0 , 0.5 , 0 , 0.5 , 0,
            1/3 , 0 , 1/3 , 0 , 1/3,
            0 , 0.5 , 0 , 0.5 , 0), ncol = 5, byrow = TRUE)

(IrW = I - rho*W)
     [,1] [,2] [,3] [,4] [,5]
[1,]  1.0 -0.3  0.0 -0.3  0.0
[2,] -0.2  1.0 -0.2  0.0 -0.2
[3,]  0.0 -0.3  1.0 -0.3  0.0
[4,] -0.2  0.0 -0.2  1.0 -0.2
[5,]  0.0 -0.3  0.0 -0.3  1.0
# (I - rho*W)^-1
(M = solve(IrW))
       [,1]    [,2]   [,3]    [,4]   [,5]
[1,] 1.1875 0.46875 0.1875 0.46875 0.1875
[2,] 0.3125 1.28125 0.3125 0.28125 0.3125
[3,] 0.1875 0.46875 1.1875 0.46875 0.1875
[4,] 0.3125 0.28125 0.3125 1.28125 0.3125
[5,] 0.1875 0.46875 0.1875 0.46875 1.1875

The diagonal elements of M indicate how each unit i influences itself (change of xi on change of yi), and each off-diagonal elements in column j represents the effect of j on each other unit i (change of xj on change of yi).

(1.18750.468750.18750.468750.18750.31251.281250.31250.281250.31250.18750.468751.18750.468750.18750.31250.281250.31251.281250.31250.18750.468750.18750.468751.1875).

For instance, W12 indicates that unit 2 has an influence of 0.46875 on unit 1. On the other hand, W53 indicates that unit 3 has an influence of magnitude 0.1875 on unit 5.

Question

Why does unit 3 have any effect o unit 5? According to W those two units are no neighbours w53=0!

10.2 Global and local spillovers

The kind of indirect spillover effects in SAR, SAC, and SDM models differs from the kind of indirect spillover effects in SLX and SDEM models: while the first three specifications represent global spillover effects, the latter three represent local spillover effects (; ; ).

10.2.1 Local spillovers

In case of SLX and SDEM the spatial spillover effects can be interpreted as the effect of a one unit change of xk in the spatially weighted neighbouring observations on the dependent variable of the focal unit: the weighted average among neighbours; when using a row-normalised contiguity weights matrix, Wxk is the mean value of xk in the neighbouring units.

Assume we have k=2 covariates, then

WN×NXN×2θ2×1=(00.500.500.3300.3300.3300.500.500.3300.3300.3300.500.50)(3100414012007705250)(θ1θ2)=(61053190610531906105)(θ1θ2)

X <- cbind(x1 = c(3,4,1,8,5),
           x2 = c(100,140,200,70,270))
(WX <-  W %*% X)
     x1  x2
[1,]  6 105
[2,]  3 190
[3,]  6 105
[4,]  3 190
[5,]  6 105

Thus, only direct neighbours – as defined in W – contribute to those local spillover effects. The θ^ coefficients only estimate how my direct neighbour’s X values influence my own outcome y.

There are no higher order neighbours involved (as long as we do not model them), nor are there any feedback loops due to interdependence.

10.2.2 Global spillovers

In contrast, spillover effects in SAR, SAC, and SDM models do not only include direct neighbours but also neighbours of neighbours (second order neighbours) and further higher-order neighbours. This can be seen by rewriting the inverse (INρW)1 as power series:A power series of k=0Wk converges to (IW)1 if the maximum absolute eigenvalue of W<1, which is ensured by standardizing W.}

(INρW)1βk=(IN+ρW+ρ2W2+ρ3W3+...)βk=(IN+h=1ρhWh)βk,

where the identity matrix represents the direct effects and the sum represents the first and higher order indirect effects and the above mentioned feedback loops. This implies that a change in one unit i does not only affect the direct neighbours but passes through the whole system towards higher-order neighbours, where the impact declines with distance within the neighbouring system. Global indirect impacts thus are `multiplied’ by influencing direct neighbours as specified in W and indirect neighbours not connected according to W, with additional feedback loops between those neighbours.

(INN×Nρ=^0.6WN×N)1N×Nβk=(1.18750.468750.18750.468750.18750.31251.281250.31250.281250.31250.18750.468751.18750.468750.18750.31250.281250.31251.281250.31250.18750.468750.18750.468751.1875)(β1+β2).

All diagonal elements of diag(W)=wii=0. However, diagonal elements of higher order neighbours are not zero diag(W2)=diag(WW)0.

Intuitively, ρW only represents the effects between direct neighbours (and the focal unit is not a neighbour of the focal unit itself), whereas ρ2W2 contains the effects of second order neighbours, where the focal unit is a second order neighbour of the focal unit itself. Thus, (INρW)1βk includes feedback effects from ρ2W2 on (they are part of the direct impacts according to the summary measures below). This is way the diagonal above 1.

In consequence, local and global spillover effects represent two distinct kinds of spatial spillover effects (). The interpretation of local spillover effects is straightforward: it represents the effect of all neighbours as defined by W (the average over all neighbours in case of a row-normalised weights matrix).

For instance, the environmental quality in the focal unit itself but also in neighbouring units could influence the attractiveness of a district and its house prices. In this example it seems reasonable to assume that we have local spillover effects: only the environmental quality in directly contiguous units (e.g. in walking distance) is relevant for estimating the house prices.

In contrast, interpreting global spillover effects can be a bit more difficult. Intuitively, the global spillover effects can be seen as a kind of diffusion process. For example, an exogenous event might increase the house prices in one district of a city, thus leading to an adaptation of house prices in neighbouring districts, which then leads to further adaptations in other units (the neighbours of the neighbours), thereby globally diffusing the effect of the exogenous event due to the endogenous term.

Yet, those processes happen over time. In a cross-sectional framework, the global spillover effects are hard to interpret. Anselin () proposes an interpretation as an equilibrium outcome, where the partial impact represents an estimate of how this long-run equilibrium would change due to a change in xk ().

10.3 Summary impact measures

Note that the derivative in SAR, SAC, and SDM is a N×N matrix, returning individual effects of each unit on each other unit, differentiated in direct, indirect, and total impacts.

(INρW)1β=(1.18750.468750.18750.468750.18750.31251.281250.31250.281250.31250.18750.468751.18750.468750.18750.31250.281250.31251.281250.31250.18750.468750.18750.468751.1875)β

However, the individual effects (how i influences j) mainly vary because of variation in W.

Do not interpret these as “estimated” individual impacts

We estimate two scalar parameters in a SAR model: β for the direct coefficient and rho for the auto-regressive parameter.

All variation in the effects matrix (INρW)1 comes from the relationship in W which we have given a-priori!

Since reporting the individual partial effects is usually not of interest, LeSage and Pace () proposed to average over these effect matrices. While the average diagonal elements of the effects matrix (INρW)1 represent the so called direct impacts of variable xk, the average column-sums of the off-diagonal elements represent the so called indirect impacts (or spatial spillover effects).

direct impacts refer to an average effect of a unit change in xi on yi, and the indirect (spillover) impacts indicate how a change in xi, on average, influences all neighbouring units yj.

Though previous literature (; ) has established the notation of direct and indirect impacts, it is important to note that also the direct impacts comprise a spatial `multiplier’ component if we specify an endogenous lagged depended variable, as a change in xi influences yi, which influences yj, which in turn influences yi.

Usually, one should use summary measures to report effects in spatial models (). Halleck Vega and Elhorst () provide a nice summary of the impacts for each model:

Model Direct Impacts Indirect Impacts type
OLS/SEM βk
SAR/SAC Diagonal elements of (IρW)1βk Off-diagonal elements of (IρW)1βk global
SLX/SDEM βk θk local
SDM Diagonal elements of (IρW)1[βk+Wθk] Off-diagonal elements of (IρW)1[βk+Wθk] global

(INρW)1β=(1.18750.468750.18750.468750.18750.31251.281250.31250.281250.31250.18750.468751.18750.468750.18750.31250.281250.31251.281250.31250.18750.468750.18750.468751.1875)β

The different indirect effects / spatial effects mean conceptually different things:

  • Global spillover effects: SAR, SAC, SDM

  • Local spillover effects: SLX, SDEM

Common ratio between direct and indirect impacts in SAR and SAC

Note that impacts in SAR only estimate one single spatial multiplier coefficient. Thus direct and indirect impacts are bound to a common ratio, say ϕ, across all covariates.

if β1direct=ϕβ1indirect, then β2direct=ϕβ2indirect, βkdirect=ϕβkindirect.

We can calculate these impacts using impacts() with simulated distributions, e.g. for the SAR model:

mod_1.sar.imp <- impacts(mod_1.sar, listw = queens.lw, R = 300)
summary(mod_1.sar.imp, zstats = TRUE, short = TRUE)
Impact measures (lag, exact):
                  Direct     Indirect        Total
log(no2)     0.447853184  0.754466618  1.202319802
log(POPDEN) -0.062973027 -0.106086209 -0.169059236
per_mixed    0.020884672  0.035182931  0.056067603
per_asian   -0.002575602 -0.004338934 -0.006914536
per_black   -0.014253206 -0.024011369 -0.038264575
per_other   -0.001820705 -0.003067212 -0.004887917
========================================================
Simulation results ( variance matrix):
========================================================
Simulated standard errors
                  Direct     Indirect       Total
log(no2)    0.0441693298 0.0877111209 0.119711919
log(POPDEN) 0.0149810657 0.0277539086 0.041947946
per_mixed   0.0068608095 0.0122281335 0.018911333
per_asian   0.0005286191 0.0008898791 0.001381428
per_black   0.0010155041 0.0021948439 0.002705126
per_other   0.0030080857 0.0051864728 0.008183300

Simulated z-values:
                 Direct   Indirect       Total
log(no2)     10.1690813   8.616034  10.0648584
log(POPDEN)  -4.2088999  -3.835768  -4.0409931
per_mixed     2.9428632   2.788271   2.8705422
per_asian    -4.9651579  -4.947168  -5.0868084
per_black   -13.9952569 -10.891794 -14.0910349
per_other    -0.6753702  -0.668458  -0.6719186

Simulated p-values:
            Direct     Indirect   Total     
log(no2)    < 2.22e-16 < 2.22e-16 < 2.22e-16
log(POPDEN) 2.5662e-05 0.00012517 5.3225e-05
per_mixed   0.0032519  0.00529901 0.0040977 
per_asian   6.8645e-07 7.5301e-07 3.6414e-07
per_black   < 2.22e-16 < 2.22e-16 < 2.22e-16
per_other   0.4994406  0.50384127 0.5016356 
# Alternative with traces (better for large W)
W <- as(queens.lw, "CsparseMatrix")
trMatc <- trW(W, type = "mult",
              m = 30) # number of powers
mod_1.sar.imp2 <- impacts(mod_1.sar, 
                          tr = trMatc, # trace instead of listw
                          R = 300, 
                          Q = 30) # number of power series used for approximation
summary(mod_1.sar.imp2, zstats = TRUE, short = TRUE)
Impact measures (lag, trace):
                  Direct     Indirect        Total
log(no2)     0.447853101  0.754459497  1.202312598
log(POPDEN) -0.062973015 -0.106085208 -0.169058223
per_mixed    0.020884668  0.035182599  0.056067267
per_asian   -0.002575601 -0.004338893 -0.006914494
per_black   -0.014253203 -0.024011142 -0.038264346
per_other   -0.001820704 -0.003067183 -0.004887888
========================================================
Simulation results ( variance matrix):
========================================================
Simulated standard errors
                 Direct     Indirect       Total
log(no2)    0.050264654 0.1010752967 0.138217388
log(POPDEN) 0.015213135 0.0291848258 0.043543021
per_mixed   0.006790766 0.0122043557 0.018787624
per_asian   0.000562782 0.0009300614 0.001445988
per_black   0.001038811 0.0022797630 0.002713717
per_other   0.003095462 0.0053549649 0.008437334

Simulated z-values:
                 Direct    Indirect       Total
log(no2)      8.9319119   7.5616213   8.7778577
log(POPDEN)  -4.1796882  -3.7272331  -3.9584946
per_mixed     3.0032711   2.8562081   2.9409088
per_asian    -4.6615569  -4.7835265  -4.8910600
per_black   -13.7029911 -10.6234400 -14.1701374
per_other    -0.4453985  -0.4552648  -0.4523516

Simulated p-values:
            Direct     Indirect   Total     
log(no2)    < 2.22e-16 3.9746e-14 < 2.22e-16
log(POPDEN) 2.9191e-05 0.00019359 7.5424e-05
per_mixed   0.0026709  0.00428734 0.0032725 
per_asian   3.1383e-06 1.7225e-06 1.0029e-06
per_black   < 2.22e-16 < 2.22e-16 < 2.22e-16
per_other   0.6560318  0.64891877 0.6510157 

The indirect effects in SAR, SAC, and SDM refer to global spillover effects. This means a change of x in the focal units flows through the entire system of neighbours (direct nieightbours, neighbours of neighbours, …) influencing ‘their y’. One can think of this as diffusion or a change in a long-term equilibrium.

If Log NO2 increases by one unit, this increases the house price in the focal unit by 0.448 units. Overall, a one unit change in log NO2 increases the house prices in the entire neighbourhood system (direct and higher order neighbours) by 0.754.

For SLX models, nothing is gained from computing the impacts, as they equal the coefficients. Again, it’s the effects of direct neighbours only.

print(impacts(mod_1.slx, listw = queens.lw))
Impact measures (SlX, glht):
                  Direct     Indirect        Total
log(no2)    -0.440727458  0.993602103  0.552874645
log(POPDEN) -0.076839828  0.113262218  0.036422390
per_mixed   -0.033042221  0.126068686  0.093026466
per_asian   -0.002380698 -0.003828126 -0.006208824
per_black   -0.016229407 -0.018053503 -0.034282910
per_other   -0.020391354  0.048139008  0.027747654

10.4 Examples

Boillat, Ceddia, and Bottazzi ()

The paper investigates the effects of protected areas and various land tenure regimes on deforestation and possible spillover effects in Bolivia, a global tropical deforestation hotspot.

Protected areas – which in Bolivia are all based on co-management schemes - also protect forests in adjacent areas, showing an indirect protective spillover effect. Indigenous lands however only have direct forest protection effects.

Fischer et al. ()

The focus of this paper is on the role of human capital in explaining labor productivity variation among 198 European regions within a regression framework.

A ceteris paribus increase in the level of human capital is found to have a significant and positive direct impact. But this positive direct impact is offset by a significant and negative indirect (spillover) impact leading to a total impact that is not significantly different from zero.

The intuition here arises from the notion that it is relative regional advantages in human capital that matter most for labor productivity, so changing human capital across all regions should have little or no total impact on (average) labor productivity levels.

Rüttenauer ()

This study investigates the presence of environmental inequality in Germany - the connection between the presence of foreign-minority population and objectively measured industrial pollution.

Results reveal that the share of minorities within a census cell indeed positively correlates with the exposure to industrial pollution. Furthermore, spatial spillover effects are highly relevant: the characteristics of the neighbouring spatial units matter in predicting the amount of pollution. Especially within urban areas, clusters of high minority neighbourhoods are affected by high levels of environmental pollution.