`statspai.spatial`¶

spatial ¶

Spatial econometrics — StatsPAI's answer to R's spatialreg + spdep + mgwr and Python's PySAL.

Three layers, flat-imported for convenience:

Weights

W, queen_weights, rook_weights, knn_weights, distance_band, kernel_weights, block_weights.

Exploratory spatial data analysis (ESDA)

moran, moran_local, geary, getis_ord_g, getis_ord_local, join_counts, moran_plot, lisa_cluster_map.

Regression

ML estimators: sar, sem, sdm, slx, sac
Diagnostics: lm_tests, moran_residuals
Effects: impacts (LeSage-Pace direct / indirect / total)

All regression estimators accept a :class:W object, a scipy.sparse matrix, or an (n, n) ndarray for the weights argument.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 100
>>> coords = rng.uniform(0, 10, size=(n, 2))
>>> df = pd.DataFrame({"x1": rng.normal(size=n), "x2": rng.normal(size=n)})
>>> df["y"] = 1 + 0.5 * df["x1"] - 0.3 * df["x2"] + rng.normal(size=n)
>>> w = sp.knn_weights(coords, k=6); w.transform = "R"
>>> moran = sp.moran(df["y"], w, permutations=99, seed=0)
>>> result = sp.sar(w, data=df, formula='y ~ x1 + x2')
>>> list(result.params.index)
['const', 'x1', 'x2', 'rho']
>>> eff = sp.impacts(result, n_sim=200, seed=0)
>>> list(eff.columns)
['Direct', 'SE_Direct', 'Indirect', 'SE_Indirect', 'Total', 'SE_Total']
>>> lms = sp.lm_tests("y ~ x1 + x2", df, w)
>>> sorted(lms)
['LM_err', 'LM_lag', 'Robust_LM_err', 'Robust_LM_lag', 'SARMA']

SpatialDiDResult `dataclass` ¶

Bases: ResultProtocolMixin

Result object for :func:spatial_did.

The object intentionally exposes both a spatial-DiD-specific surface (direct_effect, spillover_effect, total_effect) and a broom/modelsummary-compatible surface (tidy(), glance(), params, std_errors).

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n_units, n_time = 12, 6
>>> W = np.zeros((n_units, n_units))  # ring spatial weights over units
>>> for i in range(n_units):
...     W[i, (i + 1) % n_units] = 1
...     W[i, (i - 1) % n_units] = 1
>>> rows = []
>>> for u in range(n_units):
...     fe = rng.normal()
...     for t in range(n_time):
...         d = 1 if (u < 6 and t >= 3) else 0  # half treated from t>=3
...         y = 1.0 + fe + 0.3 * t + 1.5 * d + rng.normal(scale=0.3)
...         rows.append({"unit": u, "time": t, "treat": d, "y": y})
>>> df = pd.DataFrame(rows)
>>> res = sp.spatial_did(df, y="y", treat="treat", unit="unit",
...                     time="time", W=W)
>>> isinstance(res, sp.SpatialDiDResult)
True
>>> res.n_obs
72

plot ¶

plot(kind: str = 'coef', ax: Any = None, figsize: tuple[int, int] = (8, 5), **kwargs: Any) -> Any

Plot coefficient, exposure, or spatial event-study diagnostics.

SpatialIVResult `dataclass` ¶

Bases: ResultProtocolMixin

Container for :func:spatial_iv Kelejian-Prucha S2SLS estimates.

Attributes:

Name	Type	Description
`rho`	`float`	Spatial autoregressive coefficient on `WY` (`nan` when `include_WY=False`).
`rho_se`	`float`	Standard error of `rho`.
`coefficients`	`DataFrame`	Full coefficient table with `variable`, `coef`, `se` (HC0), `z`, `p` (two-sided normal) and the `1 - alpha` normal confidence bounds `ci_lower` / `ci_upper`.
`n_obs`	`int`	Number of observations after dropping missing rows.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 60
>>> W = np.zeros((n, n))  # ring network, row-normalised
>>> for i in range(n):
...     W[i, (i + 1) % n] = 1
...     W[i, (i - 1) % n] = 1
>>> W = W / W.sum(1, keepdims=True)
>>> x = rng.normal(size=n)
>>> z = rng.normal(size=n)  # excluded instrument
>>> d = 0.5 * x + 0.8 * z + rng.normal(size=n)  # endogenous regressor
>>> y = 1.0 + 1.5 * d + 0.3 * x + rng.normal(size=n)
>>> df = pd.DataFrame({"y": y, "d": d, "x": x, "z": z})
>>> res = sp.spatial_iv(df, y="y", endog=["d"], exog=["x"], W=W,
...                     instruments=["z"])
>>> isinstance(res, sp.SpatialIVResult)
True
>>> res.n_obs
60

SpatialModel ¶

Unified spatial regression estimator.

Uses concentrated ML: profile out β and σ², optimise over ρ (or λ) in one dimension using eigenvalues of W.

The engine behind :func:sp.sar / :func:sp.sem / :func:sp.sdm. Construct it with a dense (n, n) weights matrix W, the data and a formula, then call :meth:fit to obtain an EconometricResults carrying the spatial parameter (rho for SAR / SDM, lambda for SEM).

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 50
>>> W = np.zeros((n, n))            # nearest-neighbour chain weights
>>> for i in range(n - 1):
...     W[i, i + 1] = W[i + 1, i] = 1.0
>>> df = pd.DataFrame({"x": rng.normal(size=n)})
>>> df["y"] = 1.0 + 0.5 * df["x"] + rng.normal(size=n)
>>> model = sp.SpatialModel(W, df, "y ~ x", model_type="sar")
>>> res = model.fit()
>>> type(res).__name__
'EconometricResults'
>>> model.model_type
'sar'

fit ¶

fit() -> EconometricResults

Fit by maximum likelihood with :func:sp.sar / :func:sp.sem / :func:sp.sdm, so the class and the functions report one set of estimates and standard errors.

Notes

