Tutorial: Spatial econometrics with R 4

Load package

require("spdep") 
require("spatialreg")
require("raster") 
require("cartography")
require("spData") 
require("rgdal")
require("GeoXp")
require("geodata")
require("sp")
require("sf")

library(rgdal) 
columbus <- readOGR(system.file("shapes/columbus.shp", package = "spData")[1])

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\TRINHHUONG\Documents\R\win-library\4.0\spData\shapes\columbus.shp", layer: "columbus"
## with 49 features
## It has 20 fields
## Integer64 fields read as strings:  COLUMBUS_ COLUMBUS_I POLYID

Wcontig.nb=poly2nb(columbus,queen=FALSE) 
Wcontig.listw=nb2listw(Wcontig.nb) 


olm = lm(CRIME ~ INC + HOVAL,data = columbus) #OLS,HOVAL housing value,INC #household income (in \$1,000), CRIME residential burglaries and vehicle thefts #per thousand households 
INC_lag=lag.listw(Wcontig.listw,columbus@data$INC) #lagInc, lag_hoval
HOVAL_lag= lag.listw(Wcontig.listw,columbus@data$HOVAL)
xlm=lm(CRIME ~ INC + HOVAL + INC_lag + HOVAL_lag,data = columbus) 
xlm

## 
## Call:
## lm(formula = CRIME ~ INC + HOVAL + INC_lag + HOVAL_lag, data = columbus)
## 
## Coefficients:
## (Intercept)          INC        HOVAL      INC_lag    HOVAL_lag  
##     73.9619      -1.1418      -0.2918      -1.2353       0.1828

lm.morantest(xlm,Wcontig.listw)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = CRIME ~ INC + HOVAL + INC_lag + HOVAL_lag, data =
## columbus)
## weights: Wcontig.listw
## 
## Moran I statistic standard deviate = 2.9426, p-value = 0.001627
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.245198298     -0.042105489      0.009532957

################### March 8 #################################
formula=CRIME ~ INC 
xlm2 = spatialreg::lmSLX(CRIME ~ INC + HOVAL,listw = Wcontig.listw, data = columbus, Durbin = formula) #SLM Spatial Lag of X

summary(xlm2)

## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.752  -7.253   0.160   8.595  27.002 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 76.47908    5.93353  12.889  < 2e-16 ***
## INC         -1.21129    0.37228  -3.254  0.00216 ** 
## HOVAL       -0.26275    0.09981  -2.633  0.01157 *  
## lag.INC     -0.93565    0.45003  -2.079  0.04334 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.04 on 45 degrees of freedom
## Multiple R-squared:  0.5916, Adjusted R-squared:  0.5644 
## F-statistic: 21.73 on 3 and 45 DF,  p-value: 7.527e-09

xlm3=lmSLX(CRIME ~ INC + HOVAL,listw=Wcontig.listw,data = columbus) #lagX
isSymmetric(nb2mat(Wcontig.nb,style="B"),check.attributes = FALSE) 

## [1] TRUE

isSymmetric(nb2mat(Wcontig.nb,style="W"),check.attributes = FALSE)

## [1] FALSE

coord <- coordinates(columbus)[, 1:2] 

Wnear.knn=knearneigh(coord, k=5)
Wnear.mat=nb2mat(knn2nb(Wnear.knn),style= "B")
eig=eigen(nb2mat(Wcontig.nb))$values
xlm3

## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Coefficients:
## (Intercept)          INC        HOVAL      lag.INC    lag.HOVAL  
##     73.9619      -1.1418      -0.2918      -1.2353       0.1828

laglm=spatialreg::lagsarlm(formula = CRIME ~ INC + HOVAL, data = columbus, listw = Wcontig.listw) #SAR spatial autoregressive model with WY
summary(laglm)

## 
## Call:spatialreg::lagsarlm(formula = CRIME ~ INC + HOVAL, data = columbus, 
##     listw = Wcontig.listw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -36.92792  -5.56906  -0.92424   6.21002  24.68336 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 45.264975   7.175796  6.3080 2.827e-10
## INC         -1.036346   0.305252 -3.3950 0.0006862
## HOVAL       -0.259418   0.088797 -2.9215 0.0034837
## 
## Rho: 0.42281, LR test value: 9.7192, p-value: 0.0018235
## Asymptotic standard error: 0.11558
##     z-value: 3.6582, p-value: 0.00025398
## Wald statistic: 13.383, p-value: 0.00025398
## 
## Log likelihood: -182.5176 for lag model
## ML residual variance (sigma squared): 95.723, (sigma: 9.7838)
## Number of observations: 49 
## Number of parameters estimated: 5 
## AIC: 375.04, (AIC for lm: 382.75)
## LM test for residual autocorrelation
## test value: 0.50522, p-value: 0.47722

