Add spatial early warning signals (Moran's I) and Empirical Brown's Method p-value combiner

[C95234]

Summary

A spatial early-warning-signal module, developed while building Hélios, a research/outreach project that replicates early-warning-signal methods on new domains and found itself needing this while doing so:

• ewstools/spatial.py: morans_i(), morans_i_permutation_test(), and a SpatialEWS class following the existing MultiTimeSeries conventions (data -> state -> ews, transition-aware). ewstools covers the temporal branch of the critical-slowing-down literature thoroughly but had no equivalent for the spatial branch (Dakos et al., 2010) — confirmed absent by direct inspection of core.py/helpers.py before writing any code.

(An earlier version of this PR also included ewstools/pvalues.py, a correlated-p-value combiner. Removed per review: out of scope for this package, and it had its own bugs — see the review thread below for the full story.)

Testing

• tests/test_spatial.py: hand-computed Moran's I examples (verified independently against esda.moran.Moran), known limiting cases (checkerboard, smooth gradient, constant field), the classical permutation-test property E[I] -> -1/(N-1), a p-value floor test, a NaN-handling warning test, and a negative control confirming that a strengthening spatial gradient raises Moran's I with no actual critical slowing down (see the SpatialEWS docstring's trend-control recipe).
• Full suite, current state of this branch: 50 passed, 1 skipped — 19 new tests, 0 regressions on the 32 tests that existed before this PR.
• CHANGELOG.md updated to match.

Happy to adjust API shape, naming, or scope based on further maintainer feedback.

🤖 Generated with Claude Code

[energyscholar]

Hi @C95234 — thanks for this, and apologies that it sat quiet for a few days. That's on me.

I've read it properly now. You did the things that make a contribution easy to take seriously: you checked whether the gap was real before writing anything, the tests are substantive rather than decorative, you matched the conventions already in the package, and the design follows Dakos et al. (2010) closely. Moran's coefficient at distance 1, plus Kendall tau against the control parameter, is exactly that protocol.

I had a look at Hélios as well. Pushing early-warning methods out into domains they weren't developed on is a particular interest of mine, too, so the cross-domain framing is one I actively care about.

There are a few things I'd want to work through with you — a couple of them scientific rather than cosmetic — but none of them change my interest in the direction you're taking. I'm going through it carefully and talking it over with Tom, and I'll come back early next week with something specific.

Thanks for taking the trouble to contribute this.

— Bruce

[C95234]

Hi @energyscholar — happy to help however's useful here.

[energyscholar]

Sorry this is later than "early next week" — I wanted to check the implementation properly rather
than take it on trust, and that took longer than I said.

What holds up

I went through this carefully, and the core of it is right.

• morans_i computes Moran's I correctly. Your hand-computed chain example (1–2–3–4 with values
  [1,2,3,4], I = (4/6)(2.5/5.0) = 1/3) is exactly the kind of test I want to see in this
  package, and it checks out. I also ran the function against an established implementation
  (esda.moran.Moran) on an 8×8 rook lattice: identical to the last bit of a float64.
• The checkerboard and smooth-gradient cases bracket the statistic from both sides, the shape and
  type guards are there, and the permutation null's mean lands where the theory says it should.
• The whole suite is green: 52 passed, 1 skipped against main at the time you opened it.
• SpatialEWS follows the data/state/ews conventions of MultiTimeSeries rather than
  inventing its own, which is the right call and is not the obvious one.
• The literature framing is right too: this is the spatial branch of the critical-slowing-down
  work, and Dakos et al. (2010) is the correct anchor for it.

I also ran a mutation check on spatial.py — deliberately break a line, see whether a test
notices. 2 of 3 mutations were caught. The survivor is the permutation p-value line, which is
the first item below.

The one thing I'd want addressed before merge: there is no trend control

This is a library gap rather than a mistake in your code, and it is worth stating carefully because
it is the part that a user cannot fix for themselves.

Moran's I is computed on the raw field at each time point. Moran's I is spatially mean-centred, so a
uniform drift in the mean does nothing to it — but a spatial pattern whose amplitude strengthens
over time drives Kendall τ of Moran's I upward with no critical slowing down at all. Measured on
three fields:

field	τ of Moran's I vs time
uniform drift of the mean, no spatial structure change	+0.05
a static spatial gradient, constant amplitude	−0.15
a spatial gradient that strengthens over time	+0.81

Dakos et al. control this by construction — they know what their model is doing. A library
user has no such guarantee.

The consequence for the PR is that test_spatialews_compute_ktau_detects_increasing_trend builds
exactly that field — a column gradient multiplied by t/T — and asserts ktau > 0.5. It is a
negative control mislabelled as a positive one: the test passes because the confound is present,
not because the indicator detected anything. I'd suggest renaming it to say what it demonstrates,
e.g. test_spatialews_ktau_rises_under_strengthening_gradient_without_trend_control, with a comment
saying it is a negative control. The assertion itself needn't change for the merge.

The good news is that the package already has the operation. MultiTimeSeries.detrend()
residuals piped into SpatialEWS take that τ from about +0.5 to about 0.0 on a 6×6 field, with no
new code — measured over five seeds:

seed	τ raw	τ on residuals
0	+0.546	−0.132
1	+0.480	+0.079
2	+0.506	+0.035
3	+0.540	+0.073
4	+0.605	+0.032

I'll add that as a documented recipe beside your PR, together with a shipped null and positive
control so a user can tell a rising Moran's I from a strengthening pattern.

Two honest qualifications, both of which I'll put in the docstring rather than leave implied:

• It removes slowly varying patterns only. A rise in the amplitude of spatially coherent
  fluctuations relative to unit-level noise also raises Moran's I, and a per-unit temporal smoother
  cannot touch it, because every unit's mean is flat by construction. On a 6×6 field over 600 rows
  with a fixed coherent pattern whose amplitude grows against fixed noise, τ is +0.337 raw and
  the detrend changes the mean τ by about 0.001 and no single seed's τ by more than 0.034. That
  needs a reference period, not another detrend; it's on my list.
• τ over the whole record is the right summary but not the whole story. The smoother reflects at
  the record's ends, so on a nearly noise-free field part of the gradient survives in the first and
  last few rows — the residual series is U-shaped and τ over the record is near zero partly by
  symmetry. The recipe's docstring says so.

And one interpretation point, which is why the docstring I'm proposing reads conditionally rather
than "higher = closer to tipping": Rietkerk et al. (2021, Science 374:eabj0359) argue that spatial
self-organisation can mark evasion of tipping rather than approach. A rising Moran's I is a prompt
to build a model of why the pattern forms, not a result.

A p-value that can be exactly zero

In morans_i_permutation_test, p = mean(null >= observed) can be exactly 0.0. Measured on an
8×8 rook lattice with a column-gradient field and n_permutations=500, seed=0: p_value = 0.0,
reported as a genuine significance.

One line fixes it:

"p_value": float((1 + np.sum(null_values >= observed)) / (1 + n_permutations)),

That floors it at 1/(n_permutations + 1) = 1/501 ≈ 0.002 on the same field — the standard
add-one correction (Davison & Hinkley 1997; North, Curtis & Sham 2002). esda uses the same
convention independently: on the same field with 999 permutations it returns
p_sim = 0.001 = 1/(999 + 1).

One note, not a request

morans_i_permutation_test takes a single seed, so SpatialEWS.compute_moran_significance uses
the same permutation set at every time point. The effect is within noise on everything I tried, so
this is a note rather than something to change — but it is worth a docstring line, because a reader
may reasonably expect the p-values across time to be independent draws.

ewstools.pvalues — I'd like to decline this half, on scope rather than quality

Let me start with the motivation, because your own docstring states it well: "are two or more of
these signals unusual at the same time?" That is a coherent question, and correlated indicators
computed from one system are exactly where Fisher's method misleads. I'm not disputing the
statistics.

The structural problem is composition. Your combiner takes one p-value per indicator with the series
behind it — and ewstools doesn't produce a p-value for any indicator: compute_ktau gives a bare τ
and significance is left to the user (surrogates, resampling, scipy.stats.kendalltau). So
everything the combiner would take in comes from outside the package, including for your own
Moran's I (its per-time-point p-value isn't the record-level input the combiner wants). That makes it
general statistical machinery rather than an indicator, which is the line CONTRIBUTING now draws —
combining significance across tests belongs in SciPy or statsmodels.

