Spatial Regression

1 Libraries and raw data

library(tidyverse)
library(magrittr)
library(sf)
library(tmap)

library(spdep) # Spatial weights, Moran's I, and LISA
library(spatialreg) # Spatial Regression

Loading required package: Matrix

Attaching package: ‘Matrix’

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack


Attaching package: ‘spatialreg’

The following objects are masked from ‘package:spdep’:

    get.ClusterOption, get.coresOption, get.mcOption, get.VerboseOption,
    get.ZeroPolicyOption, set.ClusterOption, set.coresOption, set.mcOption,
    set.VerboseOption, set.ZeroPolicyOption

load("../04_MexicoCity_HTS/Mexico_HTS.RData")
ls()

[1] "MexicoCity_HTS"

1.1 Remove districts

summary(MexicoCity_HTS$TotalPeople)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    634     968    1035    1027    1103    1333       1 

MexicoCity_HTS %>% filter(is.na(TotalPeople)) %>% select(Distrito)

Simple feature collection with 1 feature and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -99.09278 ymin: 19.41697 xmax: -99.04731 ymax: 19.45018
Geodetic CRS:  WGS 84
  Distrito                       geometry
1      034 MULTIPOLYGON (((-99.08461 1...

MexicoCity_HTS %<>% filter(Distrito != "034")

summary(MexicoCity_HTS$Total_Homes)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  272.0   285.0   293.0   293.7   301.0   352.0 

1.2 Map preparation

MexicoCity_Fringe <- st_read("../04_MexicoCity_HTS/UrbanFringe_MexicoCity.shp")

Reading layer `UrbanFringe_MexicoCity' from data source 
  `C:\Users\Orlan\Dropbox\Teaching\SpatialAnalysis\Tutorials\04_MexicoCity_HTS\UrbanFringe_MexicoCity.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 5 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -99.36492 ymin: 19.04824 xmax: -98.9403 ymax: 19.59276
Geodetic CRS:  WGS 84

BB_Districts <- st_bbox(MexicoCity_HTS)
BB_Districts[1] <- -99.8
BB_Districts

     xmin      ymin      xmax      ymax 
-99.80000  18.93534 -98.59687  20.06826 

2 OLS

2.1 Weighted Matrix

Neigh_ <- poly2nb(MexicoCity_HTS) 
WM_ <- nb2listw(Neigh_, style = "W")

Let’s keep the following piece of code in mind:

WM2 <- as(WM_, "CsparseMatrix")
trMC <- trW(WM2, type="MC")

2.2 Specification?

MexicoCity_HTS %<>%
  mutate(CarTrips_PerPopulation = Trips_Automovil/TotalPeople,
         Prop_Male = Total_Male/TotalPeople,
         Prop_HomesWithCars = Homes_With_Cars/Total_Homes,
         Prop_HomesWithBicycles = Homes_With_Bicycles/Total_Homes,
         Prop_Education_Low = TotalEducation_Low/TotalPeople,
         Prop_Education_Medium = TotalEducation_Medium/TotalPeople,
         Prop_Education_High = TotalEducation_High/TotalPeople)

ModelEquation <- "CarTrips_PerPopulation ~ 
                          Prop_Male + 
                          Prop_HomesWithCars + 
                          Prop_HomesWithBicycles +
                          # Prop_Education_Low + 
                          Age"

Model_OLS <- lm(ModelEquation, data = MexicoCity_HTS) 

2.3 Model Fitting

summary(Model_OLS)

Call:
lm(formula = ModelEquation, data = MexicoCity_HTS)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.42660 -0.13935 -0.04958  0.08185  2.39132 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)            -1.926931   0.922740  -2.088   0.0381 *  
Prop_Male              -0.116088   1.619548  -0.072   0.9429    
Prop_HomesWithCars      1.874276   0.207355   9.039  < 2e-16 ***
Prop_HomesWithBicycles -0.021345   0.173339  -0.123   0.9021    
Age                     0.051294   0.008715   5.886 1.79e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2863 on 188 degrees of freedom
Multiple R-squared:  0.6505,    Adjusted R-squared:  0.643 
F-statistic: 87.46 on 4 and 188 DF,  p-value: < 2.2e-16

2.4 VIF

regclass::VIF(Model_OLS)

             Prop_Male     Prop_HomesWithCars Prop_HomesWithBicycles 
              1.260071               1.751949               1.094062 
                   Age 
              2.076368 

2.5 Residuals - Statistical tests

2.5.1 Correlation

H0 (null hypothesis): There is no correlation among the residuals. HA (alternative hypothesis): The residuals are autocorrelated.

p-value is less than 0.05, reject the null hypothesis and conclude that the residuals in this regression model are autocorrelated

car::durbinWatsonTest(Model_OLS)

 lag Autocorrelation D-W Statistic p-value
   1       0.2576538      1.379249       0
 Alternative hypothesis: rho != 0

2.5.2 Normality

The null hypothesis is that the data are normally distributed If the p-value is less than the significance level (usually 0.05), we reject the null hypothesis and conclude that the data are not normally distributed.

tseries::jarque.bera.test(residuals(Model_OLS))

    Jarque Bera Test

data:  residuals(Model_OLS)
X-squared = 6389.5, df = 2, p-value < 2.2e-16

If the p-value is less than the significance level (usually 0.05), we reject the null hypothesis and conclude that the residuals do not follow a normal distribution.

ks.test(residuals(Model_OLS), "pnorm")

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  residuals(Model_OLS)
D = 0.34029, p-value < 2.2e-16
alternative hypothesis: two-sided

2.5.3 Heteroscedasticity

If the p-value is less than the significance level (usually 0.05), we reject the null hypothesis and conclude that there is evidence of heteroscedasticity in the residuals.

lmtest::bptest(Model_OLS)

    studentized Breusch-Pagan test

data:  Model_OLS
BP = 10.586, df = 4, p-value = 0.03163

3 OLS - Spatial autocorrelation of the residuals

lm.morantest(Model_OLS, WM_)

    Global Moran I for regression residuals

data:  
model: lm(formula = ModelEquation, data = MexicoCity_HTS)
weights: WM_

Moran I statistic standard deviate = 4.8456, p-value = 6.311e-07
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
     0.198201508     -0.014234611      0.001922028 

MexicoCity_HTS %<>% mutate(Residuals_OLS = residuals(Model_OLS)) 

tm_shape(MexicoCity_HTS, bbox = BB_Districts) +
  tm_fill("Residuals_OLS", style = "quantile", title = "Quantiles") +
  tm_borders(alpha = 0.6) +
  tm_shape(MexicoCity_Fringe) + 
  tm_polygons(alpha = 0, border.col = "black", lwd = 1.4, lty = "solid") +
  tm_compass(size = 2, type = "arrow", position = c(0.85,0.87)) + 
  tm_scale_bar(position = c(0.4,0.02)) +
  tm_layout(title = "Residuals - OLS", legend.position = c(0.01,0.15), scale = 1)

moran.test(MexicoCity_HTS$Residuals_OLS, WM_)

    Moran I test under randomisation

data:  MexicoCity_HTS$Residuals_OLS  
weights: WM_    

Moran I statistic standard deviate = 4.9297, p-value = 4.117e-07
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.198201508      -0.005208333       0.001702542 

moran.mc(MexicoCity_HTS$Residuals_OLS, WM_, nsim = 10000)

    Monte-Carlo simulation of Moran I

data:  MexicoCity_HTS$Residuals_OLS 
weights: WM_  
number of simulations + 1: 10001 

statistic = 0.1982, observed rank = 10000, p-value = 9.999e-05
alternative hypothesis: greater

MexicoCity_HTS$Variable <- MexicoCity_HTS$Residuals_OLS
ZonesAnalysis <- MexicoCity_HTS %>% select(Variable)


LocalMoran_Temp <- localmoran(ZonesAnalysis$Variable, WM_)

LocalMoran_Temp <- as.data.frame(LocalMoran_Temp)
names(LocalMoran_Temp)[5] <- "PValue"

ZonesAnalysis %<>% bind_cols(LocalMoran_Temp)

ZonesAnalysis$Lagged <- lag.listw(WM_, ZonesAnalysis$Variable)
ZonesAnalysis$Value <- ZonesAnalysis$Variable

MeanValues = mean(ZonesAnalysis$Value)
MeanLagged = mean(ZonesAnalysis$Lagged)
MeanMoran = mean(ZonesAnalysis$Ii)

ZonesAnalysis %<>%
  mutate(
    Values_Centered = Value - MeanValues,
    Lagged_Centered = Lagged - MeanLagged,
    Moran_Centered = Ii - MeanMoran)

ZonesAnalysis %<>%
  mutate(Significance = if_else(PValue <= 0.05, 1, 0)) %>% 
  mutate(Cuadrants = 5) %>%
  mutate(Cuadrants = if_else(Values_Centered > 0 & Lagged_Centered > 0, 1, Cuadrants)) %>% 
  mutate(Cuadrants = if_else(Values_Centered > 0 & Lagged_Centered < 0, 2, Cuadrants)) %>%
  mutate(Cuadrants = if_else(Values_Centered < 0 & Lagged_Centered < 0, 3, Cuadrants)) %>%
  mutate(Cuadrants = if_else(Values_Centered < 0 & Lagged_Centered > 0, 4, Cuadrants)) %>% 
  mutate(Cuadrants = if_else(Significance == 0, 5, Cuadrants))

# 1: High-High - Spatial clusters
# 2: High-Low - Outliers
# 3: Low-low - Spatial clusters
# 4: Low-High - Outliers
# 5: Not significance
# c("High-High", "High-Low", "Low-low", "Low-high", "Not significant")
ZonesAnalysis$Cuadrants <- factor(as.character(ZonesAnalysis$Cuadrants), 
                                  levels = c("1", "2", "3", "4", "5"),
                                  ordered = TRUE)

Breaks <- c(0, 1.5, 2.5, 3.5, 4.5, 5.5)
Labels <- c("High-high (clusters)", "High-low", 
            "Low-low (clusters)", "Low-high", "Not significant")
MyPalette <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a", "#edf8fb")

Map_LISA_Residuals <- tm_shape(ZonesAnalysis, bbox = BB_Districts) +
  tm_fill("Cuadrants",
           palette = MyPalette, breaks = Breaks, labels = Labels,
          title = "LISA cluster map") +
  tm_borders(alpha = 0.4) +
  tm_shape(MexicoCity_Fringe) + 
  tm_polygons(alpha = 0, border.col = "black", lwd = 1.4, lty = "solid") +
  tm_compass(size = 2, type = "arrow", position = c(0.85,0.87)) + 
  tm_scale_bar(position = c(0.4,0.02)) +
  tm_layout(title = "Residuals - OLS", legend.position = c(0.01,0.15), scale = 1)

Map_LISA_Residuals

4 OLS - What is next?

4.1 Anselin test

lm.LMtests(Model_OLS, WM_, 
           test = c("LMerr", "LMlag", "RLMerr", "RLMlag", "SARMA"))

    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = ModelEquation, data = MexicoCity_HTS)
weights: WM_

LMerr = 19.206, df = 1, p-value = 1.173e-05


    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = ModelEquation, data = MexicoCity_HTS)
