4  Modelos espaciais de corte transversal (cross-section)

4.1 Pacotes

library(spatialreg)

Loading required package: spData

The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
     (status 2 uses the sf package in place of rgdal)

To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`

Loading required package: Matrix

Loading required package: sf

Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

library(spdep)

Attaching package: 'spdep'

The following objects are masked from 'package:spatialreg':

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

library(sphet)

Attaching package: 'sphet'

The following object is masked from 'package:spatialreg':

    impacts

4.2 Shapefile

# Pacotes
library(sf)
library(sp)

# Abra o arquivo
montana.shp <- st_read("data/natregimes_montana.gml")

Reading layer `natregimes_montana' from data source 
  `/home/raphael/projects/ecoespacial/data/natregimes_montana.gml' 
  using driver `GML'
Simple feature collection with 55 features and 74 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -116.0625 ymin: 44.35373 xmax: -104.0426 ymax: 49
Geodetic CRS:  WGS 84

montana.shp <- st_make_valid(montana.shp)
montana.shp <- as_Spatial(montana.shp)

# Plotar o mapa
plot(montana.shp)

library(leaflet)
library(RColorBrewer)
library(htmltools)

qpal <- colorNumeric("OrRd", montana.shp@data$HR90, n=6) 

leaflet(montana.shp) %>%
  addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2,
              fillColor = ~qpal(HR90), popup = ~htmlEscape(paste0(NAME, ": ", round(HR90, 2))),
              highlight = highlightOptions(fillOpacity = 0.5,
                                           bringToFront = TRUE)) %>%
  addTiles()

4.3 Matriz de vizinhança

w1 <- spdep::nb2listw(spdep::poly2nb(montana.shp, queen = FALSE), style = "W")
W <- as(w1, "CsparseMatrix")
trMat <- trW(W, type="mult")

4.4 Especificação do modelo

esp1 <- HR90 ~ RD90
endog <- ~ UE90
instruments <- ~ FH90

• HR: homicide rate per 100,000 (1960, 1970, 1980, 1990)
• RD: resource deprivation 1960, 1970, 1980, 1990 (principal component)
• UE: unemployment rate 1960, 1970, 1980, 1990
• FH: % female headed households 1960, 1970, 1980, 1990

4.5 OLS

\(y = \alpha + X\beta + \varepsilon\)

mod1.mcrl <- lm(formula = esp1, data = montana.shp@data)
summary(mod1.mcrl)

Call:
lm(formula = esp1, data = montana.shp@data)

Residuals:
   Min     1Q Median     3Q    Max 
-6.546 -3.017 -1.025  2.422 16.951 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.7853     0.6129   6.176 9.49e-08 ***
RD90          2.8836     1.0923   2.640   0.0109 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.433 on 53 degrees of freedom
Multiple R-squared:  0.1162,    Adjusted R-squared:  0.09954 
F-statistic: 6.969 on 1 and 53 DF,  p-value: 0.01087

4.6 Testes do Multiplicador de Lagrange

mod1.lagrange <- lm.LMtests(model = mod1.mcrl, listw = w1,
                            test = c("LMerr","RLMerr","LMlag","RLMlag"))
mod1.lagrange

    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = esp1, data = montana.shp@data)
weights: w1

LMerr = 0.9605, df = 1, p-value = 0.3271


    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = esp1, data = montana.shp@data)
weights: w1

RLMerr = 0.13243, df = 1, p-value = 0.7159


    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = esp1, data = montana.shp@data)
weights: w1

LMlag = 0.84352, df = 1, p-value = 0.3584


    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = esp1, data = montana.shp@data)
weights: w1

RLMlag = 0.015448, df = 1, p-value = 0.9011

4.7 SAR (MV)

\(y = \rho W y + X \beta + \varepsilon\)

mod1.sar <- lagsarlm(formula = esp1, data = montana.shp@data, listw = w1)
summary(mod1.sar)

Call:lagsarlm(formula = esp1, data = montana.shp@data, listw = w1)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4712 -3.0399 -1.3250  2.5259 16.8433 

Type: lag 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept)  3.16687    0.85033  3.7243 0.0001959
RD90         2.91026    1.06059  2.7440 0.0060695

Rho: 0.1792, LR test value: 0.84239, p-value: 0.35871
Asymptotic standard error: 0.18437
    z-value: 0.97196, p-value: 0.33107
Wald statistic: 0.94472, p-value: 0.33107