Through 1.40.2 the class kept its own concentrated-likelihood fit with the information-matrix formulas the functions had replaced: the rho SE omitted tr(G'G), the beta SEs treated rho as known, the lambda SE used tr(G^2) + tr(G)^2 / n, and alpha never reached the intervals.

W ¶

Spatial weights matrix with CSR-sparse backing.

Parameters:

Name	Type	Description	Default
`neighbors`	`dict[int, list[int]]`	Mapping observation id -> list of neighbour ids.	required
`weights`	`dict[int, list[float]]`	Matching weight values. If `None`, binary (1.0) weights are used.	`None`
`id_order`	`sequence`	Explicit ordering of observation ids. Defaults to `sorted(neighbors)`.	`None`

Notes

w.transform re-weights the weights as constructed (never a previous transform): "O" original, "B" binary, "R" row standardised (spdep style "W"), "V" variance stabilising (spdep "S": w_ij / sqrt(sum_j w_ij^2) rescaled to sum n), "D" globally standardised to sum 1 (spdep "U"). Each matches spdep::nb2listw(style=, glist=) and libpysal's transform.

Examples:

>>> import statspai as sp
>>> # A 4-node chain graph 0-1-2-3
>>> neighbors = {0: [1], 1: [0, 2], 2: [1, 3], 3: [2]}
>>> w = sp.W(neighbors)
>>> w.n
4
>>> w.islands
[]
>>> w.transform = "R"          # row-standardize the weights
>>> float(w.full()[1].sum())   # each row now sums to 1
1.0

spatial_did ¶

spatial_did(data: DataFrame, y: str, treat: str, unit: str, time: str, W: Any, covariates: Optional[Sequence[str]] = None, cluster: Optional[str] = None, alpha: float = 0.05, *, se_type: str = 'cluster', coords: Optional[Tuple[str, str]] = None, lat: Optional[str] = None, lon: Optional[str] = None, distance_matrix: Optional[ndarray] = None, conley_cutoff: Optional[float] = None, conley_kernel: str = 'bartlett', unit_order: Optional[Sequence[Any]] = None, normalize_W: bool = True, event_study: bool = False, event_window: Optional[Tuple[int, int]] = None, event_base: int = -1) -> SpatialDiDResult

Spatial DiD with direct and spillover treatment effects.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`	Long panel with one row per unit-period.	required
`y`	`str`	Outcome, treatment, unit-id, and time-id columns. `treat` may be binary or a continuous exposure intensity; the spatial lag uses the same scale.	required
`treat`	`str`	Outcome, treatment, unit-id, and time-id columns. `treat` may be binary or a continuous exposure intensity; the spatial lag uses the same scale.	required
`unit`	`str`	Outcome, treatment, unit-id, and time-id columns. `treat` may be binary or a continuous exposure intensity; the spatial lag uses the same scale.	required
`time`	`str`	Outcome, treatment, unit-id, and time-id columns. `treat` may be binary or a continuous exposure intensity; the spatial lag uses the same scale.	required
`W`	ndarray, sparse matrix, or :class:`statspai.spatial.W`	Spatial weights matrix over units. If a StatsPAI `W` object carries an `id_order`, that order is respected.	required
`covariates`	`sequence of str`	Additional controls included after absorbing unit and time effects.	`None`
`cluster`	`str`	Cluster variable for `se_type='cluster'`. Defaults to `unit`.	`None`
`alpha`	`float`	Significance level for confidence intervals.	`0.05`
`se_type`	`('cluster', 'robust', 'conley')`	Covariance estimator. `conley` uses spatial-HAC dependence within each period and requires `conley_cutoff` plus coordinates or a unit-level distance matrix.	`"cluster"`
`coords`	`tuple(str, str)`	Convenience pair `(lat, lon)` in degrees.	`None`
`lat`	`str`	Latitude/longitude columns in degrees.	`None`
`lon`	`str`	Latitude/longitude columns in degrees.	`None`
`distance_matrix`	`ndarray`	Unit-by-unit distances in the same units as `conley_cutoff`.	`None`
`conley_cutoff`	`float`	Spatial cutoff for Conley-HAC standard errors.	`None`
`conley_kernel`	`('uniform', 'bartlett')`		`"uniform"`
`unit_order`	`sequence`	Explicit row/column order of `W`.	`None`
`normalize_W`	`bool`	Row-normalise weights before constructing `WD`.	`True`
`event_study`	`bool`	Estimate direct and spillover event-study paths by relative time.	`False`
`event_window`	`tuple(int, int)`	Inclusive event-time window. Defaults to observed support capped at [-5, 5] when `event_study=True`.	`None`
`event_base`	`int`	Omitted relative-time period.	`-1`

Returns:

Type	Description
`SpatialDiDResult`	Direct, spillover, total effects with diagnostics, plotting, and export helpers.

References

delgado2015difference

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n_units, n_time = 12, 6
>>> W = np.zeros((n_units, n_units))  # ring spatial weights over units
>>> for i in range(n_units):
...     W[i, (i + 1) % n_units] = 1
...     W[i, (i - 1) % n_units] = 1
>>> rows = []
>>> for u in range(n_units):
...     fe = rng.normal()
...     for t in range(n_time):
...         d = 1 if (u < 6 and t >= 3) else 0  # half treated from t>=3
...         y = 1.0 + fe + 0.3 * t + 1.5 * d + rng.normal(scale=0.3)
...         rows.append({"unit": u, "time": t, "treat": d, "y": y})
>>> df = pd.DataFrame(rows)
>>> res = sp.spatial_did(df, y="y", treat="treat", unit="unit",
...                     time="time", W=W)
>>> res.coefficients["variable"].tolist()
['treat', 'W_treat']
>>> res.n_obs
72

getis_ord_g ¶

getis_ord_g(y: Any, w: W, permutations: int = 999, seed: Optional[int] = None) -> SpatialStatistic

Global Getis-Ord G — overall concentration of high (or low) values.

The global G statistic summarises whether high values of y are globally concentrated among neighbouring units. It requires non-negative y. Larger-than-expected G indicates clustering of high values; smaller indicates clustering of low values.

Parameters:

Name	Type	Description	Default
`y`	`array_like`	Non-negative variable of interest, length `n`.	required
`w`	`W`	Spatial weights object defining the neighbourhood structure.	required
`permutations`	`int`	Number of random permutations for the pseudo p-value, empirical variance and z-score. `0` skips the permutation test.	`999`
`seed`	`int`	Seed for the permutation RNG for reproducibility.	`None`

Returns:

Type	Description
`SpatialStatistic`	Result with `value` (G), `expectation` (`S0 / (n(n-1))`), `variance`, `z_score`, `p_norm` and `p_sim` / `simulations`.

Raises:

Type	Description
`ValueError`	If any element of `y` is negative, or `y` is degenerate.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> coords = np.random.default_rng(0).random((50, 2))
>>> w = sp.knn_weights(coords, k=5)
>>> y = np.random.default_rng(4).random(50)  # non-negative
>>> res = sp.getis_ord_g(y, w, seed=0)
>>> res.value > 0
True

See Also

getis_ord_local : Local Gi / Gi* hotspot statistics.

References

Getis, A. and Ord, J. K. (1992). The analysis of spatial association by use of distance statistics. Geographical Analysis, 24(3), 189-206. [@getis1992analysis]

getis_ord_local ¶

getis_ord_local(y: Any, w: W, star: bool = True, permutations: int = 999, seed: Optional[int] = None) -> Dict[str, ndarray]

Local Getis-Ord Gi / Gi* — hotspot and coldspot detection.

Computes a standardised local statistic for each unit measuring whether it is surrounded by high values (a hotspot, large positive z) or low values (a coldspot, large negative z). With star=True the focal unit is included in its own neighbourhood (Gi*); with star=False it is excluded (Gi).

Parameters:

Name	Type	Description	Default
`y`	`array_like`	Variable of interest, length `n`.	required
`w`	`W`	Spatial weights object; densified internally to an `n x n` matrix.	required
`star`	`bool`	If True compute Gi* (self-included); otherwise Gi (self-excluded).	`True`
`permutations`	`int`	Accepted for API symmetry; the current implementation returns the analytic standardised statistic.	`999`
`seed`	`int`	Seed placeholder for API symmetry.	`None`

Returns:

Type	Description
`dict`	`{"Gs": ndarray, "z": ndarray}` — the `n` local statistics, which are already standardised (so `Gs` and `z` coincide).

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> coords = np.random.default_rng(0).random((40, 2))
>>> w = sp.knn_weights(coords, k=4)
>>> y = np.random.default_rng(5).random(40)
>>> out = sp.getis_ord_local(y, w, star=True)
>>> out["Gs"].shape
(40,)

See Also

getis_ord_g : Global G statistic. moran_local : LISA, a cross-product local statistic.

References

Getis, A. and Ord, J. K. (1992). The analysis of spatial association by use of distance statistics. Geographical Analysis, 24(3), 189-206. [@getis1992analysis]

moran_local ¶

moran_local(y: Any, w: W, permutations: int = 999, seed: Optional[int] = None) -> Dict[str, Any]

Local Moran's I (LISA) — per-location spatial association.

Decomposes global Moran's I into a contribution I_i for each spatial unit, identifying local clusters (high-high / low-low) and spatial outliers (high-low / low-high). Each I_i gets a conditional permutation pseudo p-value.

Parameters:

Name	Type	Description	Default
`y`	`array_like`	Variable of interest, length `n`.	required
`w`	`W`	Spatial weights object; `w.sparse` is the `n x n` weights matrix.	required
`permutations`	`int`	Number of conditional permutations for the per-location pseudo p-values. `0` skips the permutation test.	`999`
`seed`	`int`	Seed for the permutation RNG for reproducibility.	`None`

Returns:

Type	Description
`dict`	`{"Is": ndarray, "p_sim": ndarray \| None, "simulations": ndarray \| None}` where `Is` holds the `n` local statistics and `p_sim` the matching pseudo p-values.

Raises:

Type	Description
`ValueError`	If `y` has zero variance.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> coords = np.random.default_rng(0).random((40, 2))
>>> w = sp.knn_weights(coords, k=4)
>>> y = np.random.default_rng(2).standard_normal(40)
>>> out = sp.moran_local(y, w, permutations=499, seed=0)
>>> out["Is"].shape
(40,)

See Also

moran : Global Moran's I.

References

Anselin, L. (1995). Local indicators of spatial association---LISA. Geographical Analysis, 27(2), 93-115. [@anselin1995local]

lisa_cluster_map ¶

lisa_cluster_map(y: Any, w: W, gdf: Any, ax: Optional[Any] = None, p_threshold: float = 0.05) -> Any

Classify each observation HH/LL/HL/LH/Not significant and colour on a GDF.

Requires geopandas (install with pip install geopandas); the gdf must align row-for-row with y and w.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> coords = rng.uniform(size=(60, 2))
>>> w = sp.knn_weights(coords, k=5)
>>> y = coords[:, 0] + rng.normal(scale=0.3, size=60)
>>> import geopandas as gpd
>>> from shapely.geometry import Point
>>> pts = [Point(xy) for xy in coords]
>>> gdf = gpd.GeoDataFrame(geometry=pts)
>>> ax = sp.lisa_cluster_map(y, w, gdf)
>>> ax.figure.savefig("lisa.png")

moran_plot ¶

moran_plot(y: Any, w: W, ax: Optional[Any] = None) -> Any

Moran scatter: z vs spatial lag Wz. Slope of OLS line = Moran's I.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> coords = rng.uniform(size=(60, 2))
>>> w = sp.knn_weights(coords, k=5)
>>> y = coords[:, 0] + rng.normal(scale=0.3, size=60)
>>> ax = sp.moran_plot(y, w)
>>> ax.figure.savefig("moran.png")

gwr_bandwidth ¶

gwr_bandwidth(coords: Any, y: Any, X: Any, kernel: KernelName = 'bisquare', fixed: bool = False, criterion: Criterion = 'AICc', bw_min: Optional[float] = None, bw_max: Optional[float] = None, add_constant: bool = True, tol: float = 0.0001) -> float

Select a GWR bandwidth by golden-section search.

Minimises criterion of :func:statspai.spatial.gwr.gwr over the bandwidth, as R GWmodel::bw.gwr does.

Parameters:

Name	Type	Description	Default
`coords`	`(n, 2) array-like`	Projected point coordinates.	required
`y`	`(n,) array-like`	Outcome.	required
`X`	`(n, p) array-like`	Regressors (without a constant unless `add_constant=False`).	required
`kernel`	`('bisquare', 'gaussian', 'exponential')`	Kernel passed to :func:`gwr`.	`'bisquare'`
`fixed`	`bool`	`False`: adaptive bandwidth, a number of nearest neighbours (searched over integers). `True`: a fixed distance.	`False`
`criterion`	`('AICc', 'AIC', 'BIC', 'CV')`	Score to minimise; `'CV'` is the leave-one-out residual sum of squares. `GWmodel` offers `'AICc'` and `'CV'`.	`'AICc'`
`bw_min`	`float`	Search interval. Defaults: `[20, n]` neighbours (adaptive) or `[bw_max / 5000, bw_max]` with `bw_max` defaulting to the largest distance between two points (fixed), the reference's interval. With `n <= 20` observations the adaptive interval `[20, n]` would be empty and starts at `k + 2` instead (`k` local coefficients, the constant included). Adaptive searches can end one or two neighbours outside the interval, as the reference's do, because their points are rounded.	`None`
`bw_max`	`float`	Search interval. Defaults: `[20, n]` neighbours (adaptive) or `[bw_max / 5000, bw_max]` with `bw_max` defaulting to the largest distance between two points (fixed), the reference's interval. With `n <= 20` observations the adaptive interval `[20, n]` would be empty and starts at `k + 2` instead (`k` local coefficients, the constant included). Adaptive searches can end one or two neighbours outside the interval, as the reference's do, because their points are rounded.	`None`
`add_constant`	`bool`	Passed to :func:`gwr`.	`True`
`tol`	`float`	The search stops when two interior criterion values differ by at most `tol` or the golden step falls to `tol`.	`1e-4`

Returns:

Type	Description
`float`	The selected bandwidth.

Notes

A probed bandwidth that admits no GWR fit (one neighbour or less) scores +inf, with a RuntimeWarning; if no probe can be fitted the first error is raised.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> coords = rng.uniform(0, 10, size=(60, 2))
>>> X = rng.normal(size=(60, 1))
>>> y = 1 + (0.5 + 0.1 * coords[:, 0]) * X[:, 0] + rng.normal(0, 0.3, 60)
>>> bw = sp.gwr_bandwidth(coords, y, X)
>>> bool(bw > 0)
True

References

brunsdon1996geographically; gollini2015gwmodel; kiefer1953sequential.

spatial_iv ¶

spatial_iv(data: DataFrame, y: str, endog: Sequence[str], exog: Sequence[str], W: Any, instruments: Optional[Sequence[str]] = None, include_WY: bool = True, alpha: float = 0.05) -> SpatialIVResult

Spatial 2SLS estimator.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`		required
`y`	`str`		required
`endog`	`sequence of str`	Names of additional endogenous regressors.	required
`exog`	`sequence of str`	Names of exogenous regressors (included instruments).	required
`W`	ndarray or :class:`spatial.W`	Spatial weights matrix (n_obs x n_obs). Should be row-normalised.	required
`instruments`	`sequence of str`	Extra excluded instruments beyond the spatial lags of X.	`None`
`include_WY`	`bool`	If True, include WY as a regressor (spatial autoregressive coefficient ρ).	`True`
`alpha`	`float`	Significance level of the normal confidence intervals in `coefficients` (`ci_lower` / `ci_upper`).	`0.05`

Returns:

Type	Description
`SpatialIVResult`

References

kelejian1998generalized

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 60
>>> W = np.zeros((n, n))  # ring network, row-normalised
>>> for i in range(n):
...     W[i, (i + 1) % n] = 1
...     W[i, (i - 1) % n] = 1
>>> W = W / W.sum(1, keepdims=True)
>>> x = rng.normal(size=n)
>>> z = rng.normal(size=n)  # excluded instrument for the endogenous d
>>> d = 0.5 * x + 0.8 * z + rng.normal(size=n)
>>> y = 1.0 + 1.5 * d + 0.3 * x + rng.normal(size=n)
>>> df = pd.DataFrame({"y": y, "d": d, "x": x, "z": z})
>>> res = sp.spatial_iv(df, y="y", endog=["d"], exog=["x"], W=W,
...                     instruments=["z"])
>>> res.coefficients["variable"].tolist()
['d', 'rho (WY)', '(Intercept)', 'x']
>>> res.n_obs
60

sar ¶

sar(W: ArrayOrW, data: DataFrame, formula: str, row_normalize: bool = True, alpha: float = 0.05, weights: Optional[Any] = None) -> EconometricResults

Spatial Autoregressive (Lag) Model: Y = ρWY + Xβ + ε.

Parameters:

Name	Type	Description	Default
`W`	`ndarray, scipy.sparse matrix, or statspai.spatial.weights.W`	n x n spatial weights matrix, rows in the order of `data`.	required
`data`	`DataFrame`	Data with the variables named in `formula`.	required
`formula`	`str`	Model formula `"y ~ x1 + x2"`; a constant is added.	required
`row_normalize`	`bool`	Divide each row of `W` by its sum before fitting (rows that sum to zero are left as they are).	`True`
`weights`	`None`	Not supported, and refused with `MethodIncompatibility` rather than ignored: the likelihood is a joint one over all areas linked by `W`, not a sum over independently sampled units, so sampling weights do not define a weighted version of it.	`None`

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 30
>>> W = np.zeros((n, n))  # ring contiguity: each unit linked to neighbours
>>> for i in range(n):
...     W[i, (i - 1) % n] = 1.0
...     W[i, (i + 1) % n] = 1.0
>>> x = rng.normal(size=n)
>>> df = pd.DataFrame({"y": 1.0 + 0.8 * x + rng.normal(scale=0.5, size=n),
...                    "x": x})
>>> res = sp.sar(W, data=df, formula="y ~ x")
>>> list(res.params.index)
['const', 'x', 'rho']

References

anselin1988spatial

sdm ¶

sdm(W: ArrayOrW, data: DataFrame, formula: str, row_normalize: bool = True, alpha: float = 0.05) -> EconometricResults

Spatial Durbin Model: Y = ρWY + Xβ + WXθ + ε.

Reports LeSage & Pace direct / indirect / total effects in the diagnostics when n is small enough for the dense (I - ρW)^{-1}.

Parameters:

Name	Type	Description	Default
`W`	`ndarray, scipy.sparse matrix, or statspai.spatial.weights.W`	n x n spatial weights matrix, rows in the order of `data`.	required
`data`	`DataFrame`	Data with the variables named in `formula`.	required
`formula`	`str`	Model formula `"y ~ x1 + x2"`; a constant is added.	required
`row_normalize`	`bool`	Divide each row of `W` by its sum before fitting (rows that sum to zero are left as they are).	`True`

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 30
>>> W = np.zeros((n, n))
>>> for i in range(n):
...     W[i, (i - 1) % n] = 1.0
...     W[i, (i + 1) % n] = 1.0
>>> x = rng.normal(size=n)
>>> df = pd.DataFrame({"y": 1.0 + 0.8 * x + rng.normal(scale=0.5, size=n),
...                    "x": x})
>>> res = sp.sdm(W, data=df, formula="y ~ x")
>>> list(res.params.index)  # own + spatially-lagged covariate + rho
['const', 'x', 'W_x', 'rho']

References

lesage2009introduction

sem ¶

sem(W: ArrayOrW, data: DataFrame, formula: str, row_normalize: bool = True, alpha: float = 0.05) -> EconometricResults

Spatial Error Model: Y = Xβ + u, u = λWu + ε.

Parameters:

Name	Type	Description	Default
`W`	`ndarray, scipy.sparse matrix, or statspai.spatial.weights.W`	n x n spatial weights matrix, rows in the order of `data`.	required
`data`	`DataFrame`	Data with the variables named in `formula`.	required
`formula`	`str`	Model formula `"y ~ x1 + x2"`; a constant is added.	required
`row_normalize`	`bool`	Divide each row of `W` by its sum before fitting (rows that sum to zero are left as they are).	`True`

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 30
>>> W = np.zeros((n, n))
>>> for i in range(n):
...     W[i, (i - 1) % n] = 1.0
...     W[i, (i + 1) % n] = 1.0
>>> x = rng.normal(size=n)
>>> df = pd.DataFrame({"y": 1.0 + 0.8 * x + rng.normal(scale=0.5, size=n),
...                    "x": x})
>>> res = sp.sem(W, data=df, formula="y ~ x")
>>> list(res.params.index)
['const', 'x', 'lambda']

References

anselin1988spatial

lm_tests ¶

lm_tests(formula: str, data: DataFrame, W: Any, row_normalize: bool = True) -> Dict[str, Tuple[float, float]]

Anselin (1988) Lagrange multiplier battery.

Parameters:

Name	Type	Description	Default
`formula`	`str`	`"y ~ x1 + x2"` style formula (constant added automatically).	required
`data`	`DataFrame`		required
`W`	ndarray, scipy.sparse, or ``statspai.spatial.weights.W``		required
`row_normalize`	`bool`		`True`

Returns:

Type	Description
`dict`	`{name: (statistic, p_value)}` for the five tests.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> import pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 80
>>> coords = rng.uniform(size=(n, 2))
>>> w = sp.knn_weights(coords, k=5)
>>> x1 = rng.normal(size=n)
>>> y = 1.0 + 0.5 * x1 + rng.normal(scale=0.5, size=n)
>>> df = pd.DataFrame({"y": y, "x1": x1})
>>> out = sp.lm_tests("y ~ x1", df, w)
>>> sorted(out.keys())
['LM_err', 'LM_lag', 'Robust_LM_err', 'Robust_LM_lag', 'SARMA']

moran_residuals ¶

moran_residuals(residuals: ndarray, W: Any, row_normalize: bool = True, X: Optional[ndarray] = None) -> Tuple[float, float]

Moran's I applied to regression residuals (quick LM-err companion).

Parameters:

Name	Type	Description	Default
`residuals`	`ndarray`		required
`W`	ndarray, scipy.sparse, or ``statspai.spatial.weights.W``		required
`row_normalize`	`bool`		`True`
`X`	`ndarray`	The regression's design matrix, including its constant. When supplied, the p-value uses Cliff and Ord's null for regression residuals -- the one `spdep::lm.morantest` uses -- which depends on `X` through the hat matrix. Without it the function falls back to the raw-variable null, which is not the residuals' null. .. versionchanged:: 1.27.0 The statistic was always exact (5e-16 against `lm.morantest`). The p-value was computed from the observed-variable null, which is the wrong reference distribution for a projection of `y`; on a 10x10 lattice it read 3.64e-06 where `lm.morantest` reports 1.56e-06 (two-sided 3.13e-06). Pass `X` for the correct null.	`None`

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> coords = rng.uniform(size=(80, 2))
>>> w = sp.knn_weights(coords, k=5)
>>> resid = rng.normal(size=80)
>>> I, p = sp.moran_residuals(resid, w)
>>> bool(np.isfinite(I) and np.isfinite(p))
True

gns ¶

gns(W: ArrayOrW, data: DataFrame, formula: str, row_normalize: bool = True, alpha: float = 0.05) -> EconometricResults

General nesting spatial model (Manski-type): spatial lag of y, of the regressors, and of the error.

y = rho W y + Xb + WXt + u, u = lambda W u + e by maximum likelihood, OIM standard errors as Stata spregress, ml dvarlag(W) errorlag(W) ivarlag(W: x...). The model nests SAR, SEM, SDM, SDEM and SAC; it is weakly identified in small samples (rho and lambda trade off), so read the standard errors of the two before the point estimates.

Parameters:

Name	Type	Default
`W`	`ndarray, scipy.sparse matrix, or statspai.spatial.weights.W`	required
`data`	`DataFrame`	required
`formula`	`str`	required
`row_normalize`	`bool`	`True`
`alpha`	`float`	`0.05`

Returns:

Type	Description
`EconometricResults`	`params` `[const, x..., W_x..., rho, lambda]`.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 40
>>> W = np.zeros((n, n))
>>> for i in range(n):
...     W[i, (i - 1) % n] = W[i, (i + 1) % n] = 1.0
>>> x = rng.normal(size=n)
>>> df = pd.DataFrame({"y": 1 + x + rng.normal(size=n), "x": x})
>>> list(sp.gns(W, df, "y ~ x").params.index)
['const', 'x', 'W_x', 'rho', 'lambda']

sdem ¶

sdem(W: ArrayOrW, data: DataFrame, formula: str, row_normalize: bool = True, alpha: float = 0.05) -> EconometricResults

Spatial Durbin error model: y = Xb + WXt + u, u = lambda W u + e.

Parameters:

Name	Type	Description	Default
`W`	`ndarray, scipy.sparse matrix, or statspai.spatial.weights.W`	n x n spatial weights, rows in the order of `data`.	required
`data`	`DataFrame`		required
`formula`	`str`	`"y ~ x1 + x2"`; a constant is added, and the spatial lag of each regressor (not of the constant) enters as `W_<name>`.	required
`row_normalize`	`bool`		`True`
`alpha`	`float`		`0.05`

Returns:

Type	Description
`EconometricResults`	`params` `[const, x..., W_x..., lambda]` with OIM standard errors (Stata `spregress, ml errorlag(W) ivarlag(W: x...)`); `diagnostics['sigma2']` and its SE `se_sigma2`.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 40
>>> W = np.zeros((n, n))
>>> for i in range(n):
...     W[i, (i - 1) % n] = W[i, (i + 1) % n] = 1.0
>>> x = rng.normal(size=n)
>>> df = pd.DataFrame({"y": 1 + x + rng.normal(size=n), "x": x})
>>> list(sp.sdem(W, df, "y ~ x").params.index)
['const', 'x', 'W_x', 'lambda']

sar_gmm ¶

sar_gmm(W: Any, data: DataFrame, formula: str, row_normalize: bool = True, robust: Optional[str] = None, w_lags: int = 1) -> EconometricResults

Kelejian-Prucha (1998) 2SLS for SAR with spatial-lag instruments.

Instruments: [X, W X, …, W^w_lags X] (dropping constant duplicates). w_lags=1 (default) matches spreg.GM_Lag default. Endogenous regressor: W Y.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 40
>>> W = np.zeros((n, n))  # ring contiguity: each unit linked to neighbours
>>> for i in range(n):
...     W[i, (i - 1) % n] = 1.0
...     W[i, (i + 1) % n] = 1.0
>>> x = rng.normal(size=n)
>>> df = pd.DataFrame({"y": 1.0 + 0.8 * x + rng.normal(scale=0.5, size=n),
...                    "x": x})
>>> res = sp.sar_gmm(W, data=df, formula="y ~ x")
>>> list(res.params.index)
['const', 'x', 'rho']

References

kelejian1998generalized

sarar_gmm ¶

sarar_gmm(W: Any, data: DataFrame, formula: str, row_normalize: bool = True, robust: Optional[str] = None, sig2n_k: bool = False, w_lags: int = 2) -> EconometricResults

SARAR (spatial lag + spatial error) by generalized spatial 2SLS.

The Kelejian-Prucha (1998) GS2SLS estimator for y = rho W y + X beta + u, u = lambda W u + e, in the three steps of Kelejian and Prucha (1998):

2SLS of y on [W y, X] with instruments H = [X, W X, W^2 X] (lags of the non-constant columns).
lambda from the Kelejian-Prucha (1999) generalized moments on the 2SLS residuals (unweighted nonlinear least squares).
Spatial Cochrane-Orcutt: 2SLS of y - lambda W y on [(I - lambda W) W y, (I - lambda W) X], instruments [(I - lambda W) X, W X, W^2 X]. With a row-standardised W these span the same space as H, the instruments Kelejian and Prucha (1998) use in the third step; only the regressors are filtered.

beta, rho and their standard errors come from step 3. The GM step does not deliver a standard error for lambda; it is reported as NaN (R spatialreg::gstsls likewise reports no standard error for lambda).

Parameters:

Name	Type	Description	Default
`W`	array-like, sparse matrix or :class:`W`	Spatial weights.	required
`data`	`DataFrame`		required
`formula`	`str`	`"y ~ x1 + x2"`; a constant is always added.	required
`row_normalize`	`bool`	Row-standardise `W` (`spdep` style `"W"`).	`True`
`robust`	`(None, 'het')`	`"het"`: HC0 sandwich in both 2SLS steps (`gstsls(robust = TRUE)`). The `lambda` moments are still the homoskedastic Kelejian-Prucha (1999) ones; this is not the Arraiz et al. (2010) heteroskedastic GM estimator of `sphet::spreg(het = TRUE)`.	`None`
`sig2n_k`	`bool`	Divide the residual sum of squares by `n - k` instead of `n` in the homoskedastic variance (`gstsls(sig2n_k = TRUE)`).	`False`
`w_lags`	`int`	Spatial lags of `X` used as instruments: `[W X, ..., W^w_lags X]`. `2` is Kelejian-Prucha's `H = [X, WX, W^2 X]` and what R `spatialreg::gstsls` uses (run as a black box on the columbus data, its estimates and standard errors are identical with `W2X = TRUE` and `W2X = FALSE`); `1` is PySAL `spreg.GM_Combo`'s default.	`2`

Returns:

Type	Description
`EconometricResults`	`params` `[const, x..., rho, lambda]`.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 40
>>> W = np.zeros((n, n))  # ring contiguity: each unit linked to neighbours
>>> for i in range(n):
...     W[i, (i - 1) % n] = 1.0
...     W[i, (i + 1) % n] = 1.0
>>> x = rng.normal(size=n)
>>> df = pd.DataFrame({"y": 1.0 + 0.8 * x + rng.normal(scale=0.5, size=n),
...                    "x": x})
>>> res = sp.sarar_gmm(W, data=df, formula="y ~ x")
>>> list(res.params.index)
['const', 'x', 'rho', 'lambda']

References

kelejian1998generalized

sem_gmm ¶

sem_gmm(W: Any, data: DataFrame, formula: str, row_normalize: bool = True, robust: Optional[str] = None) -> EconometricResults

Kelejian-Prucha (1999) GMM for the spatial-error parameter λ.

Stage 1 — OLS on y = Xβ + u ⇒ residuals u. Stage 2 — minimise sum of squared KP moment residuals over (λ, σ²). Stage 3 — GLS with (I - λ̂ W) → final β̂.

Parameters:

Name	Type	Description	Default
`robust`	`(None, 'het')`	When `"het"`, the final β covariance uses the heteroscedasticity- robust sandwich (Arraiz et al. 2010 / `spreg.GM_Error_Het`).	`None`

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 40
>>> W = np.zeros((n, n))  # ring contiguity: each unit linked to neighbours
>>> for i in range(n):
...     W[i, (i - 1) % n] = 1.0
...     W[i, (i + 1) % n] = 1.0
>>> x = rng.normal(size=n)
>>> df = pd.DataFrame({"y": 1.0 + 0.8 * x + rng.normal(scale=0.5, size=n),
...                    "x": x})
>>> res = sp.sem_gmm(W, data=df, formula="y ~ x")
>>> list(res.params.index)
['const', 'x', 'lambda']

impacts ¶

impacts(result: Any, n_sim: int = 1000, seed: Optional[int] = None) -> DataFrame

Direct, indirect and total impacts of a spatial regression, with SEs.

Parameters:

Name	Type	Description	Default
`result`	`EconometricResults`	Output of `sp.sar`, `sp.sdm`, `sp.sac`, `sp.gns` (simulated SEs) or `sp.slx`, `sp.sdem` (exact SEs).	required
`n_sim`	`int`	Number of Monte-Carlo draws for simulated SEs. Not used for SLX / SDEM.	`1000`
`seed`	`int`	RNG seed. Not used for SLX / SDEM.	`None`

Returns:

Type	Description
`DataFrame with one row per non-constant regressor and columns`
``Direct``, ``SE_Direct``, ``Indirect``, ``SE_Indirect``, ``Total``,
``SE_Total``. ``attrs['se_method']`` is ``'simulation'`` or
``'analytic'``; ``attrs['n_draws_discarded']`` counts the draws whose
`spatial parameter fell outside the admissible interval`
``(1 / omega_min, 1 / omega_max)`` of the eigenvalues of ``W`` (real
parts), which are dropped as ``spatialreg`` drops them.

Notes

The total impact is (1/n) 1' S_k 1 (LeSage and Pace 2009). For a row-standardised W it equals beta_k / (1 - rho) in a SAR model; for any other W it does not. spatialreg::impacts on a lagsarlm fit with weights that are not row-standardised reports beta_k / (1 - rho) all the same (spatialreg 1.4.3, checked on a binary ring), so its indirect and total impacts differ from these there; its direct impacts agree.

Through 1.40.2 the admissible interval was fixed at (-0.999, 0.999), right only for a row-standardised W with smallest eigenvalue -1: draws past 1 / omega_max of an unstandardised W were kept and inflated the SEs, and valid draws below -0.999 were dropped.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(7)
>>> n = 120
>>> coords = rng.uniform(size=(n, 2))
>>> w = sp.knn_weights(coords, k=6)
>>> w.transform = "R"
>>> W = w.sparse.toarray()
>>> x = rng.standard_normal(n)
>>> y = np.linalg.solve(np.eye(n) - 0.5 * W,
...                     1 + 2 * x + rng.standard_normal(n))
>>> res = sp.sar(w, pd.DataFrame({"y": y, "x": x}), "y ~ x")
>>> imp = sp.impacts(res, n_sim=200, seed=0)
>>> list(imp.columns)
['Direct', 'SE_Direct', 'Indirect', 'SE_Indirect', 'Total', 'SE_Total']
>>> list(imp.index)
['x']

sac ¶

sac(W: ArrayOrW, data: DataFrame, formula: str, row_normalize: bool = True, alpha: float = 0.05) -> EconometricResults

SAC / SARAR: combined spatial lag + spatial error model.

y = rho W y + X b + u, u = lambda W u + e, e ~ N(0, sigma^2 I). (b, sigma^2) are profiled out and (rho, lambda) maximise the concentrated log-likelihood. Standard errors invert the observed information of the full log-likelihood in (b, rho, lambda, sigma^2) (Stata spregress, ml dvarlag(W) errorlag(W), default vce(oim)), so they account for the covariance between the regression and the spatial parameters. Through 1.39 the regression SEs treated rho and lambda as known and the spatial SEs read the diagonal of a numerical Hessian of the concentrated likelihood; both were too small.

Parameters:

Name	Type	Description	Default
`W`	`ndarray, scipy.sparse matrix, or statspai.spatial.weights.W`	n x n spatial weights, rows in the order of `data`.	required
`data`	`DataFrame`		required
`formula`	`str`	`"y ~ x1 + x2"` with plain column names; a constant is added.	required
`row_normalize`	`bool`		`True`
`alpha`	`float`		`0.05`

Returns:

Type	Description
`EconometricResults`	`params` `[const, x..., rho, lambda]` with OIM standard errors; `model_info['vcov']` (their covariance), `model_info['sigma2']` and `model_info['se_sigma2']`, `model_info['loglik']`.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 30
>>> W = np.zeros((n, n))
>>> for i in range(n):
...     W[i, (i - 1) % n] = 1.0
...     W[i, (i + 1) % n] = 1.0
>>> x = rng.normal(size=n)
>>> df = pd.DataFrame({"y": 1.0 + 0.8 * x + rng.normal(scale=0.5, size=n),
...                    "x": x})
>>> res = sp.sac(W, data=df, formula="y ~ x")
>>> list(res.params.index)  # both spatial coefficients estimated jointly
['const', 'x', 'rho', 'lambda']

References

anselin1988spatial

slx ¶

slx(W: ArrayOrW, data: DataFrame, formula: str, row_normalize: bool = True, alpha: float = 0.05) -> EconometricResults

Spatial Lag of X (SLX) model: Y = Xβ + WXθ + ε.

No autoregressive term. Estimated by ordinary least squares on the design matrix augmented with spatially-lagged covariates (skip the constant).

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 30
>>> W = np.zeros((n, n))
>>> for i in range(n):
...     W[i, (i - 1) % n] = 1.0
...     W[i, (i + 1) % n] = 1.0
>>> x = rng.normal(size=n)
>>> df = pd.DataFrame({"y": 1.0 + 0.8 * x + rng.normal(scale=0.5, size=n),
...                    "x": x})
>>> res = sp.slx(W, data=df, formula="y ~ x")
>>> list(res.params.index)
['const', 'x', 'W_x']

References

anselin1988spatial

spatial_panel ¶

spatial_panel(data: DataFrame, formula: str, entity: str, time: str, W: Any, model: ModelKind = 'sar', effects: EffectKind = 'fe', row_normalize: bool = True, vce: str = 'information', twoways_lag: str = 'splm') -> SpatialPanelResult

Fit a spatial panel model by concentrated ML (Elhorst 2014).

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`	Long-format panel; must be balanced.	required
`formula`	`str`	`"y ~ x1 + x2"`. Constant is dropped (absorbed by the entity FE).	required
`entity`	`str`	Column names identifying entity and time dimensions.	required
`time`	`str`	Column names identifying entity and time dimensions.	required
`W`	sparse matrix, ndarray, or :class:`statspai.spatial.weights.W`	N × N spatial weights matrix aligned with `sorted(data[entity].unique())`.	required
`model`	`('sar', 'sem', 'sdm', 'sac')`	`"sac"` (spatial lag and spatial error) needs `effects="re"`.	`"sar"`
`effects`	`('fe', 'twoways', 're')`	`"fe"` = entity FE only. `"twoways"` = entity + time demeaning. `"re"` = random unit effects `c_i ~ N(0, sigma_u^2)` by full maximum likelihood with a constant, as Stata `spxtregress, re` (`dvarlag` = SAR, `errorlag` = SEM, both = SAC, `dvarlag` + `ivarlag` = SDM): estimates and OIM standard errors equal Stata 18 to 1e-6. `vce` and `twoways_lag` do not apply; the standard deviations `sigma_u` / `sigma_e` and their SEs are in `result.variance_components`. Stata's `dvarlag` coefficient is `rho` here and its `errorlag` coefficient `lambda`.	`"fe"`
`row_normalize`	`bool`	Row-standardise `W`.	`True`
`vce`	`('information', 'oim')`	`"information"`: inverse of the analytic information matrix of `(sigma^2, rho, beta)` evaluated at the estimates, as `splm::spml(model = "within")` reports. `"oim"`: inverse of the observed (analytic) Hessian of the full log-likelihood -- the `vce(oim)` default of Stata's `xsmle`. The two differ in the `rho`-`rho` term (expected `tr(WA'WA)` + `beta'X'(WA)'WAX beta` against observed `(Wy)'(Wy)`) and so in every SE.	`"information"`
`twoways_lag`	`('splm', 'within')`	Spatial lag regressor of SAR / SDM under `effects="twoways"` (ignored otherwise: under entity effects the two coincide). `"splm"`: `W (Q y)`, the lag of the two-way demeaned `y`, as `splm::spml` computes it (and as `xsmle`'s log-likelihood evaluator does). `"within"`: `Q (W y)`, the two-way within transform of the lag, which is what the concentrated likelihood of the model with unit and time dummies (Lee & Yu 2010, eq. 21, the "direct approach") contains. For a row-standardised `W` the two differ by a per-period constant `(1/N) 1' W Q y_t` unless `W` is also column-stochastic, so `"splm"`'s SSE carries an extra `rho^2 N sum_t m_t^2` term. Both are direct-approach estimators with an O(1/N) bias in `rho` that does not vanish as `T` grows (see `docs/dev/campaign_phase3`). `"within"` supports `vce="oim"` only.	`"splm"`

Returns:

Type	Description
`SpatialPanelResult`	`params` / `std_errors` indexed `[x..., rho]` (`[x..., W_x..., rho]` for SDM, `[x..., lambda]` for SEM), `sigma2 = e'e / NT` and the full log-likelihood.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> N, T = 8, 6
>>> W = np.zeros((N, N))  # ring contiguity weights
>>> for i in range(N):
...     W[i, (i - 1) % N] = 1.0
...     W[i, (i + 1) % N] = 1.0
>>> rows = []
>>> for i in range(N):
...     fe = rng.normal()
...     for t in range(T):
...         x = rng.normal()
...         y = fe + 0.5 * x + rng.normal(0, 0.3)
...         rows.append({"id": i, "year": t, "y": y, "x": x})
>>> df = pd.DataFrame(rows)
>>> res = sp.spatial_panel(
...     df, "y ~ x", entity="id", time="year", W=W, model="sar",
... )
>>> res.spatial_param
'rho'
>>> (res.N, res.T)
(8, 6)
>>> bool(np.isfinite(res.params["x"]))
True

distance_to_feature ¶

distance_to_feature(points_gdf: Any, target_gdf: Any, unit: str = 'km', allow_geographic: bool = False) -> Series

Distance from each point to the nearest target feature.

The canonical use is "distance from each county seat to the canal", the continuous treatment-intensity variable in a historical-GIS design.

Parameters:

Name	Type	Description	Default
`points_gdf`	`GeoDataFrame`	Point features.	required
`target_gdf`	`GeoDataFrame`	Target layer of any geometry type; the nearest distance to the union of its features is returned. Points inside a target polygon get `0`.	required
`unit`	`('km', 'm')`	Output unit.	`'km'`
`allow_geographic`	`bool`	Opt out of the projected-CRS requirement.	`False`

Returns:

Type	Description
`Series`	Named `distance_<unit>`, indexed like `points_gdf`.

Raises:

Type	Description
`ValueError`	On any of the CRS / geometry problems listed in :func:`line_length_in_polygon`.

Examples:

>>> import statspai as sp
>>> sp.distance_to_feature(county_seats, canal, unit="km")

line_length_in_polygon ¶

line_length_in_polygon(lines_gdf: Any, polygons_gdf: Any, polygon_id: Optional[str] = None, unit: str = 'km', allow_geographic: bool = False) -> DataFrame

Length of linear features falling inside each polygon.

The canonical use is "kilometres of canal (or railway, or road) per county", the exposure variable in a historical-GIS panel.

Parameters:

Name	Type	Description	Default
`lines_gdf`	`GeoDataFrame`	Linear features. Must be `LineString` / `MultiLineString`.	required
`polygons_gdf`	`GeoDataFrame`	Areal units. Must be `Polygon` / `MultiPolygon`.	required
`polygon_id`	`str`	Column in `polygons_gdf` identifying each unit. Defaults to the index of `polygons_gdf`.	`None`
`unit`	`('km', 'm')`	Output unit.	`'km'`
`allow_geographic`	`bool`	Opt out of the projected-CRS requirement. Only set this if you have verified the CRS units are metres.	`False`

Returns:

Type	Description
`DataFrame`	One row per polygon with the id column and `line_length_<unit>`. Polygons with no intersecting line get `0.0`, not `NaN`.

Raises:

Type	Description
`ValueError`	If either frame has no CRS, a geographic CRS (and `allow_geographic` is False), a non-metre projected CRS, mismatched CRS, or missing / empty / invalid geometries.

Examples:

>>> import statspai as sp
>>> canal = canal.to_crs("EPSG:4479")
>>> counties = counties.to_crs("EPSG:4479")
>>> sp.line_length_in_polygon(canal, counties, polygon_id="county_id")
...

share_within_buffer ¶

share_within_buffer(points_gdf: Any, lines_gdf: Any, buffer_km: float = 10.0, group_col: Optional[str] = None, allow_geographic: bool = False) -> Any

Share of point features lying within a buffer of the linear features.

The canonical use is "share of market towns within 10 km of the canal", either overall or per county.

Parameters:

Name	Type	Description	Default
`points_gdf`	`GeoDataFrame`	Point features (e.g. market towns).	required
`lines_gdf`	`GeoDataFrame`	Linear features to buffer (e.g. the canal).	required
`buffer_km`	`float`	Buffer radius in kilometres. Must be strictly positive.	`10.0`
`group_col`	`str`	Column in `points_gdf` to group by (e.g. `'county_id'`). If omitted the share is computed over all points.	`None`
`allow_geographic`	`bool`	Opt out of the projected-CRS requirement.	`False`

Returns:

Type	Description
`float or DataFrame`	Without `group_col`, the overall share in `[0, 1]`. With it, a DataFrame with `group_col`, `n_points`, `n_within`, `share`.

Raises:

Type	Description
`ValueError`	On non-positive `buffer_km`, or any of the CRS / geometry problems listed in :func:`line_length_in_polygon`.

Examples:

>>> import statspai as sp
>>> sp.share_within_buffer(towns, canal, buffer_km=10, group_col="county_id")
...

block_weights ¶

block_weights(regimes: Sequence[Any]) -> W

Block (regime) spatial weights — full connectivity within each group.

Builds a weights object in which every unit is a neighbour of every other unit sharing the same regime label, and of no unit outside it. Useful for encoding membership in a region, cluster, or treatment block as a spatial structure.

Parameters:

Name	Type	Description	Default
`regimes`	`array_like`	Length-`n` array of group labels (any hashable dtype). Units with equal labels become mutual neighbours.	required

Returns:

Type	Description
`W`	Spatial weights object with block-diagonal connectivity.

Examples:

>>> import statspai as sp
>>> w = sp.block_weights([0, 0, 1, 1, 1])
>>> sorted(w.neighbors[2])
[3, 4]

distance_band ¶

distance_band(coords: ndarray, threshold: float, binary: bool = True) -> W

Distance-band spatial weights — connect units within a fixed radius.

Every unit is connected to all other units within Euclidean distance threshold. Weights are either binary (1 for every neighbour) or inverse-distance (1/d).

Parameters:

Name	Type	Description	Default
`coords`	`ndarray`	`(n, d)` array of coordinates.	required
`threshold`	`float`	Distance band radius. Units within this distance become neighbours. Choose at least the maximum nearest-neighbour distance to avoid islands (units with no neighbours).	required
`binary`	`bool`	If True all neighbour weights are `1.0`; if False they are inverse Euclidean distance `1/d` (coincident points are dropped via an infinite distance).	`True`

Returns:

Type	Description
`W`	Spatial weights object; neighbour counts vary by local density.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> coords = np.random.default_rng(0).random((30, 2))
>>> w = sp.distance_band(coords, threshold=0.3, binary=False)
>>> w.n
30

See Also

knn_weights : Fixed neighbour count instead of fixed radius.

kernel_weights ¶

kernel_weights(coords: ndarray, bandwidth: float, kernel: Literal['gaussian', 'bisquare', 'triangular'] = 'gaussian', fixed: bool = True) -> W

Kernel (distance-decay) spatial weights from point coordinates.

Assigns each neighbour a continuous weight that decays with distance according to a kernel function. Supports a fixed bandwidth (same radius everywhere) or an adaptive bandwidth (the bandwidth-th nearest neighbour distance, varying with local density).

Parameters:

Name	Type	Description	Default
`coords`	`ndarray`	`(n, d)` array of coordinates.	required
`bandwidth`	`float`	With `fixed=True`, the bandwidth distance `h` (for the bisquare and triangular kernels, neighbours beyond `h` get zero weight; the gaussian kernel has infinite support so all units are included). With `fixed=False`, the integer number `k` of nearest neighbours: each unit is connected to its `k` nearest other units and `u = d / d_k` with `d_k` the distance to the `k`-th of them (so the bisquare weight of the `k`-th neighbour is 0; this is `GWmodel::gw.weight(adaptive = TRUE)` with `bw = k + 1`, which counts the unit itself). For every kernel, including the Gaussian, only those `k` units are connected -- unlike the adaptive Gaussian of :func:`statspai.gwr`, which weights all points.	required
`kernel`	`('gaussian', 'bisquare', 'triangular')`	Kernel shape. `gaussian` `exp(-u^2/2)`; `bisquare` `(1-u^2)^2` for `u<1`; `triangular` `1-u` for `u<1`, where `u = d / bandwidth`.	`"gaussian"`
`fixed`	`bool`	Fixed (True) versus adaptive (False) bandwidth.	`True`

Returns:

Type	Description
`W`	Spatial weights object with continuous kernel weights. The fixed bisquare is spdep's `nb2listwdist(type = "dpd", alpha = 2)` on `dnearneigh(0, bandwidth)`; `w.transform = "R"` row-standardises the kernel values themselves.

Raises:

Type	Description
`ValueError`	If `kernel` is not one of the supported names.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> coords = np.random.default_rng(0).random((30, 2))
>>> w = sp.kernel_weights(coords, bandwidth=0.4, kernel="bisquare")
>>> w.n
30

See Also

distance_band : Binary / inverse-distance fixed-radius weights. knn_weights : Uniform k-nearest-neighbour weights.

knn_weights ¶

knn_weights(coords: ndarray, k: int = 5) -> W

k-nearest-neighbour spatial weights from point coordinates.

Each unit is connected to its k closest units by Euclidean distance (the unit itself is excluded). The resulting graph is generally asymmetric — being someone's nearest neighbour is not reciprocal.

Parameters:

Name	Type	Description	Default
`coords`	`ndarray`	`(n, d)` array of coordinates (typically `d=2`).	required
`k`	`int`	Number of nearest neighbours per unit. Must satisfy `0 < k < n`.	`5`

Returns:

Type	Description
`W`	Spatial weights object with exactly `k` neighbours per unit.

Raises:

Type	Description
`ValueError`	If `k` is not in the open interval `(0, n)`.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> coords = np.random.default_rng(0).random((30, 2))
>>> w = sp.knn_weights(coords, k=4)
>>> len(w.neighbors[0])
4

See Also

distance_band : Connect all units within a fixed radius. kernel_weights : Continuous distance-decay weights.

queen_weights ¶

queen_weights(gdf: Any) -> W

Queen-contiguity spatial weights from polygon geometries.

Two areal units are queen-contiguous if their geometries share at least one boundary point (an edge or a single vertex). Requires geopandas.

Parameters:

Name	Type	Description	Default
`gdf`	`GeoDataFrame`	Layer of polygon geometries (accessed via `gdf.geometry` and the spatial index `gdf.sindex`).	required

Returns:

Type	Description
`W`	Spatial weights object whose neighbour lists encode queen contiguity.

Raises:

Type	Description
`ImportError`	If `geopandas` is not installed.

Examples:

>>> import statspai as sp
>>> import geopandas as gpd
>>> gdf = gpd.read_file("counties.shp")
>>> w = sp.queen_weights(gdf)

See Also

rook_weights : Edge-only contiguity (excludes shared vertices).

rook_weights ¶

rook_weights(gdf: Any) -> W

Rook-contiguity spatial weights from polygon geometries.

Two areal units are rook-contiguous if their geometries share a boundary segment of positive length (a shared edge); units touching only at a single vertex are not neighbours. Requires geopandas.

Parameters:

Name	Type	Description	Default
`gdf`	`GeoDataFrame`	Layer of polygon geometries.	required

Returns:

Type	Description
`W`	Spatial weights object whose neighbour lists encode rook contiguity.

Raises:

Type	Description
`ImportError`	If `geopandas` is not installed.

Examples:

>>> import statspai as sp
>>> import geopandas as gpd
>>> w = sp.rook_weights(gpd.read_file("counties.shp"))

See Also

queen_weights : Contiguity that also counts shared vertices.

statspai.spatial¶

spatial ¶

SpatialDiDResult dataclass ¶

plot ¶

SpatialIVResult dataclass ¶

SpatialModel ¶

fit ¶

W ¶

spatial_did ¶

getis_ord_g ¶

getis_ord_local ¶

moran_local ¶

lisa_cluster_map ¶

moran_plot ¶

gwr_bandwidth ¶

spatial_iv ¶

sar ¶

sdm ¶

sem ¶

lm_tests ¶

moran_residuals ¶

gns ¶

sdem ¶

sar_gmm ¶

sarar_gmm ¶

sem_gmm ¶

impacts ¶

sac ¶

slx ¶

spatial_panel ¶

distance_to_feature ¶

line_length_in_polygon ¶

share_within_buffer ¶

block_weights ¶

distance_band ¶

kernel_weights ¶

knn_weights ¶

queen_weights ¶

rook_weights ¶

`statspai.spatial`¶

SpatialDiDResult `dataclass` ¶

SpatialIVResult `dataclass` ¶