Choosing Tolerance Thresholds for Moran’s I Gates

A Moran’s I fidelity gate that compares synthetic output against the source with a hand-picked absolute tolerance such as abs(I_syn - I_src) < 0.05 will reject statistically faithful output at large sample sizes and wave through structurally broken output at small ones, because the sampling variability of Moran’s I is a function of sample size and the weights matrix, not a constant.

Part of Realism Metrics & Evaluation: that stage scores how closely a synthetic layer reproduces the source’s spatial structure, and this page fixes the single most common miscalibration in that stage — treating the acceptance width of an autocorrelation statistic as a fixed number instead of deriving it from the source distribution itself.

Root Cause: A Fixed Tolerance Ignores the Sampling Distribution of I

Moran’s I measures global spatial autocorrelation — the degree to which nearby features carry similar attribute values — against a spatial weights matrix $\mathbf{W}$:

$$I = \frac{N}{W} \cdot \frac{\sum_{i}\sum_{j} w_{ij}\,(x_i - \bar{x})(x_j - \bar{x})}{\sum_{i}(x_i - \bar{x})^2}, \qquad \mathbb{E}[I] = -\frac{1}{N-1}$$

where $N$ is the number of features, $w_{ij}$ the weight between features $i$ and $j$, and $W = \sum_{i}\sum_{j} w_{ij}$ the sum of all weights. Two properties of this statistic break any fixed tolerance.

First, the expected value under the null hypothesis of no autocorrelation is not zero — it is $-1/(N-1)$. For a small synthetic sample of 200 features that offset is about $-0.005$; for a source of 50,000 features it is $-0.00002$. A tolerance centred on zero already carries a size-dependent bias before any real difference is measured.

Second, and decisively, the variance of $I$ under the null shrinks as $N$ grows and depends on the connectivity structure of $\mathbf{W}$ — the number of neighbours each feature has, and how skewed the attribute distribution is (the kurtosis term). Under randomization the null variance scales roughly as $\mathrm{Var}[I] \sim 1/N$, so the standard deviation falls as $1/\sqrt{N}$. A tolerance of $0.05$ that is comfortably inside one standard deviation at $N = 300$ can be twenty standard deviations wide at $N = 100{,}000$. The same number is simultaneously too loose to catch real distortion in a large layer and too tight to permit honest sampling noise in a small one.

The consequence in a generation pipeline is a gate that fails for the wrong reason. It fires on runs whose only “defect” is that they were sampled at a different cardinality than the reference, and it stays green on runs whose autocorrelation has genuinely collapsed but happens to land within a wide fixed band. The fix is to stop guessing the width and instead bootstrap the null band from the source, so the acceptance region is calibrated to the sample size and weights the synthetic layer will actually be evaluated at. This is the same discipline the sibling metric evaluating spatial realism with Wasserstein distance applies to distributional distance: the threshold is derived, never assumed.

Prerequisite Check: Confirm the Fixed Tolerance Is Miscalibrated

Before rewriting the gate, demonstrate that the sampling spread of I depends on N. Pin the toolchain so the weights construction and permutation streams are identical across machines.

pythonCopy
# requirements.txt — pin majors, let patch releases float
numpy==1.26.*
scipy==1.11.*
geopandas==0.14.*     # source/synthetic GeoDataFrames
libpysal==4.*         # spatial weights (KNN, Queen, distance band)
esda==2.*             # Moran's I estimator

pythonCopy
import numpy as np
from libpysal.weights import KNN
from esda.moran import Moran

def moran_i(gdf, value_col: str, k: int = 8) -> float:
    """Global Moran's I under a fixed KNN weights scheme, row-standardized."""
    w = KNN.from_dataframe(gdf, k=k)
    w.transform = "r"                         # row-standardize: W sums to N
    return Moran(gdf[value_col].to_numpy(), w, permutations=0).I

# Resample the SAME source at two cardinalities and watch the spread of I collapse.
rng = np.random.default_rng(20260610)
for n in (300, 30_000):
    vals = [moran_i(source.sample(n, random_state=int(rng.integers(1e9))), "attr")
            for _ in range(200)]
    print(n, f"mean={np.mean(vals):.4f}  sd={np.std(vals):.4f}")
# 300    mean=0.4123  sd=0.0361     <- a 0.05 tolerance is ~1.4 sd: reasonable here
# 30000  mean=0.4488  sd=0.0037     <- the SAME 0.05 tolerance is ~13 sd: far too loose

The standard deviation falls by roughly $\sqrt{100}=10\times$ when N grows $100\times$, exactly as $1/\sqrt{N}$ predicts. A single fixed number cannot serve both runs.

Fix: Bootstrap a Null Band From the Source, Then Gate Against It

Calibrate the acceptance region in three parts: fix the weights scheme so it is identical for source and synthetic, bootstrap the sampling distribution of I from the source at the synthetic’s cardinality, then accept only if the synthetic I lands inside a percentile band of that distribution. The weights scheme must be frozen in a contract exactly as the scoping rules and data contracts stage freezes the CRS — a KNN gate and a Queen-contiguity gate produce different I values on the same data, so the scheme is part of the metric definition.

Step 1 — Freeze the weights scheme

pythonCopy
from dataclasses import dataclass

@dataclass(frozen=True)
class MoranGateSpec:
    value_col: str
    k: int = 8                 # KNN neighbours; part of the metric definition
    transform: str = "r"       # row-standardization
    lo_pct: float = 2.5        # calibrated band edges
    hi_pct: float = 97.5
    n_boot: int = 500          # bootstrap replicates
    base_seed: int = 20260610  # versioned, reproducible

Step 2 — Bootstrap the null band at the synthetic’s sample size

