Section 12.1 Spatial regimes
Subsection 12.1.1 Fitting the constrained and the unconstrained models
To illustrate this solution, we will go back to the "ncovr" data we explored last week and use the same spatial weight matrix as in the previous chapter. If you have read the [Baller et al. (2001)] paper that we are in a way replicating here, you could see they decided that they needed to run separate models for the South and the North. Let’s explore this here.
data("ncovr")
# We coerce the sf object into a new sp object
ncovr_sp <- as(ncovr, "Spatial")
# Then we create a list of neighbours using 10 nearest neighbours criteria
xy <- coordinates(ncovr_sp)
nb_k10 <- knn2nb(knearneigh(xy, k = 10))
rwm <- nb2mat(nb_k10, style='W')
rwm <- mat2listw(rwm, style='W')
# remove redundant objects
rm(list = c("ncovr_sp", "xy", "nb_k10"))
We can start by running a model for the whole country to use as a baseline model:
fit1_70 <- lm(HR70 ~ RD70 + PS70 + MA70 + DV70 + UE70, data = ncovr)
summary(fit1_70)
## ## Call: ## lm(formula = HR70 ~ RD70 + PS70 + MA70 + DV70 + UE70, data = ncovr) ## ## Residuals: ## Min 1Q Median 3Q Max ## -21.12 -3.14 -1.02 1.90 69.75 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 11.1376 0.7602 14.65 <2e-16 *** ## RD70 4.1541 0.1147 36.23 <2e-16 *** ## PS70 1.1449 0.1165 9.82 <2e-16 *** ## MA70 -0.2080 0.0233 -8.94 <2e-16 *** ## DV70 1.4066 0.1064 13.22 <2e-16 *** ## UE70 -0.4268 0.0497 -8.59 <2e-16 *** ## --- ## Signif. codes: ## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 6.03 on 3079 degrees of freedom ## Multiple R-squared: 0.328, Adjusted R-squared: 0.327 ## F-statistic: 301 on 5 and 3079 DF, p-value: <2e-16
A simple visualisation of the residuals suggests the model may not be optimal for the entire study region:
ncovr$res_fit1 <- fit1_70$residuals
ggplot(ncovr, aes(x = res_fit1, colour = as.factor(SOUTH))) +
geom_density() +
xlab("Residuals") +
scale_colour_brewer(palette = "Set1", name = "Region (1 = South)") +
theme_bw()
To be able to fit these spatial regimes, that will allow for separate coefficients for counties in the North and those in the South, first we need to define a dummy for the Northern counties, and we create a vector of 1.
ncovr$NORTH <- recode(ncovr$SOUTH, `1` = 0, `0` = 1)
ncovr$ONE <- 1
Having created these new columns, we can fit a standard OLS model interacting the inputs with our dummies and setting the intercept to 0, so that each regime is allowed its own intercept:
fit2_70_sr <- lm(HR70 ~ 0 + (ONE + RD70 + PS70 + MA70 + DV70 + UE70):(SOUTH + NORTH),
data = ncovr)
summary(fit2_70_sr)
## ## Call: ## lm(formula = HR70 ~ 0 + (ONE + RD70 + PS70 + MA70 + DV70 + UE70):(SOUTH + ## NORTH), data = ncovr) ## ## Residuals: ## Min 1Q Median 3Q Max ## -17.47 -2.79 -0.76 1.71 65.60 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## ONE:SOUTH 15.0305 1.0714 14.03 <2e-16 *** ## ONE:NORTH 6.6485 1.1308 5.88 4.6e-09 *** ## RD70:SOUTH 3.2364 0.1884 17.18 <2e-16 *** ## RD70:NORTH 2.9523 0.2745 10.75 <2e-16 *** ## PS70:SOUTH 0.9953 0.2114 4.71 2.6e-06 *** ## PS70:NORTH 0.8428 0.1402 6.01 2.0e-09 *** ## MA70:SOUTH -0.1656 0.0338 -4.90 1.0e-06 *** ## MA70:NORTH -0.1749 0.0325 -5.37 8.2e-08 *** ## DV70:SOUTH 0.6271 0.1750 3.58 0.00034 *** ## DV70:NORTH 1.4263 0.1313 10.86 <2e-16 *** ## UE70:SOUTH -0.7916 0.0817 -9.69 <2e-16 *** ## UE70:NORTH -0.0293 0.0631 -0.46 0.64265 ## --- ## Signif. codes: ## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 5.81 on 3073 degrees of freedom ## Multiple R-squared: 0.645, Adjusted R-squared: 0.644 ## F-statistic: 465 on 12 and 3073 DF, p-value: <2e-16
We can see a higher intercept for the South, reflecting the higher level of homicide in these counties. Aside from this, the two more notable differences are the insignificant effect of unemployment in the North and the higher impact of the divorce rate in the North.
Subsection 12.1.2 Running the Chow test
The Chow test is often used in econometrics to determine whether the independent variables have different impacts on different sub-groups of the population. [Anselin (2007)] provides an implementation for this test that extracts the residuals and the degrees of freedom from the constrained (the simpler OLS model) and the unconstrained models (the less parsimonious spatial regime), and then computes the test based on the F distribution.
chow.test <- function(rest, unrest) {
er <- residuals(rest)
eu <- residuals(unrest)
er2 <- sum(er ^ 2)
eu2 <- sum(eu ^ 2)
k <- rest$rank
n2k <- rest$df.residual - k
c <- ((er2 - eu2) / k) / (eu2 / n2k)
pc <- pf(c, k, n2k, lower.tail = FALSE)
list(c, pc, k, n2k)
}
With the function in our environment we can then run the test, which will yield the test statistic, the p value, and two degrees of freedom.
chow.test(fit1_70,fit2_70_sr)
## [[1]] ## [1] 40.16 ## ## [[2]] ## [1] 2.674e-47 ## ## [[3]] ## [1] 6 ## ## [[4]] ## [1] 3073
How we interpret this test? When the Chow test is significant, as here, this provides evidence of spatial heterogeneity and in favour of fitting a spatial regime model, allowing for non-constant coefficients across the discrete subsets of data that we have specified.
Subsection 12.1.3 Spatial dependence with spatial heterogeneity
So far we have only adjusted for a form of spatial heterogeneity in our model, but we may still have to deal with spatial dependence. Here what we covered in Chapter 11, still applies. We can run the various specification tests and then select a model that adjusts for spatial dependence.
summary(lm.LMtests(fit2_70_sr, rwm, test = c("LMerr","LMlag","RLMerr","RLMlag")))
## Lagrange multiplier diagnostics for spatial ## dependence ## data: ## model: lm(formula = HR70 ~ 0 + (ONE + RD70 + ## PS70 + MA70 + DV70 + UE70):(SOUTH + NORTH), data ## = ncovr) ## weights: rwm ## ## statistic parameter p.value ## LMerr 227.20 1 <2e-16 *** ## LMlag 252.00 1 <2e-16 *** ## RLMerr 5.25 1 0.022 * ## RLMlag 30.05 1 4.2e-08 *** ## --- ## Signif. codes: ## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The Lagrange multiplier tests, according to Anselin’s suggestions for how to interpret these, point to a spatial lag model (since RLMlag is several orders more significant than the error alternative).
We could then fit this model using
sp:
fit3_70_sr_sar <- lagsarlm(HR70 ~ 0 +
(ONE +
RD70 +
PS70 +
MA70 +
DV70 +
UE70):(SOUTH + NORTH),
data = ncovr,
rwm)
summary(fit3_70_sr_sar)
## ## Call: ## lagsarlm(formula = HR70 ~ 0 + (ONE + RD70 + PS70 + MA70 + DV70 + ## UE70):(SOUTH + NORTH), data = ncovr, listw = rwm) ## ## Residuals: ## Min 1Q Median 3Q Max ## -17.15273 -2.66726 -0.67209 1.49219 65.86918 ## ## Type: lag ## Coefficients: (asymptotic standard errors) ## Estimate Std. Error z value Pr(>|z|) ## ONE:SOUTH 8.391260 1.103507 7.6042 2.864e-14 ## ONE:NORTH 4.903963 1.091437 4.4931 7.019e-06 ## RD70:SOUTH 2.355243 0.191445 12.3025 < 2.2e-16 ## RD70:NORTH 2.669956 0.266521 10.0178 < 2.2e-16 ## PS70:SOUTH 0.653640 0.203791 3.2074 0.0013394 ## PS70:NORTH 0.778913 0.135054 5.7674 8.050e-09 ## MA70:SOUTH -0.100568 0.032699 -3.0756 0.0021009 ## MA70:NORTH -0.138235 0.031328 -4.4126 1.022e-05 ## DV70:SOUTH 0.636867 0.168184 3.7867 0.0001527 ## DV70:NORTH 1.157985 0.128559 9.0074 < 2.2e-16 ## UE70:SOUTH -0.522353 0.079210 -6.5945 4.267e-11 ## UE70:NORTH -0.039129 0.060706 -0.6446 0.5192148 ## ## Rho: 0.3957, LR test value: 181.8, p-value: < 2.22e-16 ## Asymptotic standard error: 0.02846 ## z-value: 13.9, p-value: < 2.22e-16 ## Wald statistic: 193.3, p-value: < 2.22e-16 ## ## Log likelihood: -9710 for lag model ## ML residual variance (sigma squared): 31.21, (sigma: 5.587) ## Number of observations: 3085 ## Number of parameters estimated: 14 ## AIC: 19448, (AIC for lm: 19628) ## LM test for residual autocorrelation ## test value: 0.7147, p-value: 0.39789
We see that the AIC is better than for the non-spatial model and that there is no more residual autocorrelation. We could explore other specifications (i.e., spatial Durbin model) and the same caveats about interpreting coefficients in a spatial lag model we raised in Chapter 11 would apply here. If the Lagrange multiplier tests had suggested a spatial error model, we could also have fitted the model in the way we have also already covered simply adjusting the regression formula as in this section.
We can now also run a Chow test; only in this case, we need to adjust the previous one to account for spatial dependence. [Anselin (2007)] also provides for an implementation of this test in R with the following ad-hoc function:
spatialchow.test <- function(rest, unrest) {
lrest <- rest$LL
lunrest <- unrest$LL
k <- rest$parameters - 2
spchow <- -2.0 * (lrest - lunrest)
pchow <- pchisq(spchow, k, lower.tail = FALSE)
list(spchow, pchow, k)
}
This function will test a spatial lag model for the country as a whole with the spatial lag model with our data partitions. So first we need to estimate the appropriate lag model to be our constrained model:
fit4_70_sar <-
lagsarlm(HR70 ~ RD70 + PS70 + MA70 + DV70 + UE70,
data = ncovr,
rwm)
Which we can then use in the newly created function:
spatialchow.test(fit4_70_sar, fit3_70_sr_sar)
## [[1]] ## [,1] ## [1,] 58.59 ## ## [[2]] ## [,1] ## [1,] 8.712e-11 ## ## [[3]] ## [1] 6
The test is again highly significant providing evidence for the spatial lag model with the spatial regimes.