Log likelihood: -158.4951 for lag model
ML residual variance (sigma squared): 18.519, (sigma: 4.3034)
Number of observations: 55 
Number of parameters estimated: 4 
AIC: 324.99, (AIC for lm: 323.83)
LM test for residual autocorrelation
test value: 0.54172, p-value: 0.46172

summary(spatialreg::impacts(mod1.sar, tr=trMat, R=1000), zstats=TRUE, short=TRUE)

Impact measures (lag, trace):
       Direct  Indirect    Total
RD90 2.930558 0.6150693 3.545627
========================================================
Simulation results ( variance matrix):
========================================================
Simulated standard errors
       Direct Indirect    Total
RD90 1.123352 1.093339 1.842267

Simulated z-values:
       Direct Indirect    Total
RD90 2.556033 0.761637 2.010594

Simulated p-values:
     Direct   Indirect Total   
RD90 0.010587 0.44628  0.044368

4.8 SAR (STSLS)

mod1.sar_stsls <- spreg(
  formula = esp1, data = montana.shp@data, listw = w1,
  model = "lag", het = TRUE, endog = endog, instruments = instruments
)
summary(mod1.sar_stsls)

 Generalized stsls

Call:
spreg(formula = esp1, data = montana.shp@data, listw = w1, endog = endog, 
    instruments = instruments, model = "lag", het = TRUE)

Residuals:
    Min.  1st Qu.   Median  3rd Qu.     Max. 
-6.48374 -2.44405 -0.88756  1.14961 16.99149 

Coefficients:
            Estimate Std. Error t-value  Pr(>|t|)    
(Intercept) -1.22654    2.07993 -0.5897 0.5553893    
RD90         0.25793    1.11656  0.2310 0.8173100    
UE90         0.74020    0.20802  3.5583 0.0003733 ***
lambda       0.05640    0.58238  0.0968 0.9228503    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mod1.sac_stsls2 <- spreg(
  formula = esp1, data = montana.shp@data, listw = w1,
  model = "lag", het = TRUE
)
summary(mod1.sac_stsls2)

 Generalized stsls

Call:
spreg(formula = esp1, data = montana.shp@data, listw = w1, model = "lag", 
    het = TRUE)

Residuals:
   Min. 1st Qu.  Median 3rd Qu.    Max. 
-6.5149 -3.0147 -1.1715  2.4648 16.9064 

Coefficients:
            Estimate Std. Error t-value Pr(>|t|)  
(Intercept) 3.529675   2.405383  1.4674  0.14227  
RD90        2.894640   1.215182  2.3821  0.01722 *
lambda      0.074062   0.613635  0.1207  0.90393  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(sphet::impacts(mod1.sac_stsls2, tr=trMat, R=1000), zstats=TRUE, short=TRUE)

Impact measures (lag, trace):
       Direct  Indirect    Total
RD90 2.897931 0.2282375 3.126169
========================================================
Simulation results (IV variance matrix):
========================================================
Simulated standard errors
       Direct Indirect    Total
RD90 1.361766 7.614692 8.265529

Simulated z-values:
       Direct  Indirect     Total
RD90 2.318299 0.2932021 0.6520604

Simulated p-values:
     Direct   Indirect Total  
RD90 0.020433 0.76937  0.51436

4.9 SEM (MV)

\(y = X \beta + u \\ u = \lambda Wu + \varepsilon\)

mod1.sem <- errorsarlm(formula = esp1, data = montana.shp@data, listw = w1)
summary(mod1.sem)

Call:errorsarlm(formula = esp1, data = montana.shp@data, listw = w1)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4097 -3.0084 -1.1481  2.5629 16.6659 

Type: error 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept)  3.75130    0.73537  5.1012 3.375e-07
RD90         2.89021    1.03890  2.7820  0.005402

Lambda: 0.20192, LR test value: 0.98506, p-value: 0.32095
Asymptotic standard error: 0.18804
    z-value: 1.0738, p-value: 0.28291
Wald statistic: 1.153, p-value: 0.28291

Log likelihood: -158.4237 for error model
ML residual variance (sigma squared): 18.436, (sigma: 4.2937)
Number of observations: 55 
Number of parameters estimated: 4 
AIC: NA (not available for weighted model), (AIC for lm: 323.83)

4.10 SEM (GMM)

mod1.sem_gmm <- spreg(
  formula = esp1, data = montana.shp@data, listw = w1,
  model = "error", step1.c = FALSE, het = TRUE, endog = endog, instruments = instruments
)

summary(mod1.sem_gmm)

 Generalized stsls

Call:
spreg(formula = esp1, data = montana.shp@data, listw = w1, endog = endog, 
    instruments = instruments, model = "error", het = TRUE, step1.c = FALSE)

