Spatial dependence and LISA

Who counts as a neighbor, whether the map clusters, and where

Every spatial statistic answers a question relative to a definition of neighborhood. This article walks the exploratory spatial data analysis (ESDA) arc on the bundled India case study — 520 districts observed by satellite nighttime lights, 1996–2010: encode the neighborhood as a weights matrix, audit it visually, test whether luminosity clusters globally (Moran’s I), locate the clusters locally (LISA), and track the clustering over time.

import warnings

warnings.filterwarnings("ignore")

import geometrics as gm

gdf, df, df_dict = gm.data.load_india()
df = gm.set_labels(df, df_dict, set_panel=True)

# The paper's weights: 6 nearest neighbors on plain lon/lat centroids (crs=None)
w = gm.make_weights(gdf, method="knn", k=6, crs=None)

What a weights matrix is

The spatial weights matrix W is the formal answer to “who is whose neighbor”: cell (i, j) is nonzero when district j counts as a neighbor of district i. It is not a nuisance parameter — Moran’s I, LISA, and every spatial regression are defined conditional on the chosen W. The built-in explainer says it best:

print(gm.explain("spatial_weights").to_markdown()[:800], "...")

### Spatial weights (W)

**What it is.** A spatial weights matrix **W** encodes who counts as whose *neighbor*: cell (i, j) is nonzero when unit j is a neighbor of unit i. Contiguity criteria (**queen**: shared border or corner; **rook**: shared border only) suit administrative polygons; **k-nearest-neighbors** guarantees every unit exactly k neighbors; **distance bands** and **inverse-distance** weights encode decay with separation. Almost every spatial statistic — Moran's I, LISA, spatial regression, the spatial Markov chain — is defined *relative to a chosen W*: the weights are part of the hypothesis, not a nuisance detail.

**When to use it.** Build one W per analysis with `make_weights` and reuse it everywhere so results are comparable. Queen contiguity is the parameter-free default f ...

Building weights: make_weights

make_weights covers the standard families with the library’s conventions baked in: entity ids as the weight ids, row standardization, and a human-readable one-liner recorded on w.geometrics_meta["spec"] (every spatial result carries it forward as w_spec, so figures and tables always document their W). Three variants on the India map:

variants = {
    "queen": gm.make_weights(gdf, method="queen"),
    "knn6": w,
    "inverse distance": gm.make_weights(gdf, method="inverse_distance", power=2),
}
for name, wx in variants.items():
    print(f"{name:>18}: {wx.geometrics_meta['spec']}")

             queen: queen contiguity (1 island(s) attached to nearest neighbor), row-standardized, n=520
              knn6: 6-nearest-neighbor (geographic centroids), row-standardized, n=520
  inverse distance: inverse distance (power 2, band threshold 163388), row-standardized, n=520

Each spec string tells a small story:

• Queen contiguity connects districts that share a border or corner. One district — Mumbai, an island city that shares no boundary in this topology — has no queen neighbor at all, and make_weights attaches it to its nearest neighbor automatically (with a GeometricsWarning naming it), because a zero-neighbor unit has no spatial lag and breaks estimators downstream. The fix is recorded in the spec.
• k-nearest neighbors guarantees every district exactly k neighbors — useful when polygons vary wildly in size. crs=None reproduces the source paper’s lon/lat-centroid construction; the geometrics default (crs="auto") would project to a metric CRS first.
• Inverse distance encodes smooth decay, \(1/d^p\) within a radius. With threshold=None the band is the smallest radius that leaves no district isolated — about 163 km here.

Audit the graph before trusting it

explore_connectivity_map draws the neighbor graph over the map and summarizes its health — the standard visual check of a W before any statistic is computed on it:

conn = gm.explore_connectivity_map(gdf, w=w)
print(f"{conn.n_units} units, {conn.mean_neighbors:.0f} neighbors each "
      f"({conn.min_neighbors}-{conn.max_neighbors}); "
      f"{conn.pct_nonzero:.2f}% of pairs connected; "
      f"{conn.n_components} connected component(s); islands: {list(conn.islands)}")