##understanding what lagsarlm is doing 
A_rho <- function(rho, my_listw) 
  { 
  W <- listw2mat(my_listw) 
  n <- nrow(W) 
  diag(n) - rho * listw2mat(my_listw)
} 

hat_beta_rho <- function(rho, x, y, my_listw) 
{ 
  solve(crossprod(x), crossprod(x, A_rho(rho, my_listw) %*% y))
  } 

hat_sigma_rho <- function(rho, x, y, my_listw) 
{ 
  n <- length(y) 
  my_A <- A_rho(rho, my_listw) 
  my_beta <- hat_beta_rho(rho, x, y, my_listw) 
  my_A_inv <- solve(my_A, x %*% my_beta) 
  1 / n * t(y - my_A_inv) %*% 
    crossprod(my_A) %*% (y - my_A_inv) 
} 

my_like <- function(rho, x, y, my_listw) 
{ 
  n <- length(y) 
  my_A <- A_rho(rho, my_listw) 
  determinant(my_A, logarithm = TRUE)$modulus 
  n / 2 * log(hat_sigma_rho(rho, x, y, my_listw)) 
}

columbus_sf <-  read_sf(system.file("shapes/columbus.shp",
                                    package = "spData")[1]) 
y <- st_drop_geometry(columbus_sf)$CRIME

x <- cbind(1, st_drop_geometry(columbus_sf)$HOVAL, st_drop_geometry(columbus_sf)$INC) 
rho_possible <- seq(-0.99, 0.99, 0.01) 
ll_rho <- sapply(rho_possible, my_like, x, y, Wcontig.listw) 
plot(rho_possible, ll_rho, type = "l", xlab = "valeurs de rho", ylab = "log-vraisemblance") 

optimize(my_like, c(-0.99, 0.99), x, y, Wcontig.listw, maximum = T)

## $maximum
## [1] -0.9899191
## 
## $objective
##          [,1]
## [1,] 140.8946

logdet1=log(prod(1 - 0.5 * spatialreg::eigenw(Wcontig.listw))) 
logdet1

## [1] -1.795856

logdet2=determinant(A_rho(0.5, Wcontig.listw), logarithm = TRUE)$modulus

#####spatial durbin 
durbin <- spatialreg::lagsarlm(CRIME ~ INC + HOVAL, data = columbus, listw = Wcontig.listw, type = "mixed") #SARAR: mixed WY and WX
#laglm=spatialreg::lagsarlm(formula = CRIME ~ INC + HOVAL, data = columbus, #listw = Wcontig.listw) #SAR spatial autoregressive model with WY
summary(durbin)

## 
## Call:spatialreg::lagsarlm(formula = CRIME ~ INC + HOVAL, data = columbus, 
##     listw = Wcontig.listw, type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -37.31038  -6.07084  -0.51813   6.58699  23.49170 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 41.175254  12.086052  3.4068 0.0006572
## INC         -0.926376   0.333824 -2.7750 0.0055194
## HOVAL       -0.296256   0.091779 -3.2279 0.0012469
## lag.INC     -0.385647   0.553973 -0.6961 0.4863367
## lag.HOVAL    0.235127   0.186005  1.2641 0.2061962
## 
## Rho: 0.43826, LR test value: 6.0071, p-value: 0.014249
## Asymptotic standard error: 0.14885
##     z-value: 2.9444, p-value: 0.0032362
## Wald statistic: 8.6693, p-value: 0.0032362
## 
## Log likelihood: -181.7106 for mixed model
## ML residual variance (sigma squared): 92.238, (sigma: 9.6041)
## Number of observations: 49 
## Number of parameters estimated: 7 
## AIC: 377.42, (AIC for lm: 381.43)
## LM test for residual autocorrelation
## test value: 0.77351, p-value: 0.37913

###partial durbin 
formula=CRIME ~ INC 
durbin2=spatialreg::lagsarlm(CRIME ~ INC + HOVAL, data = columbus, listw = Wcontig.listw, Durbin = formula)
#durbin2

summary(durbin2)

