Fix NaN and inf in compute_dataset_statistics

docs/source/whats_new.rst

@@ -57,6 +57,7 @@ Requirements
 
 Bugs
 ~~~~
+- Fix :func:`moabb.analysis.meta_analysis.compute_dataset_statistics` returning NaN when two pipelines score identically on every subject of a dataset. :func:`moabb.analysis.meta_analysis.compute_pvals_wilcoxon` called ``scipy.stats.wilcoxon`` on paired differences that are all zero, where the test is undefined: older SciPy raised ``ValueError`` (the traceback in the report), and current SciPy returns 1.0 from its exact method but NaN from the normal approximation it switches to on larger samples, so the result silently depended on the number of subjects. The degenerate pair is now reported as a one-tailed p-value of 0.5 in both directions, which is what the exact method already yielded, and the branch keeps its p-values strictly inside ``(0, 1)`` the way the permutation branch does, so Stouffer's method in :func:`moabb.analysis.meta_analysis.combine_pvalues` cannot see an infinite z-score. :func:`moabb.analysis.meta_analysis.compute_effect` no longer divides by a zero standard deviation, and :func:`moabb.analysis.meta_analysis.find_significant_differences` now applies the NaN fallback that was written but left commented out, and no longer writes a bare ``NaN`` to stdout (:gh:`678` by `Aditya Singh`_)
 - Fix how the benchmark results page (:doc:`paper_results`) describes what its tables report. It stated that results are "mean accuracy and standard deviation across all folds for all sessions and subjects", and both halves are inaccurate: :class:`moabb.paradigms.MotorImagery` selects the metric from the number of classes, so two-class scenarios are scored with ROC-AUC rather than accuracy, and :class:`moabb.evaluations.WithinSessionEvaluation` averages the cross-validation folds within each session before returning a score, so the reported standard deviation is across (subject, session) pairs and not across individual folds (:gh:`1128` by `Bhargav Kowshik`_)
 - Fix datasets ignoring a change of download directory: datasets now inherit ``MNE_DATA`` without persisting a redundant per-dataset mirror, and :func:`moabb.utils.set_download_dir` removes legacy ``MNE_DATASETS_<SIGN>_PATH`` entries that still mirror the previous shared location while preserving explicit overrides. :class:`moabb.datasets.RomaniBF2025ERP` now uses the same path mechanism and honours ``path`` and ``force_update``. Adds isolated regression coverage across every dataset (:gh:`1115` by `Bruno Aristimunha`_).
 - Add a ``__repr__`` to :class:`moabb.datasets.base.BaseDataset` so datasets display by their code (e.g. ``BNCI2014-001``) when printed, instead of the verbose default ``<...object at 0x...>``. This declutters the output of ``print(paradigm.datasets)`` in the tutorials and of the paradigm and evaluation compatibility warnings (by `Danae`_)
@@ -923,4 +924,5 @@ API changes
 .. _Danae: https://github.com/dnplchrn
 .. _Henrique Lefundes: https://github.com/HenriqueLefundes
 .. _Paul-Adrien Graignic: https://github.com/pagraignic-yneuro
+.. _Aditya Singh: https://github.com/adityasingh2400
 .. _pre-commit-ci: https://github.com/apps/pre-commit-ci

moabb/analysis/meta_analysis.py

@@ -11,6 +11,13 @@
 
 log = logging.getLogger(__name__)
 