Residuals:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-6.4079 -2.4302 -0.8923  0.0002  1.0845 17.0259 

Coefficients:
             Estimate Std. Error t-value  Pr(>|t|)    
(Intercept) -1.052436   1.204375 -0.8738 0.3822034    
RD90         0.238720   1.062023  0.2248 0.8221516    
UE90         0.743321   0.198001  3.7541 0.0001739 ***
rho          0.071436   0.303046  0.2357 0.8136454    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.11 SAC (MV)

\(y = \rho Wy + X \beta + u \\ u = \lambda Wu + \varepsilon\)

mod1.sac <- sacsarlm(formula = esp1, data = montana.shp@data, listw = w1)
summary(mod1.sac)

Call:sacsarlm(formula = esp1, data = montana.shp@data, listw = w1)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4318 -3.0064 -1.2061  2.5673 16.7131 

Type: sac 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   3.5980     2.7176  1.3240 0.185507
RD90          2.9052     1.0866  2.6736 0.007505

Rho: 0.047078
Asymptotic standard error: 0.8077
    z-value: 0.058287, p-value: 0.95352
Lambda: 0.16245
Asymptotic standard error: 0.79051
    z-value: 0.2055, p-value: 0.83718

LR test value: 0.99918, p-value: 0.60678

Log likelihood: -158.4167 for sac model
ML residual variance (sigma squared): 18.481, (sigma: 4.2989)
Number of observations: 55 
Number of parameters estimated: 5 
AIC: 326.83, (AIC for lm: 323.83)

summary(spatialreg::impacts(mod1.sac, tr=trMat, R=1000), zstats=TRUE, short=TRUE)

Impact measures (sac, trace):
      Direct  Indirect    Total
RD90 2.90656 0.1422079 3.048768
========================================================
Simulation results ( variance matrix):
========================================================
Simulated standard errors
       Direct Indirect   Total
RD90 1.283029 10.06036 10.8258

Simulated z-values:
       Direct  Indirect     Total
RD90 2.433905 0.2743424 0.5434012

Simulated p-values:
     Direct   Indirect Total  
RD90 0.014937 0.78382  0.58685

4.12 SAC (GMM)

mod1.sac_gmm <- spreg(
  formula = esp1, data = montana.shp@data, listw = w1,
  model = "sarar", step1.c = FALSE, het = TRUE, endog = endog, instruments = instruments
)

summary(mod1.sac_gmm)

 Generalized stsls

Call:
spreg(formula = esp1, data = montana.shp@data, listw = w1, endog = endog, 
    instruments = instruments, model = "sarar", het = TRUE, step1.c = FALSE)

Residuals:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-6.5102 -2.4510 -0.8904 -0.0013  1.1410 16.9844 

Coefficients:
             Estimate Std. Error t-value  Pr(>|t|)    
(Intercept) -1.282437   2.519760 -0.5090 0.6107858    
RD90         0.240835   1.194968  0.2015 0.8402758    
UE90         0.743794   0.207021  3.5928 0.0003271 ***
lambda       0.065972   0.627655  0.1051 0.9162890    
rho          0.302036   0.557139  0.5421 0.5877362    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wald test that rho and lambda are both zero:
 Statistics: 0.65153 p-val: 0.41956 

mod1.sac_gmm2 <- spreg(
  formula = esp1, data = montana.shp@data, listw = w1,
  model = "sarar", step1.c = FALSE, het = TRUE
)

summary(mod1.sac_gmm2)

 Generalized stsls

Call:
spreg(formula = esp1, data = montana.shp@data, listw = w1, model = "sarar", 
    het = TRUE, step1.c = FALSE)

Residuals:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -6.521  -3.003  -1.089   0.021   2.463  16.953 

Coefficients:
            Estimate Std. Error t-value Pr(>|t|)  
(Intercept) 3.651563   2.415507  1.5117  0.13061  
RD90        2.897586   1.212861  2.3890  0.01689 *
lambda      0.032985   0.615118  0.0536  0.95723  
rho         0.012935   0.619914  0.0209  0.98335  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wald test that rho and lambda are both zero:
 Statistics: 0.040424 p-val: 0.84065 

summary(sphet::impacts(mod1.sac_gmm2, tr=trMat, R=1000), zstats=TRUE, short=TRUE)

Impact measures (lag, trace):
       Direct   Indirect    Total
RD90 2.898229 0.09819482 2.996424
========================================================
Simulation results (GSTSLS variance matrix):
========================================================
Simulated standard errors
       Direct Indirect    Total