## 
## Call:spatialreg::lagsarlm(formula = CRIME ~ INC + HOVAL, data = columbus, 
##     listw = Wcontig.listw, Durbin = formula)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -37.00348  -5.71198  -0.86229   6.20462  24.69181 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 46.413766  11.584119  4.0067 6.158e-05
## INC         -1.027956   0.330425 -3.1110  0.001864
## HOVAL       -0.259126   0.088908 -2.9145  0.003562
## lag.INC     -0.058573   0.510691 -0.1147  0.908688
## 
## Rho: 0.41092, LR test value: 5.2372, p-value: 0.022109
## Asymptotic standard error: 0.14845
##     z-value: 2.7681, p-value: 0.0056379
## Wald statistic: 7.6625, p-value: 0.0056379
## 
## Log likelihood: -182.5115 for mixed model
## ML residual variance (sigma squared): 95.993, (sigma: 9.7976)
## Number of observations: 49 
## Number of parameters estimated: 6 
## AIC: 377.02, (AIC for lm: 380.26)
## LM test for residual autocorrelation
## test value: 2.552, p-value: 0.11015

#####Bayesian estimation
Bayesmod=spatialreg::spBreg_lag(CRIME ~ INC + HOVAL, data=columbus, listw=Wcontig.listw)
summary(Bayesmod)

## 
## Iterations = 501:2500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 2000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                 Mean       SD Naive SE Time-series SE
## (Intercept)  46.0711  8.35952 0.186925       0.186925
## INC          -1.0492  0.35168 0.007864       0.007864
## HOVAL        -0.2633  0.09394 0.002101       0.002101
## rho           0.4086  0.12756 0.002852       0.002852
## sige        108.7965 24.17370 0.540540       0.577264
## 
## 2. Quantiles for each variable:
## 
##                2.5%     25%      50%      75%    97.5%
## (Intercept) 29.3672 40.5417  45.8271  51.7303  62.6306
## INC         -1.7265 -1.2899  -1.0560  -0.8124  -0.3505
## HOVAL       -0.4525 -0.3254  -0.2616  -0.1987  -0.0880
## rho          0.1606  0.3207   0.4137   0.4967   0.6509
## sige        71.2679 90.9087 105.9367 122.4151 163.7571

####stsls estimation
s2sls_est=spatialreg::stsls(CRIME ~ INC + HOVAL, data=columbus, listw=Wcontig.listw)
summary(s2sls_est)

## 
## Call:spatialreg::stsls(formula = CRIME ~ INC + HOVAL, data = columbus, 
##     listw = Wcontig.listw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -36.90706  -5.49372  -0.92969   6.21749  24.71632 
## 
## Coefficients: 
##              Estimate Std. Error t value  Pr(>|t|)
## Rho          0.419294   0.187585  2.2352  0.025403
## (Intercept) 45.459092  11.191512  4.0619 4.867e-05
## INC         -1.041009   0.388612 -2.6788  0.007389
## HOVAL       -0.259538   0.092406 -2.8087  0.004975
## 
## Residual variance (sigma squared): 104.33, (sigma: 10.214)

######################## March 9 ###############################
####impacts in Durbin model 
imp_durbin=spatialreg::impacts(durbin, listw=Wcontig.listw,R=200) 
summary(imp_durbin, zstats=TRUE)

## Impact measures (mixed, exact):
##           Direct   Indirect      Total
## INC   -1.0354954 -1.3001390 -2.3356344
## HOVAL -0.2817423  0.1729232 -0.1088192
## ========================================================
## Simulation results ( variance matrix):
## Direct:
## 
## Iterations = 1:200
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 200 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##          Mean      SD Naive SE Time-series SE
## INC   -1.0756 0.30662  0.02168       0.021682
## HOVAL -0.2699 0.09532  0.00674       0.004968
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%     50%     75%    97.5%
## INC   -1.6337 -1.3181 -1.0805 -0.8378 -0.51697
## HOVAL -0.4391 -0.3299 -0.2759 -0.2060 -0.06951
## 
## ========================================================
## Indirect:
## 
## Iterations = 1:200
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 200 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##          Mean     SD Naive SE Time-series SE
## INC   -1.3041 0.8344   0.0590         0.0590
## HOVAL  0.1761 0.3592   0.0254         0.0254
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%     50%     75%  97.5%
## INC   -2.8866 -1.8091 -1.2149 -0.8255 0.2982
## HOVAL -0.3767 -0.0163  0.1617  0.3370 0.8241
## 
## ========================================================
## Total:
## 
## Iterations = 1:200
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 200 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean     SD Naive SE Time-series SE
## INC   -2.37968 0.8915  0.06304        0.06304
## HOVAL -0.09387 0.3902  0.02759        0.02759
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%     50%      75%   97.5%
## INC   -4.4077 -2.8199 -2.3286 -1.85931 -0.7254
## HOVAL -0.6709 -0.3461 -0.1164  0.08175  0.6279
## 
## ========================================================
## Simulated standard errors
##           Direct  Indirect     Total
## INC   0.30662495 0.8344063 0.8914519
## HOVAL 0.09531741 0.3591790 0.3902294
## 
## Simulated z-values:
##          Direct   Indirect     Total
## INC   -3.507864 -1.5628887 -2.669446
## HOVAL -2.832100  0.4902235 -0.240552
## 
## Simulated p-values:
##       Direct     Indirect Total    
## INC   0.00045172 0.11808  0.0075976
## HOVAL 0.00462433 0.62398  0.8099024

