Anselin Kelejian Test

[JosiahParry]

After using S2SLS, it is recommended to check for remaining spatial autocorrelation using the Anselin-Kelejian (AK) Test. When the test statistic is significant it suggests that HAC standard errors should be reported to rather than white standard errors.

Is the AK test available yet in spatialreg?

[JosiahParry]

The test is defined here for pysal: https://github.com/pysal/spreg/blob/3ff33bcadbf70f0e1ad61099f0ea8e47966fa03c/spreg/diagnostics_sp.py#L921-L954

[rsbivand]

@JosiahParry My understanding is that the AK test only applies when there are no spatial endogenous variables on the RHS, so non-spatial IV (Anselin & Rey p. 142):

... for the special case where only non-spatial endoogenous variables are included in the model specification (i.e., specifically excluding a spatially lagged dependent variable on the RHS), ...

There are no non-spatial endogenous varables in s2sls, so that AK test is not relevant.

On p. 166-7, the AK_L test looks more like the actual implementation (no equation numbers there), and looks similar to Anselin & Kelejian (1997) p. 163, eq. 11-12. The next question would be whether this might not be better placed in sphet compared to here, since Luc's implementation is almost twenty years old. @gpiras what do you think?

[gpiras]

For what I remember (but I read Anselin and Kelejian a few years ago) the specification includes both a spatial lag and additional endogenous variable. If I recall correctly, they include the separate specification as special cases. As for the suggestion of including the test in sphet I am not really sure it would be a good idea. The user can use the spatial HAC to produce robust standard errors by default, and if they believe that the error follows a spatial autoregressive process specify a SARAR model.

[rsbivand]

@JosiahParry Could you provide a worked example from pysal/spreg for example through reticulate? I'm having trouble with not using q, the example should use wlags=2. It should be possible to calculate AK for special case 1 in s2sls and report it in the summary method.

[JosiahParry]

@rsbivand
here you go!

library(reticulate)

# reticulate::virtualenv_create("spatial", packages = c("numpy", "spreg", "libpysal"))
use_virtualenv("spatial")

lp <- import("libpysal")
np <- import("numpy")
spreg <- import("spreg")

columbus <- lp$examples$get_path("columbus.dbf")
columbus_sf <- sf::st_read(columbus)

# oddly passing in `columbus` here doesn't work! 
w <- lp$weights$contiguity$Rook$from_shapefile(lp$examples$get_path("columbus.shp"))
w$transform = "r"

reg <- spreg$GM_Lag(np$array(columbus_sf$HOVAL), np$array(cbind(columbus_sf$INC, columbus_sf$CRIME)), w = w)

reg$ak_test
#> [[1]]
#> [1] 2.075273
#> 
#> [[2]]
#> [1] 0.1497032

[rsbivand]

reticulate couldn't find my local updated python packages. I have:

import numpy as np
import libpysal
from spreg import GM_Lag, AKtest
db = libpysal.io.open(libpysal.examples.get_path("columbus.dbf"),'r')
y = np.array(db.by_col("CRIME"))
y = np.reshape(y, (49,1))
X = []
X.append(db.by_col("INC"))
X.append(db.by_col("HOVAL"))
X = np.array(X).T
w = libpysal.weights.Rook.from_shapefile(libpysal.examples.get_path("columb\
us.shp"))
w.transform = 'r'
reg_lag = GM_Lag(y, X, w=w)
ak_sp = AKtest(reg_lag, w, case='gen')
print('AK test: %f\tP-value: %f'%(ak_sp.ak, ak_sp.p))
# AK test: 0.147798	P-value: 0.700648

The case='gen' seems to matter; https://github.com/pysal/spreg/blob/3ff33bcadbf70f0e1ad61099f0ea8e47966fa03c/spreg/diagnostics_sp.py#L303-L307.

[JosiahParry]

Ah, sorry about the python environment issue :/ classic problem.

Good eye! So it seems there is an AK test for the SLM and SARAR scenarios respectively. I was unaware!

[rsbivand]

Well, I think that 'nosp' is for the case where no \rho y is included but some regressors are endogenous and instrumented, and 'gen' is for the case when \rho y is present and some regressors are endogenous and instrumented. How either perform in our setting where no regressors are endogenous and instrumented is a bit open - s2sls does not admit endogenous regressors. Also:

reg_lag1 = GM_Lag(y, X, w=w, w_lags=1, sig2n_k=True)
reg_lag1.betas.T
# array([[45.94746499, -1.05273965, -0.25984195,  0.4104519 ]])
se_betas(reg_lag1)
# array([11.51533771,  0.39427781,  0.09253382,  0.19386571])

matches:

library(spatialreg)
col <- st_read("/home/rsb/.local/lib/python3.13/site-packages/libpysal/examples/columbus/columbus.shp")
nb <- spdep::poly2nb(col, queen=FALSE)
lw <- spdep::nb2listw(nb)
lag.stsls1 <- stsls(CRIME ~ INC + HOVAL, data=col, lw, W2X=FALSE)
coef(lag.stsls1)
#        Rho (Intercept)         INC       HOVAL 
#  0.4104519  45.9474650  -1.0527396  -0.2598419 
sqrt(diag(lag.stsls1$var))
# [1]  0.19386571 11.51533771  0.39427781  0.09253382