weights: WM_

LMlag = 12.008, df = 1, p-value = 0.0005298


    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = ModelEquation, data = MexicoCity_HTS)
weights: WM_

RLMerr = 7.2671, df = 1, p-value = 0.007023


    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = ModelEquation, data = MexicoCity_HTS)
weights: WM_

RLMlag = 0.068685, df = 1, p-value = 0.7933


    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = ModelEquation, data = MexicoCity_HTS)
weights: WM_

SARMA = 19.275, df = 2, p-value = 6.524e-05

4.2 Likelihood ratio tests

To do this we neet to fit the other models

# Spatial Durbin Model SDM:
Spatial_SDM <- lagsarlm(ModelEquation, data = MexicoCity_HTS, WM_, type = "mixed") 

# Spatial Durbin Error Model SDEM:
Spatial_SDEM <- errorsarlm(ModelEquation, data = MexicoCity_HTS, WM_, etype="emixed") 

# Spatially Lagged X:
Spatial_Lagged_X <- lmSLX(ModelEquation, data = MexicoCity_HTS, WM_) 

# Spatially Lagged Y:
Spatial_Lagged_Y <- lagsarlm(ModelEquation, data = MexicoCity_HTS, WM_) 

# Spatial Error Model SEM:
Spatial_SEM <- errorsarlm(ModelEquation, data = MexicoCity_HTS, WM_) 