RD90 1.399784 9.161207 9.879057

Simulated z-values:
       Direct Indirect     Total
RD90 2.265946 0.253405 0.5560581

Simulated p-values:
     Direct   Indirect Total  
RD90 0.023455 0.79996  0.57817

4.13 SLX (MV)

\(y = X \beta + WX \theta + \varepsilon\)

mod1.slx <- lmSLX(formula = esp1, data = montana.shp@data, listw = w1)
summary(mod1.slx)

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

Residuals:
   Min     1Q Median     3Q    Max 
-6.511 -3.005 -1.008  2.432 16.870 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.7377     0.7334   5.097 4.91e-06 ***
RD90          2.8635     1.1151   2.568   0.0131 *  
lag.RD90     -0.2870     2.3742  -0.121   0.9043    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.474 on 52 degrees of freedom
Multiple R-squared:  0.1165,    Adjusted R-squared:  0.08248 
F-statistic: 3.427 on 2 and 52 DF,  p-value: 0.03998

summary(spatialreg::impacts(mod1.slx, tr=trMat, R=1000), zstats=TRUE, short=TRUE)

Impact measures (SlX, glht, n-k):
      Direct  Indirect    Total
RD90 2.86351 -0.286971 2.576539
========================================================
Standard errors:
       Direct Indirect   Total
RD90 1.115113 2.374184 2.76964
========================================================
Z-values:
      Direct   Indirect     Total
RD90 2.56791 -0.1208714 0.9302793

p-values:
     Direct   Indirect Total  
RD90 0.010231 0.90379  0.35223

4.14 SLX (STSLS)

mod1.slx_stsls <- spreg(
  formula = esp1, data = montana.shp@data, listw = w1,
  lag.instr = FALSE, Durbin = TRUE, model = "ols", step1.c = TRUE, 
  het = TRUE, endog = endog, instruments = instruments
)

summary(mod1.slx_stsls)

 Stsls with Spatial HAC standard errors

Call:
spreg(formula = esp1, data = montana.shp@data, listw = w1, endog = endog, 
    instruments = instruments, lag.instr = FALSE, model = "ols", 
    het = TRUE, step1.c = TRUE, Durbin = TRUE)

Residuals:
    Min.  1st Qu.   Median  3rd Qu.     Max. 
-6.40823 -2.46165 -0.84138  1.08488 17.06315 

Coefficients:
            Estimate Std. Error t-value  Pr(>|t|)    
(Intercept) -1.05506    1.15580 -0.9128 0.3613275    
RD90         0.23420    1.04646  0.2238 0.8229104    
lag_RD90     0.13117    2.06873  0.0634 0.9494441    
UE90         0.74709    0.19354  3.8602 0.0001133 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mod1.slx_stsls2 <- spreg(
  formula = esp1, data = montana.shp@data, listw = w1,
  lag.instr = FALSE, Durbin = TRUE, model = "ols", step1.c = TRUE, 
  het = TRUE
)

summary(sphet::impacts(mod1.slx_stsls2, tr=trMat, R=1000), zstats=TRUE, short=TRUE)

Impact measures (SlX, glht, n-k):
      Direct  Indirect    Total
RD90 2.86351 -0.286971 2.576539
========================================================
Standard errors:
       Direct Indirect    Total
RD90 1.213777 2.216879 2.652204
========================================================
Z-values:
       Direct   Indirect     Total
RD90 2.359174 -0.1294482 0.9714708

p-values:
     Direct   Indirect Total  
RD90 0.018316 0.897    0.33131

4.15 SDM (MV)

\(y = \rho Wy + X \beta + WX \theta + \varepsilon\)

mod1.sdm <- lagsarlm(formula = esp1, data = montana.shp@data, listw = w1, Durbin = TRUE)
summary(mod1.sdm)

Call:lagsarlm(formula = esp1, data = montana.shp@data, listw = w1, 
    Durbin = TRUE)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.3395 -2.9909 -1.0769  2.5768 16.5460 

Type: mixed 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  2.91618    0.98963  2.9467 0.003212
RD90         2.84391    1.06989  2.6581 0.007857
lag.RD90    -0.99818    2.33320 -0.4278 0.668785

Rho: 0.20387, LR test value: 1.002, p-value: 0.31682
Asymptotic standard error: 0.18781
    z-value: 1.0855, p-value: 0.27771
Wald statistic: 1.1783, p-value: 0.27771

