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

CHANGELOG.md

@@ -3,6 +3,25 @@
 All notable changes to `ewstools` are documented here.
 This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+
+### Added
+
+- **Spatial early warning signals** (`ewstools.spatial`): `morans_i()`,
+  `morans_i_permutation_test()`, and a new `SpatialEWS` class (mirroring the
+  `data`/`state`/`ews` conventions of `MultiTimeSeries`) for computing Moran's I
+  and its significance across space at every time point — the spatial branch of
+  the critical-slowing-down literature (Dakos et al., 2010; MacLaren, Aihara &
+  Masuda, 2025) had no equivalent in this package. See
+  `CONTRIBUTION_spatial_significance.md` for the full write-up.
+- **Combining correlated p-values** (`ewstools.pvalues`): `combine_pvalues_ebm()`,
+  an implementation of the Empirical Brown's Method (Poole et al., 2016) for
+  combining several early-warning indicators' p-values without assuming they are
+  independent, which Fisher's method requires but indicators computed from the
+  same underlying system rarely satisfy.
+- 21 new tests (`tests/test_spatial.py`, `tests/test_pvalues.py`), 0 regressions
+  on the existing suite.
+
 ## [2.1.3] — 2026-07-28
 
 Maintenance release. **No API changes and no behavioural changes** — existing code and

CONTRIBUTION_spatial_significance.md

@@ -0,0 +1,79 @@
+# Contribution: spatial early warning signals + combining correlated p-values
+
+## What was missing
+
+`ewstools` implements the temporal branch of the critical-slowing-down
+literature thoroughly — variance, autocorrelation, skew, kurtosis,
+spectral indicators, all tracked over a rolling window in time (`TimeSeries`,
+`MultiTimeSeries`). It has no equivalent for the **spatial** branch of the
+same literature: for systems observed at many locations at once (a
+lattice, a sensor network, an administrative territory grid...), spatial
+autocorrelation between neighbouring units is a separate, well-established
+early warning signal (Dakos et al., 2010, "Spatial correlation as leading
+indicator of catastrophic shifts", *Theoretical Ecology*; still an active
+area — see MacLaren, Aihara & Masuda, 2025, on generalising these methods
+from idealised regular lattices, where almost all validation work has been
+done, to irregular real-world networks). Searching `ewstools/core.py` and
+`ewstools/helpers.py` for "moran", "spatial", "fisher", "brown", or
+"combine_p" turns up nothing — confirmed by direct inspection of both
+files before writing any code here, not assumed.
+
+A related, separate gap: `ewstools` never combines its own indicators.
+If a user computes both a temporal indicator (say, variance) and a spatial
+one (Moran's I) on the same system, there is no tool in the package to ask
+"are both unusual at once?" — which needs Fisher's method at minimum, and
+a **correlation-aware** version of it if the two indicators (as is typical,
+since they are computed from the same underlying system) are not
+independent. Fisher's method, applied naively to correlated inputs, is
+biased. Confirmed absent by the same search.
+
+## What this contribution adds
+
+- **`ewstools/spatial.py`** — `morans_i()` and `morans_i_permutation_test()`
+  as standalone functions, plus a `SpatialEWS` class that follows the
+  existing `MultiTimeSeries` conventions (`data` → `state` → `ews`,
+  transition-aware, a `compute_ktau()` for trend-testing the resulting
+  indicator series exactly like the rest of the package). Computes Moran's
+  I — and its permutation-test significance — at every time point from a
+  DataFrame with one column per spatial unit and a fixed spatial weights
+  matrix.
+- **`ewstools/pvalues.py`** — `combine_pvalues_ebm()`, a real implementation
+  of the Empirical Brown's Method (Poole, Gibbs, Shmulevich, Bernard &
+  Knijnenburg, 2016, *Bioinformatics*), which extends Brown's (1975)
+  covariance-corrected version of Fisher's method by estimating the
+  correlation between tests empirically from data rather than requiring it
+  to be known or assumed.
+- 21 new tests (`tests/test_spatial.py`, `tests/test_pvalues.py`): hand-computed
+  Moran's I examples, known limiting cases (checkerboard → strong negative
+  autocorrelation, smooth gradient → strong positive autocorrelation,
+  constant field → undefined), the classical permutation-test property
+  `E[I] → -1/(N-1)` under the null, and — for the p-value combiner — that it
+  reduces to plain Fisher when the underlying data are independent, and
+  produces a *less* extreme combined p-value than naive Fisher when they are
+  positively correlated (the whole point of the correction), verified
+  against a direct evaluation of Brown's polynomial at a known correlation.
+  Full suite run locally: 52 passed, 1 skipped (0 regressions on the 32
+  tests that existed before this change).
+
+## Usage angles
+
+- **Research**: any critical-transition study with spatially resolved data
+  (ecology, epidemiology, social/territorial data, sensor networks, plasma
+  diagnostics) gains a spatial EWS tool in the same package already used for
+  the temporal one, instead of having to hand-roll Moran's I separately.
+- **Industry / infrastructure monitoring**: sensor networks (power grids,
+  structural health monitoring, industrial process control) already produce
+  exactly the "one column per unit, one row per time step" shape this class
+  expects.
+- **Public services**: territorial statistics agencies publishing
+  regularly-updated indicators across administrative units can use
+  `SpatialEWS` directly on their own panel data without depending on a
+  domain-specific tool.
+
+## Status
+
+Not yet submitted as a pull request to `ThomasMBury/ewstools`. See
+`/positionnement/ewstools` on the Hélios site for the live, honestly
+maintained status once submission happens — this repository is a local
+fork/branch pending that decision, and is not yet available even via a
+public fork URL.