5 SAR (SLY)

Spatial_Lagged_Y <- lagsarlm(ModelEquation, data = MexicoCity_HTS, WM_) 

# summary(impacts(Spatial_Lagged_Y, listw = WM_, R = 500), zstats = TRUE)
summary(impacts(Spatial_Lagged_Y, tr = trMC, R = 300), zstats = TRUE)

Impact measures (lag, trace):
                            Direct    Indirect       Total
Prop_Male              -0.40482489 -0.12715291 -0.53197780
Prop_HomesWithCars      1.78321869  0.56009758  2.34331626
Prop_HomesWithBicycles  0.08404907  0.02639928  0.11044835
Age                     0.03719539  0.01168284  0.04887823
========================================================
Simulation results ( variance matrix):
Direct:

Iterations = 1:300
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 300 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                           Mean       SD  Naive SE Time-series SE
Prop_Male              -0.28773 1.686603 0.0973761      0.0973761
Prop_HomesWithCars      1.79394 0.193549 0.0111745      0.0111745
Prop_HomesWithBicycles  0.08606 0.170480 0.0098426      0.0111856
Age                     0.03745 0.009115 0.0005262      0.0005262

2. Quantiles for each variable:

                           2.5%      25%      50%     75%  97.5%
Prop_Male              -3.54647 -1.45475 -0.32184 0.85809 3.3606
Prop_HomesWithCars      1.43983  1.64359  1.81281 1.92369 2.1495
Prop_HomesWithBicycles -0.30879 -0.01406  0.08395 0.19900 0.4166
Age                     0.01907  0.03119  0.03736 0.04392 0.0559