Two things I found while reading it, which are worth fixing wherever the code ends up living:

• The method is anti-conservative when the indicators are anti-correlated. The empirical
  covariance correction assumes the dependence goes one way. On two perfectly anti-correlated
  inputs the combined p-value is 37× too small; at ρ = −0.8 it is 11× too small. Since
  early-warning indicators can easily anti-correlate, that is a real exposure and worth a guard or
  at least a documented caveat.
• The polynomial is Kost & McDermott's, not Brown's. Poole et al. (2016) describe "Kost and
  McDermott fit a third-order polynomial … 3.263ρ + 0.710ρ² + 0.027ρ³" as the alternative to
  Empirical Brown's Method; EBM's own contribution is to estimate the covariance of the −2 ln p terms
  empirically from the ECDF-transformed data. So the function as written implements Kost &
  McDermott (2002), and a rename plus a docstring correction is all it needs. Worth checking against
  the paper yourself rather than taking my word for it. Your site's "Brown's empirical approach"
  line needs the same correction — mentioning it as a favour, not a criticism.

And the genuine offer: this deserves to exist somewhere. A small standalone package would be the
obvious home, and SciPy is the other. For what it's worth, scipy.stats.combine_pvalues' own
docstring says (verified at scipy 1.15.3): "Extensions such as Brown's method and Kost's method are
not currently implemented." — so the gap is real and acknowledged there. I can't promise anything
about SciPy's appetite for it, but the door isn't closed.