Log likelihood: -158.4075 for mixed model
ML residual variance (sigma squared): 18.422, (sigma: 4.2921)
Number of observations: 55 
Number of parameters estimated: 5 
AIC: 326.82, (AIC for lm: 325.82)
LM test for residual autocorrelation
test value: 1.0407, p-value: 0.30766

summary(spatialreg::impacts(mod1.sdm, tr=trMat, R=1000), zstats=TRUE, short=TRUE)

Impact measures (mixed, trace):
       Direct   Indirect    Total
RD90 2.825151 -0.5067765 2.318375
========================================================
Simulation results ( variance matrix):
========================================================
Simulated standard errors
      Direct Indirect    Total
RD90 1.13206 3.118891 3.679099

Simulated z-values:
       Direct   Indirect     Total
RD90 2.523144 -0.1868899 0.6179394

Simulated p-values:
     Direct   Indirect Total  
RD90 0.011631 0.85175  0.53662

4.16 SDEM (MV)

\(y = X \beta + WX \theta + u \\ u = \lambda Wu + \varepsilon\)

mod1.sdem <- errorsarlm(formula = esp1, data = montana.shp@data, listw = w1, Durbin = TRUE)
summary(mod1.sdem)

Call:errorsarlm(formula = esp1, data = montana.shp@data, listw = w1, 
    Durbin = TRUE)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.3568 -2.9878 -1.0313  2.5668 16.5735 

Type: error 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept)  3.69127    0.85423  4.3212 1.552e-05
RD90         2.83836    1.10454  2.5697   0.01018
lag.RD90    -0.33914    2.45675 -0.1380   0.89021

Lambda: 0.20244, LR test value: 0.98866, p-value: 0.32007
Asymptotic standard error: 0.18799
    z-value: 1.0769, p-value: 0.28153
Wald statistic: 1.1597, p-value: 0.28153

Log likelihood: -158.4142 for error model
ML residual variance (sigma squared): 18.429, (sigma: 4.2929)
Number of observations: 55 
Number of parameters estimated: 5 
AIC: NA (not available for weighted model), (AIC for lm: 325.82)

4.17 SDEM (GMM)

mod1.sdem_gmm <- spreg(
  formula = esp1, data = montana.shp@data, listw = w1,
  model = "error", Durbin = TRUE, 
  step1.c = TRUE, het = TRUE, endog = endog, instruments = instruments
)

summary(mod1.sdem_gmm)

 Generalized stsls

Call:
spreg(formula = esp1, data = montana.shp@data, listw = w1, endog = endog, 
    instruments = instruments, model = "error", het = TRUE, step1.c = TRUE, 
    Durbin = TRUE)

Residuals:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-6.3805 -2.4459 -0.8510  0.0021  1.0907 17.0541 

Coefficients:
             Estimate Std. Error t-value  Pr(>|t|)    
(Intercept) -1.012074   1.163428 -0.8699 0.3843513    
RD90         0.265070   1.061400  0.2497 0.8027915    
lag_RD90     0.092892   2.131567  0.0436 0.9652397    
UE90         0.739299   0.195456  3.7824 0.0001553 ***
rho          0.064293   0.302432  0.2126 0.8316499    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mod1.sdem_gmm2 <- spreg(
  formula = esp1, data = montana.shp@data, listw = w1,
  model = "error", Durbin = TRUE, 
  step1.c = TRUE, het = TRUE
)

summary(mod1.sdem_gmm2)

 Generalized stsls

Call:
spreg(formula = esp1, data = montana.shp@data, listw = w1, model = "error", 
    het = TRUE, step1.c = TRUE, Durbin = TRUE)

Residuals:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-6.4510 -2.9793 -0.9784  0.0323  2.4635 16.8868 

Coefficients:
            Estimate Std. Error t-value Pr(>|t|)    
(Intercept)  3.69503    0.78036  4.7350 2.19e-06 ***
RD90         2.84006    1.26364  2.2475  0.02461 *  
lag_RD90    -0.33439    2.40223 -0.1392  0.88929    
rho          0.19018    0.20037  0.9491  0.34256    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(sphet::impacts(mod1.sdem_gmm2, tr=w_matrix$tr, R=1000), zstats=TRUE, short=TRUE)

Impact measures (SlX, glht, n-k):
       Direct   Indirect    Total
RD90 2.840057 -0.3343859 2.505671
========================================================
Standard errors:
       Direct Indirect   Total
RD90 1.263635 2.402231 3.03313
========================================================
Z-values:
      Direct  Indirect     Total
RD90 2.24753 -0.139198 0.8261007

p-values:
     Direct   Indirect Total  
RD90 0.024606 0.88929  0.40875