========================================================
Indirect:

Iterations = 1:300
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 300 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                           Mean      SD  Naive SE Time-series SE
Prop_Male              -0.07432 0.59820 0.0345371      0.0345371
Prop_HomesWithCars      0.57754 0.24104 0.0139165      0.0139165
Prop_HomesWithBicycles  0.03266 0.06230 0.0035967      0.0035967
Age                     0.01147 0.00416 0.0002402      0.0003261

2. Quantiles for each variable:

                            2.5%       25%      50%     75%   97.5%
Prop_Male              -1.351206 -0.408711 -0.08245 0.28791 1.32733
Prop_HomesWithCars      0.198162  0.414494  0.55670 0.70050 1.15753
Prop_HomesWithBicycles -0.084967 -0.003998  0.02424 0.05967 0.18582
Age                     0.005002  0.008232  0.01115 0.01404 0.02026

========================================================
Total:

Iterations = 1:300
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 300 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                           Mean      SD  Naive SE Time-series SE
Prop_Male              -0.36206 2.25198 0.1300181      0.1300181
Prop_HomesWithCars      2.37148 0.33410 0.0192891      0.0214045
Prop_HomesWithBicycles  0.11872 0.22805 0.0131665      0.0150203
Age                     0.04893 0.01027 0.0005927      0.0005927

2. Quantiles for each variable:

                           2.5%      25%      50%     75%  97.5%
Prop_Male              -4.43430 -1.89834 -0.39078 1.22506 4.4790
Prop_HomesWithCars      1.79700  2.14015  2.34053 2.58610 3.0714
Prop_HomesWithBicycles -0.37590 -0.01901  0.11019 0.27767 0.5710
Age                     0.02768  0.04276  0.04908 0.05562 0.0674

========================================================
Simulated standard errors
                            Direct    Indirect      Total
Prop_Male              1.686603140 0.598199365 2.25197957
Prop_HomesWithCars     0.193548815 0.241041095 0.33409742
Prop_HomesWithBicycles 0.170479525 0.062296902 0.22805046
Age                    0.009114664 0.004159547 0.01026658

Simulated z-values:
                           Direct   Indirect      Total
Prop_Male              -0.1705999 -0.1242446 -0.1607729
Prop_HomesWithCars      9.2686660  2.3960422  7.0981809
Prop_HomesWithBicycles  0.5048339  0.5242102  0.5205888
Age                     4.1087870  2.7586925  4.7654762

Simulated p-values:
                       Direct     Indirect  Total     
Prop_Male              0.86454    0.9011216 0.87227   
Prop_HomesWithCars     < 2.22e-16 0.0165732 1.2641e-12
Prop_HomesWithBicycles 0.61368    0.6001323 0.60265   
Age                    3.9774e-05 0.0058033 1.8841e-06

MexicoCity_HTS$Residuals_SLY <- residuals(Spatial_Lagged_Y)
moran.mc(MexicoCity_HTS$Residuals_SLY, WM_, 1000)

    Monte-Carlo simulation of Moran I