ewstools/pvalues.py

@@ -0,0 +1,113 @@
+"""
+Combining correlated p-values.
+
+`ewstools` computes several early warning indicators from the same
+underlying data (variance, autocorrelation, skew, a spatial indicator...).
+A natural next question -- not currently answered anywhere in the
+package -- is: are two or more of these signals unusual AT THE SAME TIME?
+Fisher's (1925) classical method for combining independent p-values,
+
+    T = -2 * sum(ln p_i),  T ~ chi2(2k) under H0,
+
+is the obvious tool, but it assumes the p_i are independent. Indicators
+computed from the same underlying system are typically correlated, which
+biases Fisher's method (too liberal if the correlation is positive, too
+conservative if negative).
+
+The Empirical Brown's Method (Poole, Gibbs, Shmulevich, Bernard &
+Knijnenburg, 2016, "Combining dependent P-values with an empirical
+adaptation of Brown's method", Bioinformatics 32(17), i430-i436) fixes
+this: it keeps Fisher's statistic T, but instead of assuming the
+textbook chi2(2k) null distribution, it estimates the TRUE variance of T
+under the observed correlation structure directly from data (rather than
+requiring the covariance to be known analytically, as in Brown's 1975
+original method), then fits a rescaled chi-squared distribution to that
+corrected mean and variance (method of moments).
+"""
+from __future__ import annotations
+
+import numpy as np
+from scipy.stats import chi2
+
+
+def combine_pvalues_ebm(data, p_values) -> dict:
+    """Combine one-sided, possibly-correlated p-values with the Empirical
+    Brown's Method (Poole et al., 2016), building on Brown's (1975)
+    polynomial approximation for the covariance of `-2 ln(p)` terms.
+
+    Parameters
+    ----------
+    data : array-like, shape (n_observations, n_tests)
+        The raw data underlying each of the `n_tests` p-values (one
+        column per test) -- e.g. the k indicator series each p-value was
+        derived from. Used only to estimate the empirical correlation
+        between tests; this is what makes the method "empirical" rather
+        than requiring the covariance to be assumed or known exactly, as
+        in Brown's (1975) original method.
+    p_values : array-like, shape (n_tests,)
+        The one-sided p-value from each test (already computed
+        elsewhere, e.g. by a surrogate or permutation test). A small
+        p-value must consistently mean the same thing (e.g. "evidence of
+        an approaching transition") across all tests for the combination
+        to be meaningful.
+
+    Returns
+    -------
+    dict with keys:
+        p_combined : combined p-value.
+        chi2_fisher : Fisher's uncorrected statistic, -2*sum(ln p_i).
+        chi2_ebm : chi2_fisher rescaled by the empirical correction factor `c`.
+        df_ebm : effective degrees of freedom after correction (<= 2*n_tests
+            when p-values are positively correlated).
+        c : empirical scale factor (> 1 under positive correlation, which is
+            what makes the combined p-value LESS extreme than naive Fisher
+            when the underlying tests are positively correlated -- the
+            correction Fisher's method is missing).
+
+    References
+    ----------
+    Brown, M. B. (1975). "A method for combining non-independent,
+    one-sided tests of significance." Biometrics, 31(4), 987-992.
+    Poole, W., Gibbs, D. L., Shmulevich, I., Bernard, B., & Knijnenburg,
+    T. A. (2016). "Combining dependent P-values with an empirical
+    adaptation of Brown's method." Bioinformatics, 32(17), i430-i436.
+    """
+    data = np.asarray(data, dtype=float)
+    p_values = np.asarray(p_values, dtype=float)
+    k = len(p_values)
+
+    if data.ndim != 2 or data.shape[1] != k:
+        raise ValueError("data must have shape (n_observations, n_tests) with n_tests == len(p_values)")
+    if k < 2:
+        raise ValueError("need at least two p-values to combine")
+
+    clipped = np.clip(p_values, 1e-15, 1.0)
+    chi2_fisher = float(-2.0 * np.sum(np.log(clipped)))
+
+    corr = np.corrcoef(data, rowvar=False)
+
+    # Brown's (1975) polynomial approximation of Cov(-2 ln p_i, -2 ln p_j)
+    # as a function of the Pearson correlation rho between the underlying
+    # (approximately normal) variables i and j.
+    cov_sum = 0.0
+    for i in range(k):
+        for j in range(i + 1, k):
+            rho = corr[i, j]
+            cov_sum += rho * (3.263 + rho * (0.710 + rho * 0.027))
+
+    var_fisher = 4 * k + 2 * cov_sum
+    mean_fisher = 2 * k
+
+    # Method-of-moments fit: T ~ c * chi2(df), matching mean and variance.
+    c = var_fisher / (2 * mean_fisher)
+    df_ebm = 2 * mean_fisher**2 / var_fisher
+
+    p_combined = float(chi2.sf(chi2_fisher / c, df_ebm))
+
+    return {
+        "p_combined": p_combined,
+        "chi2_fisher": chi2_fisher,
+        "chi2_ebm": chi2_fisher / c,
+        "df_ebm": float(df_ebm),
+        "c": float(c),
+    }

