diff --git a/docs/source/whats_new.rst b/docs/source/whats_new.rst index fddffd360..932762f39 100644 --- a/docs/source/whats_new.rst +++ b/docs/source/whats_new.rst @@ -35,7 +35,7 @@ Requirements Bugs ~~~~ -- None yet. +- Fix :func:`moabb.analysis.meta_analysis.compute_pvals_wilcoxon` reporting the wrong tail when the sign of the mean paired difference disagrees with the signed-rank statistic: the one-tailed p-value is now taken directly from ``scipy.stats.wilcoxon(..., alternative="greater")`` instead of halving the two-sided value and choosing the side from the mean (:gh:`1177` by `Azra Bano`_) Code health ~~~~~~~~~~~ diff --git a/moabb/analysis/meta_analysis.py b/moabb/analysis/meta_analysis.py index 6d3655f72..f34b4e840 100644 --- a/moabb/analysis/meta_analysis.py +++ b/moabb/analysis/meta_analysis.py @@ -147,21 +147,28 @@ def compute_lowest_subject_scores(df, reference_pipeline, percentile=20): def compute_pvals_wilcoxon(df, order=None): """Compute Wilcoxon rank-sum test on aggregated results. - Returns kxk matrix of p-values computed via the Wilcoxon rank-sum test, - order defines the order of rows and columns + Returns a square matrix of p-values computed via the Wilcoxon rank-sum test, + order defines the order of rows and columns. + + Entry ``[i, j]`` is the one-sided p-value for the hypothesis that pipeline + ``order[i]`` scores higher than pipeline ``order[j]`` (SciPy's + ``alternative="greater"``); the opposite direction is tested at ``[j, i]``. + The p-values are not corrected for multiple comparisons: with ``k`` + pipelines this function performs ``k * (k - 1)`` tests, so apply a + correction (for example Bonferroni or Holm) before interpreting them. Parameters ---------- df: :class:`pandas.DataFrame` Aggregated results, samples are index, columns are pipelines, and values are scores - order: list - list of length (num algorithms) with names corresponding to df columns + order: list of length (n_pipelines) + Names corresponding to df columns Returns ------- pvals: ndarray of shape (n_pipelines, n_pipelines) - array of pvalues + pvalues """ _validate_finite_scores(df) if order is None: @@ -188,11 +195,13 @@ def compute_pvals_wilcoxon(df, order=None): # 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 - if diffs.mean() < 0: - p = 1 - p # was in the other side of the distribution + # One-tailed p-value that pipe1 scores higher than pipe2. The + # direction of the signed-rank test is given by its rank sums, + # not by the sign of the mean difference, so ask SciPy for the + # one-sided test rather than halving the two-sided p-value. + p = stats.wilcoxon( + df.loc[:, pipe1], df.loc[:, pipe2], alternative="greater" + )[1] # 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 diff --git a/moabb/tests/test_analysis.py b/moabb/tests/test_analysis.py index 65ae566f6..9281cfabb 100644 --- a/moabb/tests/test_analysis.py +++ b/moabb/tests/test_analysis.py @@ -6,6 +6,7 @@ import numpy as np import pandas as pd import pytest +import scipy.stats as stats from matplotlib.pyplot import Figure import moabb.analysis.meta_analysis as ma @@ -201,6 +202,26 @@ def test_wilcoxon_identical_pipelines(self, n_subjects): 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_tail_follows_rank_sums_not_mean(self): + # pipeline_1 beats pipeline_2 on 6 of 7 subjects by a small margin and + # loses badly on the remaining one, so the mean difference is negative + # while the signed-rank statistic favours pipeline_1. The one-tailed + # p-value must follow the rank sums, not the sign of the mean. See + # issue #1176. + diffs = np.array([0.02, 0.02, 0.02, 0.02, 0.02, 0.02, -0.20]) + df = pd.DataFrame({"pipeline_1": 0.7 + diffs, "pipeline_2": [0.7] * 7}) + assert diffs.mean() < 0 + pvals = ma.compute_pvals_wilcoxon(df) + expected_greater = stats.wilcoxon( + df["pipeline_1"], df["pipeline_2"], alternative="greater" + )[1] + expected_less = stats.wilcoxon( + df["pipeline_1"], df["pipeline_2"], alternative="less" + )[1] + assert np.isclose(pvals[0, 1], expected_greater), pvals + assert np.isclose(pvals[1, 0], expected_less), pvals + assert pvals[0, 1] < 0.5 < pvals[1, 0], 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