data:  MexicoCity_HTS$Residuals_SLY 
weights: WM_  
number of simulations + 1: 1001 

statistic = 0.071464, observed rank = 973, p-value = 0.02797
alternative hypothesis: greater

6 SEM

Spatial_SEM <- errorsarlm(ModelEquation, data = MexicoCity_HTS, WM_)

summary(Spatial_SEM)

Call:errorsarlm(formula = ModelEquation, data = MexicoCity_HTS, listw = WM_)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.417042 -0.125463 -0.037412  0.066539  2.347826 

Type: error 
Coefficients: (asymptotic standard errors) 
                        Estimate Std. Error z value  Pr(>|z|)
(Intercept)            -1.430771   0.855190 -1.6730   0.09432
Prop_Male              -0.683613   1.468913 -0.4654   0.64165
Prop_HomesWithCars      1.921347   0.224748  8.5489 < 2.2e-16
Prop_HomesWithBicycles -0.114785   0.205418 -0.5588   0.57631
Age                     0.044735   0.010148  4.4082 1.042e-05

Lambda: 0.40122, LR test value: 16.39, p-value: 5.1565e-05
Asymptotic standard error: 0.093219
    z-value: 4.3041, p-value: 1.677e-05
Wald statistic: 18.525, p-value: 1.677e-05

Log likelihood: -21.76592 for error model
ML residual variance (sigma squared): 0.070907, (sigma: 0.26628)
Number of observations: 193 
Number of parameters estimated: 7 
AIC: 57.532, (AIC for lm: 71.922)

MexicoCity_HTS$Residuals_SEM <- residuals(Spatial_SEM)
moran.mc(MexicoCity_HTS$Residuals_SEM, WM_, 1000)

    Monte-Carlo simulation of Moran I

data:  MexicoCity_HTS$Residuals_SEM 
weights: WM_  
number of simulations + 1: 1001 

statistic = 0.0017268, observed rank = 595, p-value = 0.4056
alternative hypothesis: greater

7 SLX

Spatial_Lagged_X <- lmSLX(ModelEquation, data = MexicoCity_HTS, WM_) 

summary(Spatial_Lagged_X)

Call:
lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
    data = as.data.frame(x), weights = weights)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.44924 -0.13976 -0.04251  0.07882  2.35127 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                -4.76320    2.25270  -2.114  0.03582 *  
Prop_Male                   0.07938    1.62872   0.049  0.96118    
Prop_HomesWithCars          1.98123    0.26881   7.370 5.58e-12 ***
Prop_HomesWithBicycles     -0.24415    0.28036  -0.871  0.38498    
Age                         0.04222    0.01455   2.901  0.00417 ** 
lag.Prop_Male               4.40362    3.56824   1.234  0.21873    
lag.Prop_HomesWithCars     -0.19219    0.41800  -0.460  0.64621    
lag.Prop_HomesWithBicycles  0.39702    0.36334   1.093  0.27595    
lag.Age                     0.02659    0.01917   1.387  0.16702    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2867 on 184 degrees of freedom
Multiple R-squared:  0.6569,    Adjusted R-squared:  0.642 
F-statistic: 44.04 on 8 and 184 DF,  p-value: < 2.2e-16

#summary(impacts(Spatial_Lagged_X, listw = WM, R=500), zstats = TRUE) 
summary(impacts(Spatial_Lagged_X, tr = trMC, R = 300), zstats = TRUE)

Impact measures (SlX, estimable, n-k):
                            Direct    Indirect      Total
Prop_Male               0.07937935  4.40361829 4.48299764
Prop_HomesWithCars      1.98123002 -0.19219349 1.78903653
Prop_HomesWithBicycles -0.24414664  0.39702087 0.15287424
Age                     0.04221822  0.02659187 0.06881009
========================================================
Standard errors:
                          Direct   Indirect      Total
Prop_Male              1.6287178 3.56823526 4.00323060
Prop_HomesWithCars     0.2688123 0.41800191 0.33389169
Prop_HomesWithBicycles 0.2803596 0.36333921 0.23971575
Age                    0.0145522 0.01916755 0.01380898
========================================================
Z-values:
                            Direct   Indirect     Total
