knitr::opts_chunk$set(error = TRUE)
library(nortest)
library(car)
library(lmtest)
library(classInt)
library(spdep)
library(maptools)
library(htmlTable)
library(RColorBrewer)
library(rgdal)

Spatial Regression Models

This lecture builds off the previous lecture on the Spatially Autoregressive Model (SAR) with either a lag or error specification. The lag model is written: \(Y= \rho W Y + X '\beta +e\)

Where Y is the dependent variable, X is the matrix of independent variables, \(\beta\) is the vector of regression parameters to be estimated from the data, \(\rho\) is the autoregressive coefficient, which tells us how strong the resemblance is, on average, between \(Y_i\) and it’s neighbors. The matrix W is the spatial weight matrix, describing the spatial network structure of the observations, like we described in the ESDA lecture.

In the lag model, we are specifying the spatial component on the dependent variable. This leads to a spatial filtering of the variable, where they are averaged over the surrounding neighborhood defined in W, called the spatially lagged variable. In R we use the spdep package, and the lagsarlm() function to fit this model.

The error model says that the autocorrelation is not in the outcome itself, but instead, any autocorrelation is attributable to there being missing spatial covariates in the data. If these spatially patterned covariates could be measures, the tne autocorrelation would be 0. This model is written:

\(Y= X' \beta +e\)

\(e=\lambda W e + v\)

This model, in effect, controls for the nuisance of correlated errors in the data that are attributable to an inherently spatial process, or to spatial autocorrelation in the measurement errors of the measured and possibly unmeasured variables in the model. This model is estimated in R using errorsarlm() in the spdep library.

Examination of Model Specification

To some degree, both of the SAR specifications allow us to model spatial dependence in the data. The primary difference between them is where we model said dependence.

The lag model says that the dependence affects the dependent variable only, we can liken this to a diffusion scenario, where your neighbors have a diffusive effect on you.

The error model says that dependence affects the residuals only. We can liken this to the missing spatially dependent covariate situation, where, if only we could measure another really important spatially associated predictor, we could account for the spatial dependence. But alas, we cannot, and we instead model dependence in our errors.

These are inherently two completely different ways to think about specifying a model, and we should really make our decision based upon how we think our process of interest operates.

That being said, this way of thinking isn’t necessarily popular among practitioners. Most practitioners want the best fitting model, ‘nuff said. So methods have been developed that test for alternate model specifications, to see which kind of model best summarizes the observed variation in the dependent variable and the spatial dependence.

More exotic types of spatial dependence

Spatial Durbin Model Another form of a spatial lag model is the Spatial Durbin Model (SDM). This model is an extension of the ordinary lag or error model that includes spatially lagged independent variables. If you remember, one issue that commonly occures with the lag model, is that we often have residual autocorrelation in the model. This autocorrelation could be attributable to a missing spatial covariate. We can get a kind of spatial covariate by lagging the predictor variables in the model using W. This model can be written:

\(Y= \rho W Y + X '\beta + W X \theta + e\)

Where, the \(\theta\) parameter vector are now the regression coefficients for the lagged predictor variables. We can also include the lagged predictors in an error model, which gives us the Durbin Error Model (DEM):

\(Y= X '\beta + W X \theta + e\)

\(e=\lambda W e + v\)

Generally, the spatial Durbin model is preferred to the ordinary error model, because we can include the “unspecified” spatial covariates from the error model into the Durbin model via the lagged predictor variables.

Spatially Autoregressive Moving Average Model Futher extensions of these models include dependence on both the outcome and the error process. Two models are described in LeSage and Pace. The Spatial Autocorrelation Model, or SAC model and the Spatially autoregressive moving average model (SARMA model). The SAC model is:

\(Y= \rho W_1 Y + X '\beta + e\)

\(e=\theta W_2 e + v\)

\(Y= (I_n - \rho W_1)^{-1} X '\beta + (I_n - \rho W_1)^{-1} (I_n - \theta W_2)^{-1} e\)

Where, you can potentially have two different spatial weight matrices, \(W_1\) and \(W_2\). Here, the lagged error term is taken over all orders of neighbors, leading to a more global error process, while the SARMA model has form:

\(Y= \rho W_1 Y + X '\beta + u\)

\(u=(I_n - \theta W_2) e\)

\(e \sim N(0, \sigma^2 I_n)\)

\(Y= (I_n - \rho W_1)^{-1} X '\beta + (I_n - \rho W_1)^{-1} (I_n - \theta W_2) e\)

which gives a “locally” weighted moving average to the residuals, which will avereage the residuals only in the local neighborhood, instead of over all neighbor orders.

Fitting these models in R can be done in the spdep library.

spdat<-readOGR(dsn="/Users/ozd504/Google Drive/classes/dem7263/class17/data/Sparks_Sparks_PSP_data_code", layer = "PSP_Data_CLnad83")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/ozd504/Google Drive/classes/dem7263/class17/data/Sparks_Sparks_PSP_data_code", layer: "PSP_Data_CLnad83"
## with 3071 features
## It has 81 fields
#Create a k=5 nearest neighbor set
us.nb5<-knearneigh(coordinates(spdat), k=5)
us.nb5<-knn2nb(us.nb5)
us.wt5<-nb2listw(make.sym.nb(us.nb5), style="W")


nbs<-poly2nb(spdat, queen = T)
nbs
## Neighbour list object:
## Number of regions: 3071 
## Number of nonzero links: 18016 
## Percentage nonzero weights: 0.1910288 
## Average number of links: 5.866493 
## 5 regions with no links:
## 1185 1191 1835 2896 2909
us.lw<-nb2listw(make.sym.nb(nbs), zero.policy = T)

Next, we scale the outcomes and predictors:

spdat$mortz<-scale(spdat$ARSMORT, center=T, scale=T)
spdat$densz<-scale(spdat$PDENSITY, center=T, scale=T)
spdat$rurz<-scale(spdat$PRURAL, center=T, scale=T)
spdat$blacz<-scale(spdat$PBLACK, center=T, scale=T)
spdat$hisz<-scale(spdat$PHISP, center=T, scale=T)
spdat$femz<-scale(spdat$FEMHH, center=T, scale=T)
spdat$unemz<-scale(spdat$UNEMP, center=T, scale=T)
spdat$povz<-scale(spdat$PPERSONPO, center=T, scale=T)
spdat$incz<-scale(spdat$MEDHHINC, center=T, scale=T)
spdat$hvalz<-scale(spdat$MEDHVAL, center=T, scale=T)
spdat$giniz<-scale(spdat$GINI001, center=T, scale=T)

Moran’s I for the mortality rate:

moran.test(spdat$mortz,listw = us.lw , zero.policy = T)
## 
##  Moran I test under randomisation
## 
## data:  spdat$mortz  
## weights: us.lw  
## 
## Moran I statistic standard deviate = 47.396, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.5087738773     -0.0003262643      0.0001153784

Missing data

It’s common to see some polygons with missing data for the outcome, below I show how to deal with this. First I knock out some observations to mimic random missing data.

samp<-sample(1:dim(spdat@data)[1],size = 500, replace = F)
spdat$newmort<-spdat$mortz
spdat$newmort[samp]<-NA

summary(spdat$mortz)
##        V1           
##  Min.   :-6.208043  
##  1st Qu.:-0.671315  
##  Median : 0.008418  
##  Mean   : 0.000000  
##  3rd Qu.: 0.658969  
##  Max.   : 6.743959
summary(spdat$newmort)
##        V1         
##  Min.   :-5.8223  
##  1st Qu.:-0.6733  
##  Median : 0.0076  
##  Mean   :-0.0005  
##  3rd Qu.: 0.6551  
##  Max.   : 6.7440  
##  NA's   :500

This is what you see if you try to analyze this normally:

moran.test(spdat$newmort,listw = us.lw)
## Error in na.fail.default(x): missing values in object

FAIL

So, we can tell R to not use the missing values:

moran.test(spdat$newmort,listw = us.wt5, na.action = na.omit, zero.policy = T)
## 
##  Moran I test under randomisation
## 
## data:  spdat$newmort  
## weights: us.wt5 
## omitted: 1, 6, 9, 15, 24, 37, 39, 50, 61, 69, 79, 84, 88, 92, 93, 103, 106, 113, 118, 119, 129, 143, 144, 156, 160, 161, 169, 184, 190, 205, 224, 237, 240, 258, 259, 274, 283, 287, 288, 289, 292, 299, 300, 304, 313, 325, 330, 333, 336, 345, 357, 362, 363, 368, 371, 375, 400, 408, 413, 415, 420, 421, 426, 428, 434, 443, 449, 458, 460, 467, 471, 474, 483, 487, 492, 509, 510, 522, 527, 545, 559, 561, 568, 572, 575, 577, 582, 583, 593, 598, 608, 610, 620, 633, 646, 649, 653, 672, 677, 682, 685, 688, 700, 701, 714, 720, 726, 740, 741, 751, 752, 754, 765, 778, 782, 786, 795, 800, 804, 807, 815, 818, 843, 853, 855, 856, 860, 861, 863, 867, 879, 892, 897, 903, 912, 915, 919, 920, 926, 929, 930, 940, 946, 948, 957, 972, 977, 981, 988, 992, 998, 999, 1001, 1003, 1004, 1013, 1015, 1018, 1021, 1022, 1024, 1032, 1043, 1045, 1046, 1050, 1055, 1062, 1071, 1073, 1075, 1085, 1093, 1096, 1097, 1100, 1125, 1134, 1140, 1141, 1153, 1162, 1166, 1171, 1179, 1181, 1189, 1200, 1207, 1209, 1215, 1229, 1236, 1250, 1266, 1303, 1305, 1310, 1311, 1314, 1320, 1328, 1336, 1337, 1342, 1347, 1359, 1368, 1379, 1385, 1386, 1390, 1391, 1394, 1405, 1410, 1411, 1440, 1446, 1455, 1460, 1461, 1470, 1476, 1493, 1497, 1499, 1501, 1502, 1516, 1519, 1520, 1523, 1529, 1530, 1531, 1535, 1536, 1545, 1546, 1547, 1548, 1549, 1576, 1582, 1593, 1594, 1599, 1601, 1605, 1608, 1609, 1610, 1613, 1618, 1623, 1627, 1629, 1631, 1649, 1664, 1665, 1674, 1677, 1678, 1681, 1700, 1714, 1717, 1721, 1724, 1726, 1728, 1736, 1738, 1749, 1750, 1762, 1765, 1770, 1776, 1777, 1779, 1780, 1787, 1810, 1821, 1829, 1833, 1837, 1840, 1844, 1845, 1846, 1851, 1857, 1858, 1869, 1883, 1898, 1902, 1906, 1907, 1917, 1922, 1924, 1929, 1930, 1934, 1937, 1941, 1947, 1955, 1963, 1971, 1984, 1991, 1993, 1999, 2011, 2016, 2018, 2022, 2023, 2040, 2048, 2055, 2056, 2059, 2064, 2065, 2066, 2068, 2070, 2072, 2077, 2087, 2091, 2095, 2100, 2118, 2120, 2121, 2125, 2126, 2127, 2129, 2133, 2141, 2158, 2161, 2163, 2171, 2180, 2181, 2186, 2187, 2193, 2194, 2196, 2198, 2204, 2205, 2213, 2221, 2224, 2225, 2227, 2235, 2245, 2251, 2260, 2280, 2282, 2285, 2286, 2296, 2300, 2301, 2304, 2311, 2317, 2319, 2321, 2332, 2343, 2347, 2362, 2368, 2371, 2378, 2389, 2400, 2404, 2415, 2425, 2429, 2437, 2439, 2440, 2442, 2452, 2481, 2486, 2503, 2507, 2509, 2515, 2519, 2535, 2548, 2550, 2580, 2586, 2589, 2597, 2598, 2610, 2614, 2615, 2619, 2620, 2626, 2627, 2628, 2631, 2639, 2652, 2665, 2676, 2680, 2687, 2693, 2696, 2697, 2699, 2700, 2703, 2710, 2715, 2717, 2727, 2736, 2744, 2755, 2761, 2762, 2768, 2769, 2772, 2780, 2781, 2786, 2795, 2803, 2817, 2836, 2859, 2865, 2877, 2881, 2885, 2891, 2893, 2895, 2897, 2901, 2905, 2908, 2912, 2920, 2921, 2928, 2941, 2943, 2952, 2961, 2963, 2964, 2965, 2970, 2971, 2972, 2974, 2986, 3000, 3005, 3008, 3009, 3020, 3041, 3045, 3048, 3050, 3051, 3058, 3063, 3064, 3066, 3069 
## 
## Moran I statistic standard deviate = 39.545, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.5151667220     -0.0003891051      0.0001699708

Which tells R to omit the missing values and allow for there to be empty neighbor sets zero.policy=T. Interestingly, this is also what you do if you have islands in the data which have no neighbors.

Models

Now we show how to fit the various models described above. We start with OLS and move through the various other models.

fit.ols<-lm(mortz~rurz+blacz+hisz+unemz+hvalz+densz+giniz,data=spdat)
summary(fit.ols)
## 
## Call:
## lm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz + densz + 
##     giniz, data = spdat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3588 -0.5376  0.0029  0.5319  6.9022 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.519e-16  1.542e-02   0.000 1.000000    
## rurz        -1.775e-01  1.790e-02  -9.915  < 2e-16 ***
## blacz        1.825e-01  1.762e-02  10.355  < 2e-16 ***
## hisz        -1.819e-01  1.694e-02 -10.742  < 2e-16 ***
## unemz        1.414e-01  1.709e-02   8.271  < 2e-16 ***
## hvalz       -2.749e-01  1.877e-02 -14.648  < 2e-16 ***
## densz        6.496e-02  1.728e-02   3.758 0.000174 ***
## giniz        2.288e-01  1.741e-02  13.141  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8543 on 3063 degrees of freedom
## Multiple R-squared:  0.2719, Adjusted R-squared:  0.2702 
## F-statistic: 163.4 on 7 and 3063 DF,  p-value: < 2.2e-16
lm.morantest(fit.ols, listw=us.wt5)
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz +
## densz + giniz, data = spdat)
## weights: us.wt5
## 
## Moran I statistic standard deviate = 35.747, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     0.3824920327    -0.0014972584     0.0001153862
set.seed(1234)

#SAR - Lag model
fit.lag<-lagsarlm(mortz~rurz+blacz+hisz+unemz+hvalz+densz+giniz, spdat, listw=us.wt5, type="lag", method="spam")
summary(fit.lag, Nagelkerke=T)
## 
## Call:lagsarlm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz + 
##     densz + giniz, data = spdat, listw = us.wt5, type = "lag", 
##     method = "spam")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -5.306204 -0.414023 -0.014568  0.396916  7.192214 
## 
## Type: lag 
## Coefficients: (numerical Hessian approximate standard errors) 
##               Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) -0.0007955  0.0024105  -0.3300    0.7414
## rurz        -0.1555929  0.0149062 -10.4381 < 2.2e-16
## blacz        0.0624896  0.0151577   4.1226 3.746e-05
## hisz        -0.0882377  0.0143535  -6.1475 7.873e-10
## unemz        0.0981214  0.0142580   6.8819 5.908e-12
## hvalz       -0.2077946  0.0156986 -13.2365 < 2.2e-16
## densz        0.0602777  0.0143202   4.2093 2.562e-05
## giniz        0.1378480  0.0147634   9.3371 < 2.2e-16
## 
## Rho: 0.57296, LR test value: 903.02, p-value: < 2.22e-16
## Approximate (numerical Hessian) standard error: 0.016395
##     z-value: 34.948, p-value: < 2.22e-16
## Wald statistic: 1221.4, p-value: < 2.22e-16
## 
## Log likelihood: -3418.366 for lag model
## ML residual variance (sigma squared): 0.50497, (sigma: 0.71061)
## Nagelkerke pseudo-R-squared: 0.45737 
## Number of observations: 3071 
## Number of parameters estimated: 10 
## AIC: 6856.7, (AIC for lm: 7757.8)
#SAR - Error model
fit.err<-errorsarlm(mortz~rurz+blacz+hisz+unemz+hvalz+densz+giniz, spdat, listw=us.wt5, etype="error", method="spam")
summary(fit.err, Nagelkerke=T)
## 
## Call:errorsarlm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz + 
##     densz + giniz, data = spdat, listw = us.wt5, etype = "error", 
##     method = "spam")
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -5.2681368 -0.4012841 -0.0093146  0.3899623  7.0884112 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) -0.0020848  0.0348884  -0.0598 0.9523498
## rurz        -0.1483483  0.0159388  -9.3074 < 2.2e-16
## blacz        0.0890120  0.0251400   3.5406 0.0003992
## hisz        -0.0429855  0.0255515  -1.6823 0.0925088
## unemz        0.1121305  0.0182787   6.1345 8.543e-10
## hvalz       -0.3001374  0.0218386 -13.7435 < 2.2e-16
## densz        0.0892346  0.0181584   4.9142 8.913e-07
## giniz        0.1435818  0.0160992   8.9186 < 2.2e-16
## 
## Lambda: 0.6357, LR test value: 895.93, p-value: < 2.22e-16
## Approximate (numerical Hessian) standard error: 0.017514
##     z-value: 36.298, p-value: < 2.22e-16
## Wald statistic: 1317.5, p-value: < 2.22e-16
## 
## Log likelihood: -3421.91 for error model
## ML residual variance (sigma squared): 0.49608, (sigma: 0.70433)
## Nagelkerke pseudo-R-squared: 0.45612 
## Number of observations: 3071 
## Number of parameters estimated: 10 
## AIC: 6863.8, (AIC for lm: 7757.8)
#Spatial Durbin Model
fit.durb<-lagsarlm(mortz~rurz+blacz+hisz+unemz+hvalz+densz+giniz, spdat, listw=us.wt5, type="mixed", method="spam")
summary(fit.durb, Nagelkerke=T)
## 
## Call:lagsarlm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz + 
##     densz + giniz, data = spdat, listw = us.wt5, type = "mixed", 
##     method = "spam")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -4.447116 -0.401981 -0.023728  0.393768  7.184011 
## 
## Type: mixed 
## Coefficients: (numerical Hessian approximate standard errors) 
##                Estimate  Std. Error  z value  Pr(>|z|)
## (Intercept) -0.00040223  0.00099558  -0.4040 0.6861971
## rurz        -0.15800388  0.01591124  -9.9303 < 2.2e-16
## blacz       -0.02218746  0.03103484  -0.7149 0.4746578
## hisz         0.08280240  0.03175921   2.6072 0.0091288
## unemz        0.09226346  0.01903100   4.8481 1.247e-06
## hvalz       -0.32223768  0.02391956 -13.4717 < 2.2e-16
## densz        0.10314541  0.01882519   5.4791 4.275e-08
## giniz        0.15399858  0.01629208   9.4524 < 2.2e-16
## lag.rurz     0.03243152  0.02771247   1.1703 0.2418859
## lag.blacz    0.08255284  0.03572120   2.3110 0.0208311
## lag.hisz    -0.20589230  0.03525479  -5.8401 5.216e-09
## lag.unemz   -0.05223098  0.02673870  -1.9534 0.0507740
## lag.hvalz    0.22819954  0.03207903   7.1137 1.130e-12
## lag.densz   -0.10780993  0.02779247  -3.8791 0.0001048
## lag.giniz    0.04661633  0.02829491   1.6475 0.0994519
## 
## Rho: 0.5848, LR test value: 805.91, p-value: < 2.22e-16
## Approximate (numerical Hessian) standard error: 0.017519
##     z-value: 33.381, p-value: < 2.22e-16
## Wald statistic: 1114.3, p-value: < 2.22e-16
## 
## Log likelihood: -3364.203 for mixed model
## ML residual variance (sigma squared): 0.48576, (sigma: 0.69696)
## Nagelkerke pseudo-R-squared: 0.47618 
## Number of observations: 3071 
## Number of parameters estimated: 17 
## AIC: 6762.4, (AIC for lm: 7566.3)
#Spatial Durbin Error Model
fit.errdurb<-errorsarlm(mortz~rurz+blacz+hisz+unemz+hvalz+densz+giniz, spdat, listw=us.wt5, etype="emixed", method="spam")
summary(fit.errdurb, Nagelkerke=T)
## 
## Call:errorsarlm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz + 
##     densz + giniz, data = spdat, listw = us.wt5, etype = "emixed", 
##     method = "spam")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -4.539759 -0.404916 -0.017887  0.391254  7.142038 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) -0.0011490  0.0308848  -0.0372 0.9703245
## rurz        -0.1591970  0.0167432  -9.5082 < 2.2e-16
## blacz        0.0080199  0.0277935   0.2886 0.7729240
## hisz         0.0495444  0.0295538   1.6764 0.0936573
## unemz        0.0986970  0.0181457   5.4391 5.354e-08
## hvalz       -0.3113460  0.0226300 -13.7581 < 2.2e-16
## densz        0.0969662  0.0181388   5.3458 9.002e-08
## giniz        0.1695962  0.0166093  10.2109 < 2.2e-16
## lag.rurz    -0.0262226  0.0412557  -0.6356 0.5250290
## lag.blacz    0.1656563  0.0416413   3.9782 6.945e-05
## lag.hisz    -0.2846636  0.0419083  -6.7925 1.102e-11
## lag.unemz    0.0295638  0.0382490   0.7729 0.4395644
## lag.hvalz    0.0737843  0.0440600   1.6746 0.0940063
## lag.densz   -0.0637188  0.0375589  -1.6965 0.0897906
## lag.giniz    0.1496058  0.0409451   3.6538 0.0002584
## 
## Lambda: 0.5917, LR test value: 783.31, p-value: < 2.22e-16
## Approximate (numerical Hessian) standard error: 0.016234
##     z-value: 36.449, p-value: < 2.22e-16
## Wald statistic: 1328.5, p-value: < 2.22e-16
## 
## Log likelihood: -3375.505 for error model
## ML residual variance (sigma squared): 0.48832, (sigma: 0.6988)
## Nagelkerke pseudo-R-squared: 0.47231 
## Number of observations: 3071 
## Number of parameters estimated: 17 
## AIC: 6785, (AIC for lm: 7566.3)
#SAC Model
fit.sac<-sacsarlm(mortz~rurz+blacz+hisz+unemz+hvalz+densz+giniz, spdat, listw=us.wt5, type="sac", method="MC")
summary(fit.sac, Nagelkerke=T)
## 
## Call:sacsarlm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz + 
##     densz + giniz, data = spdat, listw = us.wt5, type = "sac", 
##     method = "MC")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -4.413474 -0.342289 -0.011137  0.331609  6.392057 
## 
## Type: sac 
## Coefficients: (numerical Hessian approximate standard errors) 
##               Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) -0.0039573  0.0365220  -0.1084   0.91371
## rurz        -0.1163918  0.0136058  -8.5546 < 2.2e-16
## blacz        0.0050483  0.0124803   0.4045   0.68584
## hisz         0.0519999  0.0269410   1.9301   0.05359
## unemz        0.0866326  0.0161794   5.3545 8.579e-08
## hvalz       -0.2851716  0.0203905 -13.9855 < 2.2e-16
## densz        0.0900695  0.0161722   5.5694 2.556e-08
## giniz        0.1150812  0.0136554   8.4275 < 2.2e-16
## 
## Rho: -0.74667
## Approximate (numerical Hessian) standard error: 0.038566
##     z-value: -19.361, p-value: < 2.22e-16
## Lambda: 0.90078
## Approximate (numerical Hessian) standard error: 0.0095761
##     z-value: 94.065, p-value: < 2.22e-16
## 
## LR test value: 1092.2, p-value: < 2.22e-16
## 
## Log likelihood: -3323.8 for sac model
## ML residual variance (sigma squared): 0.36588, (sigma: 0.60488)
## Nagelkerke pseudo-R-squared: 0.48978 
## Number of observations: 3071 
## Number of parameters estimated: 11 
## AIC: 6669.6, (AIC for lm: 7757.8)
#SMA Model
fit.sma<-spautolm(mortz~rurz+blacz+hisz+unemz+hvalz+densz+giniz, spdat, listw=us.wt5, family="SMA")
summary(fit.sma)
## 
## Call: spautolm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz + 
##     densz + giniz, data = spdat, listw = us.wt5, family = "SMA")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -5.179461 -0.405557 -0.012908  0.389163  7.089521 
## 
## Coefficients: 
##               Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) -0.0011724  0.0233821  -0.0501      0.96
## rurz        -0.1563420  0.0164214  -9.5206 < 2.2e-16
## blacz        0.1438620  0.0219424   6.5563 5.514e-11
## hisz        -0.1106442  0.0217567  -5.0855 3.666e-07
## unemz        0.1271424  0.0181033   7.0232 2.169e-12
## hvalz       -0.2876579  0.0208977 -13.7651 < 2.2e-16
## densz        0.0793319  0.0182085   4.3569 1.319e-05
## giniz        0.1633269  0.0165091   9.8932 < 2.2e-16
## 
## Lambda: 0.63795 LR test value: 671.92 p-value: < 2.22e-16 
## Numerical Hessian standard error of lambda: 0.025865 
## 
## Log likelihood: -3533.917 
## ML residual variance (sigma squared): 0.62581, (sigma: 0.79108)
## Number of observations: 3071 
## Number of parameters estimated: 10 
## AIC: 7087.8

Using the Lagrange Multiplier Test (LMT)

The so-called Lagrange Multiplier (econometrician’s jargon for a score test) test. These tests compare the model fits from the OLS, spatial error, and spatial lag models using the method of the score test.

For those who don’t remember, the score test is a test based on the relative change in the first derivative of the likelihood function around the maximum likelihood. The particular thing here that is affecting the value of this derivative is the autoregressive parameter, \(\rho\) or \(\lambda\). In the OLS model \(\rho\) or \(\lambda\) = 0 (so both the lag and error models simplify to OLS), but as this parameter changes, so does the likelihood for the model, hence why the derivative of the likelihood function is used. This is all related to how the estimation routines estimate the value of \(\rho\) or \(\lambda\).

In general, you fit the OLS model to your dependent variable, then submit the OLS model fit to the LMT testing procedure.

Then you look to see which model (spatial error, or spatial lag) has the highest value for the test.

Enter the uncertainty… So how much bigger, you might say?

Well, drastically bigger, if the LMT for the error model is 2500 and the LMT for the lag model is 2480, this is NOT A BIG DIFFERENCE, only about 1%. If you see a LMT for the error model of 2500 and a LMT for the lag model of 250, THIS IS A BIG DIFFERENCE.

So what if you don’t see a BIG DIFFERENCE, HOW DO YOU DECIDE WHICH MODEL TO USE???

Well, you could think more, but who has time for that.

The econometricians have thought up a “better” LMT test, the so-called robust LMT, robust to what I’m not sure, but it is said that it can settle such problems of a “not so big difference” between the lag and error model specifications.

So what do you do? In general, think about your problem before you run your analysis, should this fail you, proceed with using the LMT, if this is inconclusive, look at the robust LMT, and choose the model which has the larger value for this test.

Here’s how we do the Lagrange Multiplier test in R:

lm.LMtests(fit.ols, listw=us.wt5, test="all")
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz +
## densz + giniz, data = spdat)
## weights: us.wt5
## 
## LMerr = 1258.7, df = 1, p-value < 2.2e-16
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz +
## densz + giniz, data = spdat)
## weights: us.wt5
## 
## LMlag = 1305.8, df = 1, p-value < 2.2e-16
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz +
## densz + giniz, data = spdat)
## weights: us.wt5
## 
## RLMerr = 46.465, df = 1, p-value = 9.325e-12
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz +
## densz + giniz, data = spdat)
## weights: us.wt5
## 
## RLMlag = 93.651, df = 1, p-value < 2.2e-16
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = mortz ~ rurz + blacz + hisz + unemz + hvalz +
## densz + giniz, data = spdat)
## weights: us.wt5
## 
## SARMA = 1352.3, df = 2, p-value < 2.2e-16

There is a 3.74% difference the regular LM test between the error and lag models, but a 101.72% difference in the Robust LM tests. In this case, I would say that either the lag model looks like the best one, using the Robust Lagrange multiplier test, or possibly the SARMA model, since it’s test is 3.56% difference between it and the lag model. Unfortunately, there is no a robust test for SARMA model.

Of course, the AIC is also your friend:

AICs<-c(AIC(fit.ols),AIC(fit.lag), AIC(fit.err), AIC(fit.durb), AIC(fit.errdurb), AIC(fit.sac), AIC(fit.sma))
plot(AICs, type="l", lwd=1.5, xaxt="n", xlab="")
axis(1, at=1:7,labels=F) #6= number of models
labels<-c("OLS", "Lag","Err", "Durbin","Err Durbin", "SAC", "SMA" )
text(1:7, par("usr")[3]-.25, srt=45, adj=1, labels=labels, xpd=T)
mtext(side=1, text="Model Specification", line=3)
symbols(x= which.min(AICs), y=AICs[which.min(AICs)], circles=1, fg=2,lwd=2,add=T)

knitr::kable(data.frame(Models=labels, AIC=round(AICs, 2)))
Models AIC
OLS 7757.75
Lag 6856.73
Err 6863.82
Durbin 6762.41
Err Durbin 6785.01
SAC 6669.60
SMA 7087.83

Which shows that the Spatial Durbin model best fits the data, although the degree of difference between it an the SAC model is small. A likelihood ratio test could be used:

anova(fit.sac, fit.durb)
##          Model df    AIC  logLik Test L.Ratio    p-value
## fit.sac      1 11 6669.6 -3323.8    1                   
## fit.durb     2 17 6762.4 -3364.2    2  80.808 2.4425e-15

Which indicates that the SAC model fits significantly better than the Durbin model. SAC it is!!

Interpreting effects in spatial lag models

In spatial lag models, interpretation of the regression effects is complicated. Each observation will have a direct effect of its predictors, but each observation will also have in indirect effect of the information of its neighbors, although Spatial Error models do not have this issue. In OLS, the impact/effect of a predictor is straight forward:

\[\frac {\delta y_i} {\delta x_{ik}} = \beta_k\] and

\[\frac {\delta y_i} {\delta x_{jk}} = 0\]

but when a model has a spatial lag of either the outcome or a predictor, this becomes more complicated, indeed:

\[\frac {\delta y_i} {\delta x_{jk}}\]

may not = 0,

or \[\frac {\delta y_i} {\delta x_{jk}} = S_r(W)\] , where \[S_r(W) = (I_n - \rho W)^{-1} \beta_k\]

This implies that a change in the ith region’s predictor can affect the jth region’s outcome

  • We have 2 situations:

  • \(S_r(W)_{ii}\), or the direct impact of an observation’s predictor on its own outcome, and:

  • \(S_r(W)_{ij}\), or the indirect impact of an observation’s neighbor’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 neighbors 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 Lesage and Pace, 2009. Here, we use Monte Carlo simulation to obtain simulated distributions of the various impacts. We are looking for the first part of the output and

W <- as(us.wt5, "CsparseMatrix")
trMC <- trW(W, type="MC")
im<-impacts(fit.sac, tr=trMC, R=100)
sums<-summary(im,  zstats=T)
data.frame(sums$res)
##        direct     indirect       total
## 1 -0.12736535  0.060739306 -0.06662605
## 2  0.00552429 -0.002634481  0.00288981
## 3  0.05690249 -0.027136249  0.02976624
## 4  0.09480045 -0.045209420  0.04959103
## 5 -0.31205798  0.148817434 -0.16324054
## 6  0.09856136 -0.047002961  0.05155840
## 7  0.12593117 -0.060055357  0.06587581
data.frame(sums$pzmat)
##             Direct     Indirect        Total
## rurz  0.000000e+00 0.000000e+00 0.000000e+00
## blacz 5.665000e-01 5.686653e-01 5.650648e-01
## hisz  8.833768e-02 8.984462e-02 8.839540e-02
## unemz 1.113372e-08 2.569246e-08 2.065345e-08
## hvalz 0.000000e+00 0.000000e+00 0.000000e+00
## densz 1.658948e-08 2.375242e-08 4.392931e-08
## giniz 6.217249e-15 3.892442e-13 9.325873e-15

We see all variables have a significant direct effect, we also see that the rural, unemployment, home value, density and Gini index all have significant direct impacts. These also have significant indirect impacts

We can likewise see the effects by order of neighbors, similar to what Yang et al(2015) do in their Table 4.

Here, I do this up to 5th order neighbors.

im2<-impacts(fit.sac, tr=trMC, R=100, Q=5)
sums2<-summary(im2,  zstats=T, reportQ=T, short=T)
sums2
## Impact measures (sac, trace):
##            Direct     Indirect       Total
## rurz  -0.12736535  0.060739306 -0.06662605
## blacz  0.00552429 -0.002634481  0.00288981
## hisz   0.05690249 -0.027136249  0.02976624
## unemz  0.09480045 -0.045209420  0.04959103
## hvalz -0.31205798  0.148817434 -0.16324054
## densz  0.09856136 -0.047002961  0.05155840
## giniz  0.12593117 -0.060055357  0.06587581
## =================================
## Impact components
## $direct
##            rurz         blacz         hisz        unemz        hvalz
## Q1 -0.116391778  0.0050483270  0.051999874  0.086632612 -0.285171611
## Q2  0.000000000  0.0000000000  0.000000000  0.000000000  0.000000000
## Q3 -0.011493277  0.0004985045  0.005134804  0.008554665 -0.028159690
## Q4  0.002660082 -0.0001153773 -0.001188434 -0.001979950  0.006517469
## Q5 -0.002796937  0.0001213131  0.001249576  0.002081813 -0.006852777
##           densz        giniz
## Q1  0.090069488  0.115081158
## Q2  0.000000000  0.000000000
## Q3  0.008894044  0.011363858
## Q4 -0.002058498 -0.002630128
## Q5  0.002164402  0.002765442
## 
## $indirect
##           rurz        blacz        hisz       unemz       hvalz
## Q1  0.00000000  0.000000000  0.00000000  0.00000000  0.00000000
## Q2  0.08690603 -0.003769425 -0.03882665 -0.06468581  0.21292855
## Q3 -0.05339668  0.002316005  0.02385582  0.03974417 -0.13082726
## Q4  0.04579118 -0.001986127 -0.02045794 -0.03408325  0.11219302
## Q5 -0.03338008  0.001447813  0.01491308  0.02484543 -0.08178456
##          densz       giniz
## Q1  0.00000000  0.00000000
## Q2 -0.06725201 -0.08592743
## Q3  0.04132089  0.05279541
## Q4 -0.03543539 -0.04527555
## Q5  0.02583109  0.03300420
## 
## $total
##           rurz        blacz        hisz       unemz       hvalz
## Q1 -0.11639178  0.005048327  0.05199987  0.08663261 -0.28517161
## Q2  0.08690603 -0.003769425 -0.03882665 -0.06468581  0.21292855
## Q3 -0.06488996  0.002814509  0.02899062  0.04829883 -0.15898695
## Q4  0.04845126 -0.002101504 -0.02164637 -0.03606320  0.11871049
## Q5 -0.03617701  0.001569126  0.01616265  0.02692724 -0.08863734
##          densz       giniz
## Q1  0.09006949  0.11508116
## Q2 -0.06725201 -0.08592743
## Q3  0.05021493  0.06415927
## Q4 -0.03749389 -0.04790568
## Q5  0.02799549  0.03576964
## 
## ========================================================
## Simulation results (numerical Hessian approximation variance matrix):
## ========================================================
## Simulated z-values:
##            Direct   Indirect       Total
## rurz   -9.3585465  9.8413224  -8.3476191
## blacz   0.3566302 -0.3547153   0.3580224
## hisz    1.6912774 -1.6832481   1.6942661
## unemz   5.9028924 -5.8099935   5.7955877
## hvalz -14.9409817 13.5302269 -13.5308850
## densz   5.7329163 -5.5828904   5.6835178
## giniz   8.1669200 -7.8597172   7.9432420
## 
## Simulated p-values:
##       Direct     Indirect   Total     
## rurz  < 2.22e-16 < 2.22e-16 < 2.22e-16
## blacz 0.721369   0.722803   0.720327  
## hisz  0.090784   0.092327   0.090215  
## unemz 3.5718e-09 6.2475e-09 6.8082e-09
## hvalz < 2.22e-16 < 2.22e-16 < 2.22e-16
## densz 9.8718e-09 2.3655e-08 1.3195e-08
## giniz 2.2204e-16 3.7748e-15 1.9984e-15
## ========================================================
## Simulated impact components z-values:
## $Direct
##         rurz      blacz      hisz     unemz      hvalz     densz     giniz
## Q1 -9.098466  0.3570218  1.692508  5.889645 -14.768410  5.738504  8.152525
## Q2       NaN        NaN       NaN       NaN        NaN       NaN       NaN
## Q3 -8.633573  0.3510692  1.660002  5.276482  -9.215440  4.999099  6.584695
## Q4  6.928893 -0.3477198 -1.631183 -4.654974   6.755963 -4.382151 -5.436598
## Q5 -5.557704  0.3442169  1.595905  4.060864  -5.257384  3.815565  4.526459
## 
## $Indirect
##         rurz      blacz      hisz     unemz     hvalz     densz     giniz
## Q1       NaN        NaN       NaN       NaN       NaN       NaN       NaN
## Q2  9.811045 -0.3541957 -1.680864 -5.766268 12.987782 -5.530625 -7.736218
## Q3 -8.633573  0.3510692  1.660002  5.276482 -9.215440  4.999099  6.584695
## Q4  6.928893 -0.3477198 -1.631183 -4.654974  6.755963 -4.382151 -5.436598
## Q5 -5.557704  0.3442169  1.595905  4.060864 -5.257384  3.815565  4.526459
## 
## $Total
##         rurz      blacz      hisz     unemz      hvalz     densz     giniz
## Q1 -9.098466  0.3570218  1.692508  5.889645 -14.768410  5.738504  8.152525
## Q2  9.811045 -0.3541957 -1.680864 -5.766268  12.987782 -5.530625 -7.736218
## Q3 -8.633573  0.3510692  1.660002  5.276482  -9.215440  4.999099  6.584695
## Q4  6.928893 -0.3477198 -1.631183 -4.654974   6.755963 -4.382151 -5.436598
## Q5 -5.557704  0.3442169  1.595905  4.060864  -5.257384  3.815565  4.526459
## 
## 
## Simulated impact components p-values:
## $Direct
##    rurz       blacz   hisz     unemz      hvalz      densz      giniz     
## Q1 < 2.22e-16 0.72108 0.090549 3.8703e-09 < 2.22e-16 9.5517e-09 4.4409e-16
## Q2 NA         NA      NA       NA         NA         NA         NA        
## Q3 < 2.22e-16 0.72554 0.096914 1.3169e-07 < 2.22e-16 5.7599e-07 4.5582e-11
## Q4 4.2415e-12 0.72805 0.102852 3.2402e-06 1.4189e-11 1.1751e-05 5.4308e-08
## Q5 2.7335e-08 0.73068 0.110510 4.8891e-05 1.4612e-07 0.00013587 5.9980e-06
## 
## $Indirect
##    rurz       blacz   hisz     unemz      hvalz      densz      giniz     
## Q1 NA         NA      NA       NA         NA         NA         NA        
## Q2 < 2.22e-16 0.72319 0.092789 8.1046e-09 < 2.22e-16 3.1909e-08 1.0214e-14
## Q3 < 2.22e-16 0.72554 0.096914 1.3169e-07 < 2.22e-16 5.7599e-07 4.5582e-11
## Q4 4.2415e-12 0.72805 0.102852 3.2402e-06 1.4189e-11 1.1751e-05 5.4308e-08
## Q5 2.7335e-08 0.73068 0.110510 4.8891e-05 1.4612e-07 0.00013587 5.9980e-06
## 
## $Total
##    rurz       blacz   hisz     unemz      hvalz      densz      giniz     
## Q1 < 2.22e-16 0.72108 0.090549 3.8703e-09 < 2.22e-16 9.5517e-09 4.4409e-16
## Q2 < 2.22e-16 0.72319 0.092789 8.1046e-09 < 2.22e-16 3.1909e-08 1.0214e-14
## Q3 < 2.22e-16 0.72554 0.096914 1.3169e-07 < 2.22e-16 5.7599e-07 4.5582e-11
## Q4 4.2415e-12 0.72805 0.102852 3.2402e-06 1.4189e-11 1.1751e-05 5.4308e-08
## Q5 2.7335e-08 0.73068 0.110510 4.8891e-05 1.4612e-07 0.00013587 5.9980e-06

So we see that, for instance, for the direct impact of %rural, -.117/-.127= 91.8% of the effect is due to a county’s own influence on itself, while (0 + -.0113 +0.002698620+-0.002815391)/-.127 = 8.99 % of the effect of %rural comes from other neighboring counties. This is confusing, but shows that the effects of a variable can come from two sources in the SAR models.