The checklist, exactly — it should be minutes, not a puzzle

(a) The p-value floor. The one-liner above, plus a docstring sentence: "The p-value has a floor
of 1/(n_permutations + 1) (Davison & Hinkley 1997; North, Curtis & Sham 2002) so it can never be
exactly zero."

(b) The test relabel, plus two docstring blocks. Rename
test_spatialews_compute_ktau_detects_increasing_trend →
test_spatialews_ktau_rises_under_strengthening_gradient_without_trend_control, assertion unchanged,
with the comment:

# NEGATIVE CONTROL: this field has no critical slowing down. Its Moran's I
# rises only because the spatial gradient strengthens over time. See the
# Trend control section of the SpatialEWS docstring.

Then two blocks for the SpatialEWS docstring. Both are suggestions — reword freely, it's your
class as much as mine. In Notes:

Trend control. Moran's I is computed on the raw field at each time point. A spatial pattern
whose amplitude changes slowly over time (for example a gradient that strengthens) produces a
trend in Moran's I with no change in the system's dynamics. Where such patterns are plausible,
compute the indicator on residuals from a per-unit temporal detrend: pass the field through
MultiTimeSeries.detrend and hand the residual columns to SpatialEWS (recipe below). The
temporal detrend removes slowly varying spatial patterns only. It does not remove a rise in the
amplitude of spatially coherent fluctuations relative to unit-level noise, which also raises
Moran's I. Dakos et al. (2010, Fig. 6) show two further ways spatial correlation rises with no
change in proximity to a transition: increased connectivity (6a) and increased environmental
heterogeneity (6b). Compare the indicator against a reference period and against the field's
spatial amplitude before reading a trend as slowing down.

Recipe (MultiTimeSeries appends columns to the frame it is given -- pass a copy)::

mts = MultiTimeSeries(df.copy(), transition=t_trans)
mts.detrend(method="Gaussian", bandwidth=0.2)
resid = mts.state[[f"{c}_residuals" for c in df.columns]].set_axis(
    df.columns, axis=1)
sews = SpatialEWS(resid, weights=W, transition=t_trans)

Rows after transition have NaN residuals and are excluded by SpatialEWS in the same way.
Rows containing NaN produce NaN Moran's I (see compute_moran).

and on the weights parameter — this is the one I'd most want in, because nothing else in the class
says it:

weights is fixed for the whole record. If the real network changes over time — units added or
removed, links strengthened — Moran's I trends for that reason alone, independently of the
system's dynamics (Dakos et al. 2010, Fig. 6a). Recompute over a sub-network that is present and
connected throughout. weights[i, j] refers to data.columns[i] and data.columns[j] by
position: reordering, merging or sorting the columns after building weights changes the value
silently; there is no check.