Prop_Male               0.04873733  1.2341166 1.1198450
Prop_HomesWithCars      7.37030960 -0.4597909 5.3581344
Prop_HomesWithBicycles -0.87083375  1.0927003 0.6377313
Age                     2.90115850  1.3873377 4.9829962

p-values:
                       Direct     Indirect Total     
Prop_Male              0.9611286  0.21716  0.26278   
Prop_HomesWithCars     1.7031e-13 0.64567  8.4086e-08
Prop_HomesWithBicycles 0.3838449  0.27453  0.52365   
Age                    0.0037179  0.16534  6.2607e-07

MexicoCity_HTS$Residuals_SLX <- residuals(Spatial_Lagged_X)
moran.mc(MexicoCity_HTS$Residuals_SLX, WM_, 1000)

    Monte-Carlo simulation of Moran I

data:  MexicoCity_HTS$Residuals_SLX 
weights: WM_  
number of simulations + 1: 1001 

statistic = 0.20062, observed rank = 1001, p-value = 0.000999
alternative hypothesis: greater

8 SDM

Spatial_SDM <- lagsarlm(ModelEquation, data = MexicoCity_HTS, WM_, type = "mixed") 

summary(Spatial_SDM)

Call:lagsarlm(formula = ModelEquation, data = MexicoCity_HTS, listw = WM_, 
    type = "mixed")

Residuals:
      Min        1Q    Median        3Q       Max 
-0.432570 -0.126697 -0.044845  0.070641  2.303348 

Type: mixed 
Coefficients: (asymptotic standard errors) 
                             Estimate Std. Error z value  Pr(>|z|)
(Intercept)                -3.8249906  2.0962277 -1.8247  0.068046
Prop_Male                  -0.1814561  1.5018323 -0.1208  0.903831
Prop_HomesWithCars          1.9702295  0.2474804  7.9612 1.776e-15
Prop_HomesWithBicycles     -0.3033832  0.2581791 -1.1751  0.239959
Age                         0.0422458  0.0134022  3.1522  0.001621
lag.Prop_Male               4.6306169  3.2869657  1.4088  0.158900
lag.Prop_HomesWithCars     -0.8992244  0.4201823 -2.1401  0.032348
lag.Prop_HomesWithBicycles  0.4401547  0.3344791  1.3159  0.188194
lag.Age                     0.0021121  0.0183514  0.1151  0.908372

Rho: 0.39539, LR test value: 16.363, p-value: 5.2283e-05
Asymptotic standard error: 0.093311
    z-value: 4.2373, p-value: 2.2622e-05
Wald statistic: 17.955, p-value: 2.2622e-05

Log likelihood: -19.97021 for mixed model
ML residual variance (sigma squared): 0.069674, (sigma: 0.26396)
Number of observations: 193 
Number of parameters estimated: 11 
AIC: 61.94, (AIC for lm: 76.304)
LM test for residual autocorrelation
test value: 0.038034, p-value: 0.84537

# summary(impacts(Spatial_SDM, listw = WM, R = 500), zstats = TRUE)
summary(impacts(Spatial_SDM, tr = trMC, R = 300), zstats = TRUE)

Impact measures (mixed, trace):
                            Direct    Indirect      Total
Prop_Male               0.22988190  7.12883141 7.35871331
Prop_HomesWithCars      1.95938245 -0.18798782 1.77139463
Prop_HomesWithBicycles -0.27449212  0.50070605 0.22621394
Age                     0.04394349  0.02942249 0.07336598
========================================================
Simulation results ( variance matrix):
Direct:

Iterations = 1:300
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 300 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                           Mean      SD  Naive SE Time-series SE
Prop_Male               0.19784 1.67720 0.0968334       0.096833
Prop_HomesWithCars      1.95145 0.22321 0.0128869       0.012887
Prop_HomesWithBicycles -0.25745 0.22575 0.0130335       0.013034
Age                     0.04438 0.01196 0.0006905       0.000794

2. Quantiles for each variable:

                           2.5%      25%      50%      75%   97.5%
Prop_Male              -3.00712 -0.95370  0.14732  1.26864 3.59994
Prop_HomesWithCars      1.52177  1.80618  1.95746  2.09797 2.39055
Prop_HomesWithBicycles -0.69425 -0.41385 -0.25665 -0.09375 0.14178
Age                     0.01945  0.03671  0.04522  0.05223 0.06738

========================================================
Indirect:

