Section 11.6 Fitting and interpreting a spatially lagged model
We saw that in the previous case a spatial error model was perhaps more appropriate, but we also saw earlier in one of the simpler models that there was an interaction between the variable indicating whether the county is in the South and one of the explanatory variables. This is suggestive of the kind of spatial heterogeneity we were discussing at the outset of the chapter. When the relationship of the explanatory variables is not constant across the study area, we need to take this into account. Fitting the same model for the whole country may not be appropriate and there are so many interaction terms you can include before things get difficult to interpret. A fairly crude approach, instead of adding multiple interaction terms, is to partition the data. We could, in this example, partition the data into Southern and Northern counties and run separate models for this using two subsets of the data. For the next illustration, we are going to focus on the Southern counties and in a different decade, in which the lag model (as we will see) is more appropriate.
# Let's look at southern counties
ncovr_s_sf <- subset(ncovr, SOUTH == 1)
# We coerce the sf object into a new sp object
ncovr_s_sp <- as(ncovr_s_sf, "Spatial")
# Then we create a list of neighbours using 10 nearest neighbours criteria
xy_s <- coordinates(ncovr_s_sp)
nb_s_k10 <- knn2nb(knearneigh(xy_s, k = 10))
rwm_s <- nb2mat(nb_s_k10, style='W')
rwm_s <- mat2listw(rwm_s, style='W')
# We create new equation, this time excluding the southern dummy
eq2_70_S <- HR70 ~ RD70 + PS70 + MA70 + DV70 + UE70
fit4_70 <- lm(eq2_70_S, data=ncovr_s_sf)
summary(lm.LMtests(fit4_70, rwm_s, test = c("LMerr","LMlag","RLMerr","RLMlag")))
## Lagrange multiplier diagnostics for spatial ## dependence ## data: ## model: lm(formula = eq2_70_S, data = ncovr_s_sf) ## weights: rwm_s ## ## statistic parameter p.value ## LMerr 1.47e+02 1 <2e-16 *** ## LMlag 1.75e+02 1 <2e-16 *** ## RLMerr 1.55e-03 1 0.97 ## RLMlag 2.86e+01 1 9e-08 *** ## --- ## Signif. codes: ## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Maximum Likelihood (ML) estimation of the spatial lag model is carried out with the
spatialreg::lagsarlm() function. The required arguments are again a regression "formula", a data set and a listw spatial weights object. The default method uses Ord’s eigenvalue decomposition of the spatial weights matrix.
fit4_70_sar <- lagsarlm(eq2_70_S, data=ncovr_s_sf, rwm_s)
summary(fit4_70_sar)
## ## Call: ## lagsarlm(formula = eq2_70_S, data = ncovr_s_sf, listw = rwm_s) ## ## Residuals: ## Min 1Q Median 3Q Max ## -17.13387 -4.09583 -0.95433 2.85297 65.93842 ## ## Type: lag ## Coefficients: (asymptotic standard errors) ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 7.437000 1.413763 5.2604 1.437e-07 ## RD70 2.247458 0.244104 9.2070 < 2.2e-16 ## PS70 0.603290 0.256056 2.3561 0.018469 ## MA70 -0.091496 0.041144 -2.2238 0.026160 ## DV70 0.633473 0.211198 2.9994 0.002705 ## UE70 -0.485272 0.099684 -4.8681 1.127e-06 ## ## Rho: 0.4523, LR test value: 114.8, p-value: < 2.22e-16 ## Asymptotic standard error: 0.04111 ## z-value: 11, p-value: < 2.22e-16 ## Wald statistic: 121, p-value: < 2.22e-16 ## ## Log likelihood: -4769 for lag model ## ML residual variance (sigma squared): 49.18, (sigma: 7.013) ## Number of observations: 1412 ## Number of parameters estimated: 8 ## AIC: 9555, (AIC for lm: 9668) ## LM test for residual autocorrelation ## test value: 7.662, p-value: 0.0056391
Remember what we said earlier in the spatial lag model that we are simply adding as an additional explanatory variable the values of \(y\) in the surrounding area. What we mean by "surrounding" will be defined by our spatial weight matrix. It’s important to emphasise again that one has to think very carefully and explore appropriate definitions of "surrounding".
In our spatial lag model, you will notice that there is a new term \(Rho\text{.}\) What is this? This is our spatial lag. It is a variable that measures the homicide rate in the counties that are defined as surrounding each county in our spatial weight matrix. We are simply using this variable as an additional explanatory variable to our model, so that we can appropriately take into account the spatial clustering detected by our Moran’s I test. You will notice that the estimated coefficient for this term is both positive and statistically significant. In other words, when the homicide rate in surrounding areas increases, so does the homicide rate in each country, even when we adjust for the other explanatory variables in our model.
You also see at the bottom further tests for spatial dependence, a likelihood ratio test. This is not a test for residual spatial autocorrelation after we introduce our spatial lag. What you want is for this test to be significant because in essence is further evidence that the spatial lag model is a good fit.
How about the coefficients? It may be tempting to look at the regression coefficients for the other explanatory variables for the original OLS model and compare them to those in the spatial lag model. But you should be careful when doing this. Their meaning now has changed:
Interpreting the substantive effects of each predictor in a spatial lag model is much more complex than in a nonspatial model (or in a spatial error model) because of the presence of the spatial multiplier that links the independent variables to the dependent. In the nonspatial model, it does not matter which unit is experiencing the change on the independent variable. The effect in the dependent variable of a change in the value of an independent variable is constant across all observations ([179], p. 107).
Remember we say, when interpreting a regression coefficient for variable \(x_i\text{,}\) that they indicate how much \(y\) goes up or down for every one unit increase in \(x_i\) when holding all other variables in the model constant. In our example, for the non-spatial model this "effect" is the same for every county in our dataset. But in the spatial lag model, things are not the same. We cannot interpret the regression coefficients for the substantive predictors in the same way because the "substantive effects of the independent variables vary by observation as a result of the different neighbours for each unit in the data" [179].
In the standard OLS regression model, the coefficients for any of the explanatory variables measure the absolute "impact" of these variables. It is a simpler scenario. We look at the "effect" of \(x\) in \(y\) within each county. So \(x\) in county A affects \(y\) in county A. In the spatial lag model there are two components to how \(x\) affects \(y\text{.}\) The predictor \(x\) affects \(y\) within each county directly, but remember we are also including the spatial lag, the measure of \(y\) in the surrounding counties (call them B, C, and D). By virtue of including the spatial lag (a measure of \(y\) in county B, C, and D) we are indirectly incorporating as well the "effects" that \(x\) has on \(y\) in counties B, C, and D. So the "effect" of a covariate (independent variable) is the sum of two particular "effects": a direct, local effect of the covariate in that unit, and an indirect, spillover "effect" due to the spatial lag.
In other words, in the spatial lag model, the coefficients only focus on the "short-run impact" of \(x_i\) on \(y_i\text{,}\) rather than the net "effect". As [202] explain, "Since the value of \(y_i\) will influence the level of homicide in other counties \(y_j\) and these \(y_j\text{,}\) in turn, feedback on to \(y_i\text{,}\) we need to take into account the additional effects that the short impact of \(x_i\) exerts on \(y_i\) through its impact on the level of homicide in other counties". You can still read the coefficients in the same way but need to keep in mind that they are not measuring the net "effect". Part of their "effect" will be captured by the spatial lag. Yet, you may still want to have a look at whether things change dramatically, particularly in terms of their significance (which is not the case in this example).
In sum, this implies that a change in the \(i\)th region’s predictor can affect the \(j\)th region’s outcome. We have two situations: (a) the direct impact of an observation’s predictor on its own outcome, and (b) the indirect impact of an observation’s neighbour’s predictor on its outcome. This leads to three quantities that we want to know:
-
Average Direct Impact, which is similar to a traditional interpretation
-
Average Total Impact, which would be the total of direct and indirect impacts of a predictor on one’s outcome
-
Average Indirect Impact, which would be the average impact of one’s neighbours on one’s outcome
These quantities can be found using the
impacts() function in the spdep library. We follow the example that converts the spatial weight matrix into a "sparse" matrix and power it up using the trW() function. This follows the approximation methods described in [204]. Here, we use Monte Carlo simulation to obtain simulated distributions of the various impacts.
W <- as(rwm_s, "CsparseMatrix")
trMC <- trW(W, type="MC")
fit4_70_sar.im <- impacts(fit4_70_sar, tr=trMC, R=100)
summary(fit4_70_sar.im, zstats=TRUE, short=TRUE)
## Impact measures (lag, trace): ## Direct Indirect Total ## RD70 2.30369 1.79942 4.1031 ## PS70 0.61838 0.48302 1.1014 ## MA70 -0.09379 -0.07326 -0.1670 ## DV70 0.64932 0.50719 1.1565 ## UE70 -0.49741 -0.38853 -0.8859 ## ======================================================== ## Simulation results ( variance matrix): ## ======================================================== ## Simulated standard errors ## Direct Indirect Total ## RD70 0.23328 0.27135 0.40515 ## PS70 0.27102 0.22932 0.49148 ## MA70 0.04071 0.03240 0.07195 ## DV70 0.21484 0.19467 0.39866 ## UE70 0.10801 0.09425 0.18968 ## ## Simulated z-values: ## Direct Indirect Total ## RD70 9.776 6.664 10.092 ## PS70 2.336 2.199 2.314 ## MA70 -2.188 -2.165 -2.213 ## DV70 3.147 2.775 3.051 ## UE70 -4.541 -4.118 -4.632 ## ## Simulated p-values: ## Direct Indirect Total ## RD70 < 2e-16 2.7e-11 < 2e-16 ## PS70 0.0195 0.0279 0.0207 ## MA70 0.0287 0.0304 0.0269 ## DV70 0.0017 0.0055 0.0023 ## UE70 5.6e-06 3.8e-05 3.6e-06
We see that all the variables have significant direct, indirect and total effects. You may want to have a look at how things differ when you just run a non-spatial model.