520 units, 6 neighbors each (6-6); 1.15% of pairs connected; 1 connected component(s); islands: []

conn.fig

The scalars are the audit: a single connected component (global statistics implicitly assume every unit can reach every other through a chain of neighbors), no islands, and a sparse graph — each spatial lag averages over a small, local set of six districts. A companion histogram of neighbor counts is available as conn.fig_hist (trivially a spike at 6 for k-NN weights, more informative for contiguity).

print(conn.interpret())

The spatial weights graph (**6-nearest-neighbor (geographic centroids), row-standardized, n=520**) links 520 units with 6 neighbors each on average (ranging from 6 to 6); 1.15% of all possible unit pairs are connected, so each spatial lag averages over a small, local set of neighbors.

The graph forms a **single connected component** — every unit can reach every other through a chain of neighbors, which is what global spatial statistics (Moran's I, spatial models) implicitly assume.

No unit is an island: every unit has at least one neighbor, so spatial lags are defined everywhere.

_These are associations, not causal effects. A causal reading needs a research design — see `explain('correlation_vs_causation')`._

Is initial luminosity clustered? The Moran scatterplot

explore_moran_plot z-standardizes the variable, plots each district’s value against the average of its neighbors’ values (the spatial lag), and colors the four quadrants: High-High and Low-Low (clustering of similars) on the diagonal, Low-High and High-Low (spatial outliers) off it. Under row-standardized weights the slope of the fitted line equals global Moran’s I, and a stat box in the corner reports I, its expectation under spatial randomness, and the permutation pseudo p-value:

moran = gm.explore_moran_plot(df, "log_ntl_pc_1996", gdf=gdf, w=w, period=1996)
moran.fig

moran.glance().round(3)

	var	moran_i	expected_i	z_sim	p_sim	permutations	n_obs
0	log_ntl_pc_1996	0.73	-0.002	30.323	0.001	999	520

Initial luminosity per capita is strongly and positively autocorrelated — Moran’s I ≈ 0.73 against an expectation near zero, the value the source paper reports. Most districts sit in the two clustering quadrants:

print(moran.df["quadrant"].value_counts().to_string())
print()
print(moran.interpret())

quadrant
HH    252
LL    195
HL     42
LH     31

Global Moran's I for **log_ntl_pc_1996** in 1996 is 0.73, against an expectation of -0.00193 under spatial randomness — statistically significant at the 1% level (pseudo p = 0.001 from 999 permutations, under 6-nearest-neighbor (geographic centroids), row-standardized, n=520).

The dependence is **positive**: similar values cluster in space — high values sit next to high values and low next to low, so the map shows contiguous patches rather than a random scatter.

447 of 520 regions (86%) fall in the clustering quadrants of the scatter (High-High or Low-Low); the rest are surrounded by neighbors unlike themselves. The slope of the fitted line equals Moran's I under row-standardized weights, so a steeper line means stronger clustering.

_These are associations, not causal effects. A causal reading needs a research design — see `explain('correlation_vs_causation')`._

Where does it cluster? LISA

Global Moran’s I says that the map clusters; local Moran statistics (LISA) say where. Each district gets a local statistic and its own permutation pseudo p-value; districts significant at alpha receive their quadrant’s cluster label, everything else is “Not significant”:

lisa = gm.explore_lisa_cluster_map(df, "log_ntl_pc_1996", gdf=gdf, w=w, period=1996)
lisa.fig

print(lisa.interpret())

Local Moran statistics (LISA) locate *where* **log_ntl_pc_1996** in 1996 clusters or stands out, under 6-nearest-neighbor (geographic centroids), row-standardized, n=520. The accompanying global Moran's I is 0.73 (pseudo p = 0.001), consistent with overall clustering of similar values.

At the 0.05 significance level, 281 of 520 regions show significant local association: **169 High-High** hot spots (high values surrounded by high neighbors) and **101 Low-Low** cold spots (low surrounded by low) mark clustering, while **9 High-Low** and **2 Low-High** regions are spatial outliers that break with their surroundings. The remaining 239 regions are not significant — their local pattern is compatible with randomness.

LISA pseudo p-values are computed region by region without a multiple-testing adjustment, so treat borderline clusters cautiously and read the map as descriptive of where dependence concentrates.

_These are associations, not causal effects. A causal reading needs a research design — see `explain('correlation_vs_causation')`._

The picture matches the literature: a bright High-High belt and a dim Low-Low block, with only a handful of spatial outliers.

How sensitive are the clusters to alpha?

The cluster labels are a significance mask, so they move with the threshold. Re-running at alpha=0.01 keeps the geography but shrinks the significant set:

import pandas as pd

lisa01 = gm.explore_lisa_cluster_map(
    df, "log_ntl_pc_1996", gdf=gdf, w=w, period=1996, alpha=0.01
)
pd.DataFrame(
    {
        "alpha = 0.05": [lisa.n_hh, lisa.n_ll, lisa.n_lh, lisa.n_hl, lisa.n_ns],
        "alpha = 0.01": [lisa01.n_hh, lisa01.n_ll, lisa01.n_lh, lisa01.n_hl, lisa01.n_ns],
    },
    index=["High-High", "Low-Low", "Low-High", "High-Low", "Not significant"],
)

	alpha = 0.05	alpha = 0.01
High-High	169	95
Low-Low	101	78
Low-High	2	0
High-Low	9	2
Not significant	239	345

Tightening alpha from 0.05 to 0.01 drops the significant count from 281 to 175 districts. The hot and cold cores survive; the borderline fringe does not — which is exactly how a LISA map should be read.

Is the clustering strengthening? Moran’s I over time

explore_moran_over_time pivots the panel wide, keeps the one entity set with complete data in every period (so the same W applies throughout and the values are comparable), and re-tests each of the six satellite years:

mot = gm.explore_moran_over_time(df, "ntl_total", gdf=gdf, w=w)
mot.df.round(3)

	period	moran_i	z_sim	p_sim	n_obs
0	1996	0.429	18.325	0.001	520
1	1999	0.405	17.567	0.001	520
2	2000	0.450	18.946	0.001	520
3	2004	0.388	16.725	0.001	520
4	2005	0.397	17.044	0.001	520
5	2010	0.416	17.361	0.001	520

mot.fig

print(mot.interpret())

Global Moran's I for **ntl_total** is tracked over 6 periods (1996 to 2010) on a fixed set of regions, under 6-nearest-neighbor (geographic centroids), row-standardized, n=520: it moves from 0.429 to 0.416.

The series is broadly **stable**: the degree to which similar values cluster in space changes little over the window.

Per-period permutation tests flag 6 of 6 periods as significant at the 5% level (filled markers in the figure); open markers are periods where the pattern is compatible with spatial randomness.

_These are associations, not causal effects. A causal reading needs a research design — see `explain('correlation_vs_causation')`._

The answer is: strong throughout, but not strengthening. Moran’s I for total luminosity fluctuates in a narrow band (roughly 0.39–0.45) and every year rejects spatial randomness at the 1% level — spatial dependence is a persistent feature of India’s economic geography over the window, not a trend.

A caveat, and where next

LISA runs one permutation test per district — 520 tests here, with no multiple-testing adjustment, so at alpha = 0.05 roughly 26 “significant” districts would be expected even under complete spatial randomness. Treat the cluster map as descriptive of where dependence concentrates (the alpha sensitivity table above is the practical check), not as 520 confirmatory hypothesis tests. The same warning is built into lisa.interpret().

Detecting dependence is the exploratory half. Modeling it — letting growth respond to neighbors’ outcomes and characteristics, and decomposing associations into direct and spillover components — is the job of the spreg suite, covered in Spatial spillovers: the spreg suite.