+# Stouffer's method, used by combine_pvalues, turns a p-value of exactly 0 or 1
+# into an infinite z-score. The permutation branch already keeps its p-values
+# strictly inside (0, 1); these are the tightest bounds that do the same for the
+# Wilcoxon branch without moving any p-value that is already valid.
+_P_FLOOR = np.nextafter(0.0, 1.0)
+_P_CEIL = np.nextafter(1.0, 0.0)
+
 
 def collapse_session_scores(df):
     """Prepare results dataframe for computing statistics.
@@ -64,13 +71,28 @@ def compute_pvals_wilcoxon(df, order=None):
             if i != j:
                 pipe1 = order[i]
                 pipe2 = order[j]
+                diffs = df.loc[:, pipe1] - df.loc[:, pipe2]
+                if (diffs == 0).all():
+                    # Wilcoxon is undefined when every paired difference is
+                    # zero. Old SciPy raised ValueError here; current SciPy
+                    # returns 1.0 from its exact method but NaN from the normal
+                    # approximation it switches to on larger samples, so the
+                    # result silently depended on the number of subjects. The
+                    # two pipelines are indistinguishable, so the one-tailed
+                    # p-value is 0.5 in both directions, which is what the exact
+                    # method already yields.
+                    out[i, j] = 0.5
+                    continue
                 p = stats.wilcoxon(df.loc[:, pipe1], df.loc[:, pipe2])[1]
                 p /= 2
                 # we want the one-tailed p-value
-                diff = (df.loc[:, pipe1] - df.loc[:, pipe2]).mean()
-                if diff < 0:
+                if diffs.mean() < 0:
                     p = 1 - p  # was in the other side of the distribution
-                out[i, j] = p
+                # Keep p strictly inside (0, 1) so Stouffer's method stays
+                # finite, as the permutation branch already does. The normal
+                # approximation can underflow to an exact 0, which the one-tailed
+                # flip then turns into an exact 1.
+                out[i, j] = min(max(p, _P_FLOOR), _P_CEIL)
     return out
 
 
@@ -217,8 +239,18 @@ def compute_effect(df, order=None):
             if i != j:
                 # for now it's just the standardized difference
                 diffs = df.loc[:, pipe1] - df.loc[:, pipe2]
-                diffs = diffs.mean() / diffs.std()
-                out[i, j] = diffs
+                mean, std = diffs.mean(), diffs.std()
+                if std == 0:
+                    # The paired differences have no spread, so the standardized
+                    # difference is 0/0 when the two pipelines score identically
+                    # and c/0 when they differ by a constant. Identical pipelines
+                    # have no effect, so report 0 rather than the NaN (plus
+                    # RuntimeWarning) that NumPy would produce. A constant offset
+                    # really is an unbounded effect, so keep the sign and make
+                    # the infinity deliberate rather than a division accident.
+                    out[i, j] = 0.0 if mean == 0 else np.sign(mean) * np.inf
+                else:
+                    out[i, j] = mean / std
     return out
 
 
@@ -345,9 +377,12 @@ def find_significant_differences(df, perm_cutoff=20):
                 t = T_full.loc[(slice(None), algs[i]), algs[j]]
                 P[i, j] = combine_pvalues(p, nsubs)
                 if np.isnan(P[i, j]):
+                    # A dataset can be missing this pipeline pair entirely, in
+                    # which case pivot_table leaves a NaN behind. Treat the
+                    # missing evidence as "no difference" instead of returning a
+                    # NaN that every downstream consumer has to special-case.
                     log.info("NaN p-value found, turned to 1")
-                    print("NaN")
-                    # P[i, j] = 1.0
+                    P[i, j] = 1.0
                 T[i, j] = combine_effects(t, nsubs)
     dfP = pd.DataFrame(index=algs, columns=algs, data=P)
     dfT = pd.DataFrame(index=algs, columns=algs, data=T)

moabb/tests/test_analysis.py

@@ -1,9 +1,12 @@
+import logging
 import os
 import shutil
+import warnings
 from pathlib import Path
 
 import numpy as np
 import pandas as pd
+import pytest
 from matplotlib.pyplot import Figure
 
 import moabb.analysis.meta_analysis as ma
@@ -166,6 +169,137 @@ def test_compute_pvals_random_cannot_be_zero(self):
         p1vsp2 = pvals[0, 1]
         assert p1vsp2 >= 1 / n_perms, f"P-values cannot be zero {pvals}"
 
+    @pytest.mark.parametrize("n_subjects", [22, 25, 40])
+    def test_wilcoxon_identical_pipelines(self, n_subjects):
+        # Wilcoxon is undefined when every paired difference is zero. SciPy
+        # returns 1.0 from its exact method but NaN from the normal
+        # approximation it uses on larger samples, so the answer used to depend
+        # on the number of subjects. See issue #678.
+        df = pd.DataFrame(
+            {"pipeline_1": [0.7] * n_subjects, "pipeline_2": [0.7] * n_subjects}
+        )
+        with warnings.catch_warnings():
+            warnings.simplefilter("error", RuntimeWarning)
+            pvals = ma.compute_pvals_wilcoxon(df)
+        assert np.isfinite(pvals).all(), f"P-values must be finite {pvals}"
+        assert pvals[0, 1] == 0.5, f"Indistinguishable pipelines give 0.5 {pvals}"
+        assert pvals[1, 0] == 0.5, f"Indistinguishable pipelines give 0.5 {pvals}"
+
+    def test_wilcoxon_stays_inside_unit_interval(self):
+        # Stouffer's method maps 0 and 1 to an infinite z-score, so the Wilcoxon
+        # branch must keep p strictly inside (0, 1), as the permutation branch
+        # already does.
+        rng = np.random.RandomState(0)
+        n = 60
+        base = rng.uniform(0.4, 0.9, size=n)
+        df = pd.DataFrame({"pipeline_1": base, "pipeline_2": base + 0.3})
+        pvals = ma.compute_pvals_wilcoxon(df)
+        offdiag = pvals[~np.eye(2, dtype=bool)]
+        assert np.all(offdiag > 0), f"P-values cannot be zero {pvals}"
+        assert np.all(offdiag < 1), f"P-values cannot be one {pvals}"
+
+    def test_compute_effect_zero_spread(self):
+        # Identical pipelines give 0/0 and a constant offset gives c/0. Neither
+        # should come back as NaN. See issue #678.
+        # 0.75 - 0.5 is exact in binary, so the paired differences really do
+        # have a standard deviation of zero rather than a rounding residue.
+        df = pd.DataFrame(
+            {
+                "pipeline_1": [0.5] * 10,
+                "pipeline_2": [0.5] * 10,
+                "pipeline_3": [0.75] * 10,
+            }
+        )
+        with warnings.catch_warnings():
+            warnings.simplefilter("error", RuntimeWarning)
+            effect = ma.compute_effect(df)
+        assert not np.isnan(effect).any(), f"Effect sizes must not be NaN {effect}"
+        assert effect[0, 1] == 0.0, f"Identical pipelines have no effect {effect}"
+        assert effect[1, 0] == 0.0, f"Identical pipelines have no effect {effect}"
+        assert effect[2, 0] == np.inf, f"A constant gain is unbounded {effect}"
+        assert effect[0, 2] == -np.inf, f"A constant loss is unbounded {effect}"
+
+    def test_dataset_statistics_no_nan_for_identical_pipelines(self):
+        # End-to-end check of the path reported in issue #678: two pipelines
+        # that score identically used to make compute_dataset_statistics emit
+        # NaN, and find_significant_differences then printed "NaN" to stdout and
+        # left the NaN in place.
+        n_subjects = 25  # above perm_cutoff, so the Wilcoxon branch is used
+        results = pd.DataFrame(
+            [
+                {"pipeline": pipeline, "dataset": "D1", "subject": subject, "score": 0.7}
+                for subject in range(n_subjects)
+                for pipeline in ("pipeline_1", "pipeline_2")
+            ]
+        )
+
+        stats_df = ma.compute_dataset_statistics(results)
+        assert not stats_df["p"].isna().any(), f"NaN p-value {stats_df}"
+        assert not stats_df["smd"].isna().any(), f"NaN effect size {stats_df}"
+
+        dfP, dfT = ma.find_significant_differences(stats_df)
+        offdiag = ~np.eye(len(dfP), dtype=bool)
+        assert np.isfinite(dfP.to_numpy(dtype=float)[offdiag]).all(), dfP
+        assert np.isfinite(dfT.to_numpy(dtype=float)[offdiag]).all(), dfT
+
+    @staticmethod
+    def _stats_with_missing_pair():
+        """Two datasets where only D1 ran pipeline_1 against pipeline_2.
+
+        ``pivot_table`` leaves a NaN in the (D2, pipeline_1) x pipeline_2 cell,
+        which ``combine_pvalues`` then turns into a NaN combined p-value.
+        """
+        pipelines = ["pipeline_1", "pipeline_2", "pipeline_3"]
+        rows = [
+            {
+                "dataset": "D1",
+                "pipe1": pipe1,
+                "pipe2": pipe2,
+                "p": 0.3,
+                "smd": 0.1,
+                "nsub": 10,
+            }
+            for pipe1 in pipelines
+            for pipe2 in pipelines
+            if pipe1 != pipe2
+        ]
+        for pipe in ("pipeline_1", "pipeline_2"):
+            rows.append(
+                {
+                    "dataset": "D2",
+                    "pipe1": pipe,
+                    "pipe2": "pipeline_3",
+                    "p": 0.4,
+                    "smd": 0.2,
+                    "nsub": 10,
+                }
+            )
+            rows.append(
+                {
+                    "dataset": "D2",
+                    "pipe1": "pipeline_3",
+                    "pipe2": pipe,
+                    "p": 0.6,
+                    "smd": -0.2,
+                    "nsub": 10,
+                }
+            )
+        return pd.DataFrame(rows)
+
+    def test_find_significant_differences_turns_nan_into_one(self, caplog):
+        # The NaN backstop was written but left commented out, with a bare
+        # print("NaN") standing in for it. See issue #678.
+        stats_df = self._stats_with_missing_pair()
+        with caplog.at_level(logging.INFO, logger="moabb.analysis.meta_analysis"):
+            dfP, _ = ma.find_significant_differences(stats_df)
+        assert dfP.loc["pipeline_1", "pipeline_2"] == 1.0, dfP
+        assert dfP.loc["pipeline_2", "pipeline_1"] == 1.0, dfP
+        assert "NaN p-value found, turned to 1" in caplog.text
+
+    def test_find_significant_differences_does_not_print(self, capsys):
+        ma.find_significant_differences(self._stats_with_missing_pair())
+        assert capsys.readouterr().out == "", "find_significant_differences must be quiet"
+
 
 class TestResults:
     def setup_method(self, method):