Draw n_boot subsamples of the source, each the size of the synthetic layer, and compute I on each. This is the distribution of I you would expect from data statistically indistinguishable from the source, observed at the exact cardinality the synthetic will be judged at — so it absorbs both the $-1/(N-1)$ offset and the $1/\sqrt{N}$ spread automatically.

pythonCopy
import numpy as np

def bootstrap_null_band(source, spec: MoranGateSpec, n_synth: int) -> tuple[float, float]:
    """Percentile band of Moran's I from source subsamples sized like the synthetic."""
    rng = np.random.default_rng(spec.base_seed)
    n = min(n_synth, len(source))                      # cannot draw more than we have
    boot = np.empty(spec.n_boot)
    for b in range(spec.n_boot):
        sub = source.sample(n, replace=(n_synth > len(source)),
                            random_state=int(rng.integers(1 << 32)))
        boot[b] = moran_i(sub, spec.value_col, spec.k)
    lo, hi = np.percentile(boot, [spec.lo_pct, spec.hi_pct])
    return float(lo), float(hi)

Step 3 — Gate the synthetic against the calibrated band

pythonCopy
from dataclasses import dataclass

@dataclass(frozen=True)
class GateResult:
    passed: bool
    i_synth: float
    band: tuple[float, float]
    reason: str

def morans_i_gate(source, synth, spec: MoranGateSpec) -> GateResult:
    """Accept only if synthetic I lies inside the source-calibrated null band."""
    lo, hi = bootstrap_null_band(source, spec, n_synth=len(synth))
    i_syn = moran_i(synth, spec.value_col, spec.k)
    ok = lo <= i_syn <= hi
    reason = "within source null band" if ok else (
        f"I={i_syn:.4f} outside [{lo:.4f}, {hi:.4f}] "
        f"({'over-clustered' if i_syn > hi else 'under-clustered'})")
    return GateResult(ok, i_syn, (lo, hi), reason)

A synthetic layer whose autocorrelation is higher than the band (i_syn > hi) is over-clustered — often a sign that a point process simulation model has collapsed onto too few attraction centres. One below the band is over-dispersed, the attribute smeared too uniformly across space. The reason string names which failure occurred so the pipeline can route the run to the right diagnostic.

Verification Step: Assert Calibration and Directionality in CI

Wire the gate as a pytest so a regression cannot promote a miscalibrated threshold. These assertions belong in the automated evaluation the CI/CD integration for spatial data stage runs on every generation commit.

pythonCopy
import numpy as np

def test_morans_i_gate_is_calibrated(source):
    spec = MoranGateSpec(value_col="attr")

    # 1. The band tracks sample size: larger N -> narrower band.
    lo_s, hi_s = bootstrap_null_band(source, spec, n_synth=300)
    lo_l, hi_l = bootstrap_null_band(source, spec, n_synth=30_000)
    assert (hi_l - lo_l) < (hi_s - lo_s), "band must narrow as N grows"

    # 2. A faithful resample of the source PASSES at both sizes.
    for n in (300, 30_000):
        faithful = source.sample(n, replace=True, random_state=7)
        assert morans_i_gate(source, faithful, spec).passed

    # 3. A structurally broken layer is REJECTED with the right direction.
    broken = source.copy()
    broken["attr"] = np.sort(broken["attr"].to_numpy())   # force over-clustering
    res = morans_i_gate(source, broken, spec)
    assert not res.passed and "over-clustered" in res.reason

    # 4. Reproducibility: the band is byte-stable for a fixed seed.
    b1 = bootstrap_null_band(source, spec, n_synth=1000)
    b2 = bootstrap_null_band(source, spec, n_synth=1000)
    assert b1 == b2, "bootstrap band not reproducible under fixed seed"

The size-monotonicity check (assertion 1) is the load-bearing one: it is the property a fixed tolerance can never satisfy, and it is what proves the gate is calibrated rather than guessed.

Edge Cases & Gotchas

Antimeridian-spanning weights. KNN and distance-band weights computed on raw longitudes treat features at $+179^\circ$ and $-179^\circ$ as thousands of kilometres apart, so a layer straddling the $\pm180^\circ$ seam gets a fragmented weights graph and an artificially low I. Build the weights in a projected or equal-area CRS whose central meridian sits inside the study area (or on great-circle distances), and freeze that projection in the same contract as k, so the source band and the synthetic are measured on identical geometry.

Synthetic larger than the source. When the synthetic layer has more features than the reference, you cannot draw a same-size subsample without replacement. Bootstrapping with replacement (as the replace= flag above handles) is correct but slightly inflates the band’s tails; if the size ratio exceeds roughly 3:1, prefer resampling the synthetic down to the source size for the comparison, and note the reduced power in the gate report rather than silently widening the band.

Near-constant attributes and the zero-variance trap. If a generator produces an attribute column with (near-)zero variance, the denominator $\sum_i (x_i-\bar{x})^2$ collapses and I is numerically unstable or undefined — yet a naive band can still return finite percentiles from degenerate replicates. Guard with an explicit variance floor before computing I, and fail the gate with a distinct “degenerate attribute” reason so it is never confused with an autocorrelation mismatch.

Frequently Asked Questions

Related

• Realism Metrics & Evaluation — the parent workflow that scores synthetic layers against the source across multiple spatial statistics.
• Evaluating Spatial Realism with Wasserstein Distance — a sibling metric that likewise derives its threshold from the source rather than assuming one.
• CI/CD Integration for Spatial Data — where the calibrated gate runs automatically on every generation commit.
• Point Process Simulation Models — the generators whose mode collapse or over-dispersion this gate is most likely to catch.
• Scoping Rules & Data Contracts — where the weights scheme and CRS that define the metric are frozen.