(c) Two lines in ewstools/__init__.py, so the module is importable as the rest of the package
is:

from . import spatial
from .spatial import SpatialEWS

(d) A block in docs/source/ewstools.rst — an ewstools.spatial submodule section with the
same automodule stanza as the others, plus a sentence of description. Without it the module builds
fine and simply never appears on Read the Docs; nothing warns you. I built the docs locally with
warnings-as-errors to check: with the block, 8 morans_i anchors on the page; without it, zero.

(e) The split. Delete CONTRIBUTION_spatial_significance.md, ewstools/pvalues.py and
tests/test_pvalues.py, and in CHANGELOG.md under Unreleased: drop the sentence "See
CONTRIBUTION_spatial_significance.md for the full write-up.", drop the "Combining correlated
p-values (ewstools.pvalues) …" bullet, and fix the new-test count. spatial.py and
test_spatial.py need no change for the split — I checked the coupling: nothing imports
pvalues, and the one "pvalues" string in test_spatial.py is a test name.

(f) The missing-unit warning. Right now a single NaN cell nulls the whole time point silently,
and a gappy sensor field is the main real use of this class. Same guard in both
compute_moran and compute_moran_significance:

        values = df_pre.apply(lambda row: morans_i(row.to_numpy(), self.weights), axis=1)
        self.ews["morans_i"] = values
        n_nan = int(values.isna().sum())
        if n_nan:
            warnings.warn(
                f"Moran's I is nan for {n_nan} of {len(df_pre)} time points. A time "
                "point yields nan if any unit is missing at that time, or if the "
                "field has no variation there. Missing units are not dropped: the "
                "whole time point is discarded.",
                RuntimeWarning,
                stacklevel=2,
            )

(import warnings at the top of spatial.py; the same three lines in
compute_moran_significance, counting NaN of the p-value Series. Assigning the pre-transition
Series into self.ews and then counting NaN over that column counts the post-transition rows
too, which is why the count is taken on values.) One test exercises both:

def test_compute_moran_warns_on_rows_with_nan():
    rng = np.random.default_rng(0)
    data = pd.DataFrame(rng.normal(size=(6, 9)), index=np.arange(6))
    data.iloc[2, 4] = np.nan
    weights = _grid_rook_weights(3)

    spatial = SpatialEWS(data, weights)
    with pytest.warns(RuntimeWarning, match=r"nan for 1 of 6 time points"):
        spatial.compute_moran()
    assert spatial.ews["morans_i"].isna().sum() == 1

    with pytest.warns(RuntimeWarning, match=r"nan for 1 of 6 time points"):
        spatial.compute_moran_significance(n_permutations=20, seed=0)
    assert spatial.ews["morans_i_pvalue"].isna().sum() == 1

Also a docstring line on compute_moran: "Computed on the raw field; for trend control see the
recipe in SpatialEWS. Rows containing NaN produce NaN and are counted in a warning."

One courtesy line, and genuinely your call — nothing about it holds up the merge. Your > 0.5
assertion on the 3×3/T=30 fixture passes at seed 5 by 0.017, and fails on about half of seeds
(526 of 1000: mean 0.486, sd 0.096). At 6×6/T=60 a > 0.30 bound holds on 1000/1000. It's the
geometry and the threshold together rather than either alone — at 3×3 even > 0.3 fails about 5%
of seeds.

Follow-ups I'm taking, explicitly not merge conditions

• The trend-control recipe as a documented recipe, with a shipped null and positive control for it.
• A lattice_weights convenience for building a grid W, plus pointers to libpysal.weights for
  everything else. To be fair to your N4 point: an arbitrary weights matrix is the
  irregular-network generalisation — what's missing is a convenience builder, and the rest is a
  library's job, not this package's.
• A reference-period baseline for the spatial indicator, since the value depends on the unit set and
  the weights as much as on the field.