ewstools/spatial.py

@@ -0,0 +1,170 @@
+"""
+Spatial early warning signals.
+
+Most of `ewstools` treats early warning signals (EWS) as a purely
+temporal phenomenon: a rolling window slides along a single time series
+(`TimeSeries`) or a handful of co-measured series (`MultiTimeSeries`),
+and indicators like variance or lag-1 autocorrelation are tracked over
+time. But a large and separate branch of the critical-slowing-down
+literature -- starting with Dakos et al. (2010), "Spatial correlation as
+leading indicator of catastrophic shifts", and still an active area (see
+e.g. MacLaren, Aihara & Masuda, 2025, on generalising spatial EWS from
+regular lattices to irregular real-world networks) -- looks instead
+across SPACE at a single instant: do neighbouring units of a system
+become more correlated with each other as a critical transition
+approaches? Moran's I is the standard statistic for this, but it is not
+implemented anywhere in `ewstools` today.
+
+This module adds that missing piece, following the existing package's
+conventions (`ewstools.core.MultiTimeSeries`): a spatial analogue that
+takes one column per spatial unit and computes a spatial indicator at
+EVERY time point (cross-sectional at each row), rather than a temporal
+indicator over a rolling window.
+"""
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+
+
+def morans_i(values, weights) -> float:
+    """Moran's I spatial autocorrelation statistic.
+
+    I = (N / S0) * [sum_ij w_ij (x_i - xbar)(x_j - xbar)] / [sum_i (x_i - xbar)^2]
+
+    Parameters
+    ----------
+    values : array-like, shape (n_units,)
+        Observed value at each spatial unit.
+    weights : array-like, shape (n_units, n_units)
+        Spatial weights matrix, w_ij > 0 if units i and j are neighbours,
+        0 otherwise (zero diagonal). Not required to be row-standardised.
+
+    Returns
+    -------
+    float
+        Moran's I, or nan if the weights or the values have no variation
+        (undefined in that case).
+    """
+    x = np.asarray(values, dtype=float)
+    w = np.asarray(weights, dtype=float)
+    n = len(x)
+    if w.shape != (n, n):
+        raise ValueError(f"weights must be an {n} x {n} matrix to match values")
+
+    deviations = x - x.mean()
+    s0 = w.sum()
+    if s0 == 0:
+        return float("nan")
+
+    denominator = (deviations**2).sum()
+    if denominator == 0:
+        return float("nan")
+
+    numerator = deviations @ w @ deviations
+    return float((n / s0) * (numerator / denominator))
+
+
+def morans_i_permutation_test(values, weights, n_permutations: int = 500, seed=None) -> dict:
+    """Significance of Moran's I by permutation (spatial analogue of a
+    surrogate test: reshuffle the values across spatial units while
+    keeping the network itself fixed, following the standard approach
+    for testing spatial autocorrelation, e.g. Dakos et al., 2010).
+
+    Returns
+    -------
+    dict with keys: observed_i, p_value, n_permutations, null_mean, null_std.
+    `p_value` is one-sided (probability of a permuted I at least as large
+    as the observed one) -- appropriate for testing an INCREASE in
+    spatial correlation as an early warning signal.
+    """
+    observed = morans_i(values, weights)
+    if np.isnan(observed):
+        return {"observed_i": None, "p_value": None, "n_permutations": 0, "null_mean": None, "null_std": None}
+
+    rng = np.random.default_rng(seed)
+    values = np.asarray(values, dtype=float)
+    null_values = np.empty(n_permutations)
+    for k in range(n_permutations):
+        null_values[k] = morans_i(rng.permutation(values), weights)
+
+    return {
+        "observed_i": observed,
+        "p_value": float(np.mean(null_values >= observed)),
+        "n_permutations": n_permutations,
+        "null_mean": float(np.mean(null_values)),
+        "null_std": float(np.std(null_values)),
+    }
+
+
+class SpatialEWS:
+    """
+    Spatially-resolved data on which to compute spatial early warning
+    signals, following the `data`/`state`/`ews` conventions of
+    `ewstools.core.MultiTimeSeries`.
+
+    Parameters
+    ----------
+    data : pandas.DataFrame
+        One column per spatial unit (grid cell, node, sensor, region...),
+        one row per time point. Index represents time and is carried
+        over.
+    weights : array-like, shape (n_units, n_units)
+        Spatial weights matrix for the units in `data.columns` (same
+        order), fixed over time -- the network itself is not assumed to
+        change, only the values observed on it.
+    transition : float, optional
+        Time value at which a transition occurs, if any. If given,
+        spatial EWS are only computed up to this point.
+    """
+
+    def __init__(self, data, weights, transition=None):
+        if not isinstance(data, pd.DataFrame):
+            raise TypeError("data must be a pandas DataFrame (one column per spatial unit)")
+        weights = np.asarray(weights, dtype=float)
+        n_units = data.shape[1]
+        if weights.shape != (n_units, n_units):
+            raise ValueError(f"weights must be an N x N matrix with N = data.shape[1] = {n_units}")
+
+        self.state = data
+        self.weights = weights
+        self.transition = float(transition) if transition else transition
+        self.var_names = data.columns
+        self.ews = pd.DataFrame(index=data.index)
+        self.ktau = dict()
+
+    def _pre_transition(self) -> pd.DataFrame:
+        if self.transition:
+            return self.state[self.state.index <= self.transition]
+        return self.state
+
+    def compute_moran(self):
+        """Compute Moran's I at every time point. Output stored in
+        `self.ews['morans_i']`.
+        """
+        df_pre = self._pre_transition()
+        self.ews["morans_i"] = df_pre.apply(lambda row: morans_i(row.to_numpy(), self.weights), axis=1)
+
+    def compute_moran_significance(self, n_permutations: int = 500, seed=None):
+        """Permutation-test p-value for Moran's I at every time point.
+        Output stored in `self.ews['morans_i_pvalue']`.
+        """
+        df_pre = self._pre_transition()
+        self.ews["morans_i_pvalue"] = df_pre.apply(
+            lambda row: morans_i_permutation_test(row.to_numpy(), self.weights, n_permutations, seed)["p_value"],
+            axis=1,
+        )
+
+    def compute_ktau(self, tmin="earliest", tmax="latest"):
+        """Kendall tau of each spatial EWS against time -- same convention
+        as `ewstools.core.TimeSeries.compute_ktau`. Output stored in the
+        `self.ktau` dict.
+        """
+        if tmin == "earliest":
+            tmin = self.ews.dropna(how="all").index[0]
+        if tmax == "latest":
+            tmax = self.ews.dropna(how="all").index[-1]
+
+        df_ews = self.ews[(self.ews.index >= tmin) & (self.ews.index <= tmax)].copy()
+        time_values = pd.Series(data=df_ews.index, index=df_ews.index)
+        self.ktau = dict(df_ews.corrwith(time_values, method="kendall", axis=0))