Iterations = 1:300
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 300 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                           Mean      SD Naive SE Time-series SE
Prop_Male               7.13448 5.66967 0.327339       0.327339
Prop_HomesWithCars     -0.15171 0.57295 0.033079       0.029871
Prop_HomesWithBicycles  0.48314 0.42863 0.024747       0.024747
Age                     0.02826 0.02552 0.001473       0.001473

2. Quantiles for each variable:

                           2.5%       25%      50%      75%    97.5%
Prop_Male              -3.79847  3.441348  7.09681 10.86772 18.39668
Prop_HomesWithCars     -1.17876 -0.491039 -0.16342  0.17461  0.99895
Prop_HomesWithBicycles -0.36278  0.192500  0.50568  0.77079  1.39399
Age                    -0.01521  0.009882  0.02793  0.04431  0.07934

========================================================
Total:

Iterations = 1:300
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 300 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                          Mean      SD Naive SE Time-series SE
Prop_Male              7.33232 6.56947 0.379289       0.379289
Prop_HomesWithCars     1.79974 0.55816 0.032226       0.032226
Prop_HomesWithBicycles 0.22570 0.39924 0.023050       0.023050
Age                    0.07264 0.02349 0.001356       0.001356

2. Quantiles for each variable:

                           2.5%      25%     50%      75%   97.5%
Prop_Male              -4.82264  3.03100 6.97300 11.83149 19.7914
Prop_HomesWithCars      0.80686  1.45521 1.78276  2.15513  2.9270
Prop_HomesWithBicycles -0.54150 -0.01329 0.21205  0.49629  1.0030
Age                     0.02988  0.05712 0.07381  0.08597  0.1175

========================================================
Simulated standard errors
                           Direct   Indirect      Total
Prop_Male              1.67720321 5.66967274 6.56947251
Prop_HomesWithCars     0.22320703 0.57294737 0.55816300
Prop_HomesWithBicycles 0.22574728 0.42862829 0.39923845
Age                    0.01196058 0.02551567 0.02348687

Simulated z-values:
                           Direct   Indirect     Total
Prop_Male               0.1179592  1.2583585 1.1161204
Prop_HomesWithCars      8.7427758 -0.2647949 3.2243904
Prop_HomesWithBicycles -1.1404232  1.1271854 0.5653166
Age                     3.7106577  1.1074642 3.0927632

Simulated p-values:
                       Direct     Indirect Total    
Prop_Male              0.90609995 0.20826  0.2643706
Prop_HomesWithCars     < 2.22e-16 0.79117  0.0012624
Prop_HomesWithBicycles 0.25411004 0.25966  0.5718584
Age                    0.00020672 0.26809  0.0019830

MexicoCity_HTS$SDM <- residuals(Spatial_SDM)
moran.mc(MexicoCity_HTS$SDM, WM_, 1000)

    Monte-Carlo simulation of Moran I

data:  MexicoCity_HTS$SDM 
weights: WM_  
number of simulations + 1: 1001 

statistic = 0.0017916, observed rank = 603, p-value = 0.3976
alternative hypothesis: greater

9 SDEM

Spatial_SDEM <- errorsarlm(ModelEquation, data = MexicoCity_HTS, WM_, etype="emixed")

summary(Spatial_SDEM)

Call:errorsarlm(formula = ModelEquation, data = MexicoCity_HTS, listw = WM_, 
    etype = "emixed")

Residuals:
      Min        1Q    Median        3Q       Max 
-0.413021 -0.132669 -0.042050  0.074065  2.305038 

Type: error 
Coefficients: (asymptotic standard errors) 
                            Estimate Std. Error z value  Pr(>|z|)
(Intercept)                -5.006409   2.587116 -1.9351 0.0529742
Prop_Male                   0.165953   1.584209  0.1048 0.9165704
Prop_HomesWithCars          1.933859   0.235438  8.2139  2.22e-16
Prop_HomesWithBicycles     -0.270395   0.243243 -1.1116 0.2662996
Age                         0.045224   0.012534  3.6082 0.0003083
lag.Prop_Male               5.076083   3.743991  1.3558 0.1751645
lag.Prop_HomesWithCars     -0.076900   0.444190 -0.1731 0.8625541
lag.Prop_HomesWithBicycles  0.426103   0.356293  1.1959 0.2317227
lag.Age                     0.019061   0.019066  0.9997 0.3174376