• On your point about reimplementing rather than wrapping: morans_i stays in the package under
  the scope rule because it is computed over time as a resilience indicator, which is what
  CONTRIBUTING's in-scope clause covers — the same reason variance and autocorrelation are in
  scope. Whether to take esda and libpysal as hard dependencies is a package-level decision, not
  one this PR should carry. So the reimplementation is accepted as it stands; the follow-up is an
  interop pointer to esda.moran.Moran for analytic inference on a single snapshot. And a new
  module rather than helpers.py is fine for a new family of indicators.

Say if you'd rather do any of these yourself — I'd rather hand them over than duplicate you.

If life intervenes

Say the word and I'll apply the checklist as separate commits on top of yours — your commit and
authorship stay exactly as they are. No hurry either way; nothing here is on a clock.

One last thing

Merging into main doesn't publish a release — it lands in main, and shipping is a separate
tagging decision. So there's no pressure from the merge itself.

Thanks for this. A first pull request that computes the right statistic, tests it against a
hand-computed answer, and cites its literature correctly is rarer than it should be.

[C95234]

Thanks for this — genuinely one of the most careful reviews I've had on anything, and the mutation testing especially is not something I'd have thought to do myself.

Applied the checklist in the latest commit:

• (a) p-value floor — done, with the docstring line.
• (b) Renamed the test to test_spatialews_ktau_rises_under_strengthening_gradient_without_trend_control, marked it explicitly as a negative control, and added the Trend control recipe + the weights caveat to the SpatialEWS docstring, both close to verbatim from your suggestion. While I was in there I also moved the fixture itself from 3x3/T=30 to 6x6/T=60 with the >0.30 bound, since you'd already done the work of validating it at 1000/1000 seeds — I reran it independently on another 500 seeds before committing, also 500/500. The assertion change wasn't required for the merge per your note, but leaving an acknowledged coin-flip in place once I knew about it didn't sit right.
• (c)/(d) ewstools/spatial.py now exported from __init__.py and added to docs/source/ewstools.rst.
• (f) NaN warning added to both compute_moran and compute_moran_significance, plus your test verbatim.
• The split — agreed, and done: ewstools/pvalues.py, tests/test_pvalues.py and CONTRIBUTION_spatial_significance.md are gone from this PR, CHANGELOG updated (test count too). I also fixed the two things you found in it (the anti-conservative negative-correlation bug — clipped to 0 now, verified it matches independent-Fisher exactly on a perfectly anti-correlated case — and the Kost & McDermott / Brown naming mix-up) in a copy I'm keeping outside this repo. Appreciate the SciPy docstring pointer — didn't know that gap was already acknowledged there.

Full suite: 46 passed, 1 skipped, 0 regressions on the 32 that existed before this PR.

Everything above I did myself rather than leaving it to you — happy to take a pass at any of the "follow-ups, not merge conditions" list too if useful, just say which.

🤖 Generated with Claude Code

[C95234]

Also went ahead and added the three follow-ups from your "not merge conditions" list, since they were small enough to fold in now rather than leave for later:

• `lattice_weights(n_rows, n_cols=None, connectivity="rook"|"queen")` — the grid convenience, with the `libpysal.weights` pointer for everything else.
• The `esda.moran.Moran` interop line in the `SpatialEWS` docstring.
• A reference-period note (didn't attempt an implementation of it — that one's still yours if you want it, this is just the docstring caveat).

Full suite: 50 passed, 1 skipped.

🤖 Generated with Claude Code

[energyscholar]

Checked the new head properly rather than on trust: 50 passed / 1 skipped locally, the
trend-control recipe runs verbatim (tau +0.68 -> -0.07 on the strengthening gradient), and the
new fixture holds on 1000/1000 seeds. Thank you for all of it, the fixture change and the
follow-ups included. Approving.

Three small things I found, none blocking. I'll carry them in my follow-up PR, so there is
nothing more for you to do here:

1. The p-value floor has no test (my checklist's omission, not yours).
2. esda's Moran row-standardises by default, so it gives a different I from morans_i on the
   same W; the pointer will say transformation='b'.
3. lattice_weights(-2) quietly returns a 4x4 zero matrix; the follow-up rejects it.

Next from me: the shipped null/positive control and the reference-period baseline, built on
your lattice_weights.