###test of rho = zero 
lm.LMtests(olm,Wcontig.listw,test=c("LMlag","RLMlag"))

## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = CRIME ~ INC + HOVAL, data = columbus)
## weights: Wcontig.listw
## 
## LMlag = 8.7599, df = 1, p-value = 0.003079
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = CRIME ~ INC + HOVAL, data = columbus)
## weights: Wcontig.listw
## 
## RLMlag = 3.0722, df = 1, p-value = 0.07964

###prediction in LAG model 
###prediction in LAG model

y_TC <- predict(laglm, listw = Wcontig.listw, pred.type = "TC", all.data = TRUE, zero.policy = FALSE)

y_TS <- predict(laglm, listw = Wcontig.listw, pred.type = "TS", all.data = TRUE, zero.policy = FALSE)

y_BP <- predict(laglm, listw = Wcontig.listw, pred.type = "BP", all.data = TRUE, zero.policy = FALSE)

pred<-list("TC"=y_TC,"TS"=y_TS,"BP"=y_BP)

pred <- data.frame(TC = as.numeric(y_TC), TS = as.numeric(y_TS), BP = as.numeric(y_BP))
sapply(pred, function(x) mean((columbus_sf$CRIME - x) ^ 2))

##        TC        TS        BP 
## 117.22186  95.72350  89.06412

###fitting the SEM model 
sem=spautolm(CRIME ~ INC + HOVAL, data = columbus, listw = Wcontig.listw) # drop spatialreg::
summary(sem)

## 
## Call: 
## spautolm(formula = CRIME ~ INC + HOVAL, data = columbus, listw = Wcontig.listw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -34.4004  -6.0196  -1.3928   7.0105  24.0611 
## 
## Coefficients: 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 60.375187   5.325070 11.3379 < 2.2e-16
## INC         -0.961044   0.331146 -2.9022  0.003706
## HOVAL       -0.303198   0.092641 -3.2728  0.001065
## 
## Lambda: 0.54847 LR test value: 8.1273 p-value: 0.0043603 
## Numerical Hessian standard error of lambda: 0.15002 
## 
## Log likelihood: -183.3136 
## ML residual variance (sigma squared): 94.968, (sigma: 9.7451)
## Number of observations: 49 
## Number of parameters estimated: 5 
## AIC: 376.63

sem2=errorsarlm(CRIME ~ INC + HOVAL, data = columbus, listw = Wcontig.listw)
summary(sem2)

## 
## Call:
## errorsarlm(formula = CRIME ~ INC + HOVAL, data = columbus, listw = Wcontig.listw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -34.4004  -6.0196  -1.3928   7.0105  24.0611 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 60.375188   5.325070 11.3379 < 2.2e-16
## INC         -0.961044   0.331146 -2.9022  0.003706
## HOVAL       -0.303198   0.092641 -3.2728  0.001065
## 
## Lambda: 0.54847, LR test value: 8.1273, p-value: 0.0043603
## Asymptotic standard error: 0.13138
##     z-value: 4.1747, p-value: 2.9832e-05
## Wald statistic: 17.428, p-value: 2.9832e-05
## 
## Log likelihood: -183.3136 for error model
## ML residual variance (sigma squared): 94.968, (sigma: 9.7451)
## Number of observations: 49 
## Number of parameters estimated: 5 
## AIC: 376.63, (AIC for lm: 382.75)

####testing strategy 
lm.LMtests(olm,Wcontig.listw,test=c("LMerr"))

## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = CRIME ~ INC + HOVAL, data = columbus)
## weights: Wcontig.listw
## 
## LMErr = 5.8149, df = 1, p-value = 0.01589

lm.LMtests(olm,Wcontig.listw,test=c("RLMerr", "RLMerr"))

## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = CRIME ~ INC + HOVAL, data = columbus)
## weights: Wcontig.listw
## 
## RLMerr = 0.12715, df = 1, p-value = 0.7214
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = CRIME ~ INC + HOVAL, data = columbus)
## weights: Wcontig.listw
## 
## RLMerr = 0.12715, df = 1, p-value = 0.7214