Lambda: 0.39509, LR test value: 16.33, p-value: 5.3211e-05
Asymptotic standard error: 0.093687
    z-value: 4.2172, p-value: 2.4739e-05
Wald statistic: 17.784, p-value: 2.4739e-05

Log likelihood: -19.98688 for error model
ML residual variance (sigma squared): 0.06969, (sigma: 0.26399)
Number of observations: 193 
Number of parameters estimated: 11 
AIC: 61.974, (AIC for lm: 76.304)

summary(Spatial_SDEM)

Call:errorsarlm(formula = ModelEquation, data = MexicoCity_HTS, listw = WM_, 
    etype = "emixed")

Residuals:
      Min        1Q    Median        3Q       Max 
-0.413021 -0.132669 -0.042050  0.074065  2.305038 

Type: error 
Coefficients: (asymptotic standard errors) 
                            Estimate Std. Error z value  Pr(>|z|)
(Intercept)                -5.006409   2.587116 -1.9351 0.0529742
Prop_Male                   0.165953   1.584209  0.1048 0.9165704
Prop_HomesWithCars          1.933859   0.235438  8.2139  2.22e-16
Prop_HomesWithBicycles     -0.270395   0.243243 -1.1116 0.2662996
Age                         0.045224   0.012534  3.6082 0.0003083
lag.Prop_Male               5.076083   3.743991  1.3558 0.1751645
lag.Prop_HomesWithCars     -0.076900   0.444190 -0.1731 0.8625541
lag.Prop_HomesWithBicycles  0.426103   0.356293  1.1959 0.2317227
lag.Age                     0.019061   0.019066  0.9997 0.3174376

Lambda: 0.39509, LR test value: 16.33, p-value: 5.3211e-05
Asymptotic standard error: 0.093687
    z-value: 4.2172, p-value: 2.4739e-05
Wald statistic: 17.784, p-value: 2.4739e-05

Log likelihood: -19.98688 for error model
ML residual variance (sigma squared): 0.06969, (sigma: 0.26399)
Number of observations: 193 
Number of parameters estimated: 11 
AIC: 61.974, (AIC for lm: 76.304)

MexicoCity_HTS$SDEM <- residuals(Spatial_SDEM)
moran.mc(MexicoCity_HTS$SDEM, WM_, 1000)

    Monte-Carlo simulation of Moran I

data:  MexicoCity_HTS$SDEM 
weights: WM_  
number of simulations + 1: 1001 

statistic = 0.00040948, observed rank = 580, p-value = 0.4206
alternative hypothesis: greater

10 Likelihood ratio test (Lagrange Multipliers)

LR.Sarlm(Spatial_SDM, Spatial_Lagged_X)

    Likelihood ratio for spatial linear models

data:  
Likelihood ratio = 16.363, df = 1, p-value = 5.228e-05
sample estimates:
     Log likelihood of Spatial_SDM Log likelihood of Spatial_Lagged_X 
                         -19.97021                          -28.15195 

LR.Sarlm(Spatial_SDM, Spatial_Lagged_Y)

    Likelihood ratio for spatial linear models

data:  
Likelihood ratio = 9.6249, df = 4, p-value = 0.04724
sample estimates:
     Log likelihood of Spatial_SDM Log likelihood of Spatial_Lagged_Y 
                         -19.97021                          -24.78265 

LR.Sarlm(Spatial_SDM, Spatial_SEM)

    Likelihood ratio for spatial linear models

data:  
Likelihood ratio = 3.5914, df = 4, p-value = 0.4641
sample estimates:
Log likelihood of Spatial_SDM Log likelihood of Spatial_SEM 
                    -19.97021                     -21.76592 

LR.Sarlm(Spatial_SDEM, Spatial_SEM)

    Likelihood ratio for spatial linear models

data:  
Likelihood ratio = 3.5581, df = 4, p-value = 0.4691
sample estimates:
Log likelihood of Spatial_SDEM  Log likelihood of Spatial_SEM 
                     -19.98688                      -21.76592 

LR.Sarlm(Spatial_SDEM, Spatial_Lagged_X)

    Likelihood ratio for spatial linear models

data:  
Likelihood ratio = 16.33, df = 1, p-value = 5.321e-05
sample estimates:
    Log likelihood of Spatial_SDEM Log likelihood of Spatial_Lagged_X 
                         -19.98688                          -28.15195