suggesting that stsls needs a sig2n_k argument - it is default False in GM_Lag.

[rsbivand]

@JosiahParry Please try from the AK branch if you like - the test value is calculated un stsls and reported by printing the summary, commit 3faf1fe.

> library(spatialreg)
Loading required package: spData
Loading required package: Matrix
Loading required package: sf
Linking to GEOS 3.13.0, GDAL 3.10.0, PROJ 9.5.0; sf_use_s2() is TRUE
> col <- st_read("/home/rsb/.local/lib/python3.13/site-packages/libpysal/examples/columbus/columbus.shp")
Reading layer `columbus' from data source 
  `/home/rsb/.local/lib/python3.13/site-packages/libpysal/examples/columbus/columbus.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 49 features and 20 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 5.874907 ymin: 10.78863 xmax: 11.28742 ymax: 14.74245
CRS:           NA
> nb <- spdep::poly2nb(col, queen=FALSE)
> lw <- spdep::nb2listw(nb)
> lag.stsls2 <- stsls(CRIME ~ INC + HOVAL, data=col, lw)
> summary(lag.stsls2)

Call:stsls(formula = CRIME ~ INC + HOVAL, data = col, listw = lw)

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)
Anselin-Kelejian (1997) test for residial autocorrelation: 0.12991
        p-value: 0.71852

$ python
Python 3.13.0 (main, Oct  8 2024, 00:00:00) [GCC 14.2.1 20240912 (Red Hat 14.2.1-3)] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy as np
>>> import libpysal
>>> from spreg import GM_Lag, AKtest
>>> db = libpysal.io.open(libpysal.examples.get_path("columbus.dbf"),'r')
>>> y = np.array(db.by_col("CRIME"))
>>> y = np.reshape(y, (49,1))
>>> X = []
>>> X.append(db.by_col("INC"))
>>> X.append(db.by_col("HOVAL"))
>>> X = np.array(X).T
>>> w = libpysal.weights.Rook.from_shapefile(libpysal.examples.get_path("columbus.shp"))
/home/rsb/.local/lib/python3.13/site-packages/libpysal/io/iohandlers/pyShpIO.py:247: FutureWarning: Objects based on the `Geometry` class will deprecated and removed in a future version of libpysal.
  shp = self.type(vertices)
/home/rsb/.local/lib/python3.13/site-packages/libpysal/cg/shapes.py:1408: FutureWarning: Objects based on the `Geometry` class will deprecated and removed in a future version of libpysal.
  self._part_rings = [Ring(vertices)]
>>> w.transform = 'r'
>>> reg_lag = GM_Lag(y, X, w=w, w_lags=2, sig2n_k=False)
>>> ak_sp = AKtest(reg_lag, w, case='gen')
>>> print('AK test: %f\tP-value: %f'%(ak_sp.ak, ak_sp.p))
AK test: 0.129912	P-value: 0.718523
>>> reg_lag.betas.T
array([[45.45909249, -1.0410089 , -0.25953844,  0.41929355]])

[JosiahParry]

@rsbivand this looks great!

library(spatialreg)

ncfp <- system.file("shape/nc.shp", package = "sf")
nc <- sf::st_read(ncfp)

listw <- spdep::poly2nb(nc) |> 
  spdep::nb2listw()

summary(stsls(SID74 ~ NWBIR74, nc, listw))
#> Call:stsls(formula = SID74 ~ NWBIR74, data = nc, listw = listw)
#> 
#> Residuals:
#>       Min        1Q    Median        3Q       Max 
#> -10.63051  -2.04121  -0.53359   1.65751  16.71110 
#> 
#> Coefficients: 
#>               Estimate Std. Error t value Pr(>|t|)
#> Rho         0.01515452 0.09207607  0.1646  0.86927
#> (Intercept) 1.45620634 0.70064925  2.0784  0.03768
#> NWBIR74     0.00486227 0.00025487 19.0776  < 2e-16
#> 
#> Residual variance (sigma squared): 11.98, (sigma: 3.4612)
#> Anselin-Kelejian (1997) test for residial autocorrelation: 0.032647
#>         p-value: 0.85661

tmp <- tempfile(fileext = ".gal")
spdep::write.nb.gal(listw$neighbours, tmp)

library(reticulate)

use_virtualenv("spatial")
lp <- import("libpysal")
np <- import("numpy")
spreg <- import("spreg")

w <- lp$io$open(tmp)$read()
w$transform = "r"

reg <- spreg$GM_Lag(np$array(nc$SID74), np$array(cbind(nc$NWBIR74)), w = w, w_lags=2L, sig2n_k=FALSE)
reg$ak_test
#> [[1]]
#> [1] 0.03264722
#> 
#> [[2]]
#> [1] 0.8566145

[rsbivand]

Good, I'll create PR from the branch to merge to main later today, and see how CI goes there. Do you need this released soon?

[JosiahParry]

I'm in no hurry! This is sufficient for me to work against. Thank you very much :)

[rsbivand]

PR merged into main. Thanks for the suggestion. I'll add a news item and better examples.

[JosiahParry]

Brilliant!