From 6e60e380ee04f1af4acf6df9393399a8051bb78f Mon Sep 17 00:00:00 2001 From: James Kent Date: Thu, 20 Aug 2026 15:03:30 -0500 Subject: [PATCH 1/2] [FIX] concordant z is the statistic, not the inverse of its corrected p The "concordant" mode reported `norm.isf(p)` of a p-value that had already been doubled to correct for testing both tails. That is a tail mismatch rather than a convention: it shrank the statistic to absorb a multiplicity penalty that belongs to the p-value alone. Winkler et al. (2016), which the mode docstrings cite, define the concordant statistic as T = max(-2 sum ln p_k, -2 sum ln (1 - p_k)) -- the better of the two directed combinations, unshrunk -- and apply the two-test correction to its p-value. That is `norm.isf(p / 2)`, and it now reproduces the winning directed statistic exactly (4e-16 over 1650 cases across both estimators). Three symptoms had one cause. The reported z was negative wherever p > 0.5, carried the opposite sign to the effect there, and was -inf wherever the cap put p at exactly 1 -- so the least significant results had the largest magnitudes. A single input now combines to its own z; it was 2.7822 for an input of 3.0. Thresholding the old map at any z_c selected p < sf(z_c) instead of p < 2 sf(z_c), so every concordant z-map ran at exactly half its intended two-sided alpha. The fix only raises |z|, so at a fixed threshold it can add voxels but never remove them. `p` and `logp` are unchanged. The four pinned concordant values were not independent references: before 091b633 the tests computed `ss.norm.isf(p)` themselves and asserted it against a constant, and that commit moved the expression into the estimator and carried the constants across. They assert that norm.isf(p) equals norm.isf(p), and so could not detect this. This completes #125, which gave the concordant statistic its direction but left the magnitude alone. Adds a parametrized invariant that z and p describe the same tail in all three modes -- the check that would have caught it, and that the two right-tailed modes have always satisfied. Co-Authored-By: Claude Opus 5 (1M context) --- pymare/estimators/combination.py | 35 +++++++++- pymare/results.py | 15 +++++ pymare/tests/test_combination_tests.py | 88 +++++++++++++++++++++++--- 3 files changed, 126 insertions(+), 12 deletions(-) diff --git a/pymare/estimators/combination.py b/pymare/estimators/combination.py index 93f4af1..78bc0fd 100644 --- a/pymare/estimators/combination.py +++ b/pymare/estimators/combination.py @@ -97,7 +97,36 @@ def p_value(self, z, *args, **kwargs): return np.exp(self.log_p_value(z, *args, **kwargs)) def fit(self, z, *args, **kwargs): - """Fit the estimator to z-values.""" + """Fit the estimator to z-values. + + Notes + ----- + ``params_["z"]`` and ``params_["p"]`` describe the same tail, in every + mode: ``norm.sf(z) == p`` under ``"directed"`` and ``"undirected"``, + whose p-values are right-tailed, and ``2 * norm.sf(|z|) == p`` under + ``"concordant"``, whose p-value carries the correction for testing both + tails. + + .. versionchanged:: 0.0.11 + The ``"concordant"`` statistic was ``norm.isf(p)`` applied to that + already-doubled p-value, which is a tail mismatch rather than a + different convention: it shrank the statistic to absorb the + multiplicity penalty that :footcite:t:`winkler2016non` place on the + p-value alone. Their concordant statistic is + ``T = max(-2 sum ln p_k, -2 sum ln (1 - p_k))`` -- the better of the + two directed combinations, unshrunk -- which is what + ``norm.isf(p / 2)`` returns. The old form also made the reported + statistic non-monotone in the evidence: it was negative wherever + ``p > 0.5``, carried the opposite sign to the effect there, and was + ``-inf`` wherever the cap put ``p`` at exactly 1, which gave the + least significant results the largest magnitudes. A single input now + combines to its own z, as it should. + + References + ---------- + .. footbibliography:: + + """ # This resets the Estimator's dataset_ attribute. fit_dataset will overwrite if called. self.dataset_ = None @@ -112,8 +141,8 @@ def fit(self, z, *args, **kwargs): # correction for two tests is an added log(2), the cap a minimum # against log(1). log_p = np.minimum(0.0, np.log(2.0) + np.minimum(log_p1, log_p2)) - z_calc = -ndtri_exp(log_p) - z_calc[log_p2 < log_p1] *= -1 + z_calc = -ndtri_exp(log_p - np.log(2.0)) + z_calc = np.where(log_p2 < log_p1, -z_calc, z_calc) else: if self.mode == "undirected": z = np.abs(z) diff --git a/pymare/results.py b/pymare/results.py index b96767a..b89baf1 100644 --- a/pymare/results.py +++ b/pymare/results.py @@ -779,6 +779,21 @@ class CombinationTestResults: Array of right-tailed p-values. Default = None. logp : :obj:`numpy.ndarray`, optional Array of natural logarithms of the right-tailed p-values. Default = None. + + Notes + ----- + Any one of the three determines the other two, but not equally well, so the + container rebuilds from the most precise one it was given. ``logp`` is the + primitive: over a few hundred inputs a combined p-value falls below the + smallest positive double, where ``p`` flushes to zero and the ``z`` rebuilt + from that zero is ``+inf``. + + Rebuilding ``z`` assumes the p-value is right-tailed, as the parameter says. + A ``"concordant"`` p-value is not -- it carries the correction for having + tested both tails, and its statistic needs both ``norm.isf(p / 2)`` and the + direction of the tail that won, neither of which the p-value alone + determines. Pass ``z`` for that mode, as + :meth:`~pymare.estimators.combination.CombinationTest.summary` does. """ def __init__(self, estimator, dataset, z=None, p=None, logp=None): diff --git a/pymare/tests/test_combination_tests.py b/pymare/tests/test_combination_tests.py index 883e17e..eba53eb 100644 --- a/pymare/tests/test_combination_tests.py +++ b/pymare/tests/test_combination_tests.py @@ -13,16 +13,16 @@ _params = [ (StoufferCombinationTest, _z1, "directed", [4.69574]), (StoufferCombinationTest, _z1, "undirected", [4.87462819]), - (StoufferCombinationTest, _z1, "concordant", [4.55204117]), + (StoufferCombinationTest, _z1, "concordant", [4.69574275]), (StoufferCombinationTest, _z2, "directed", [4.69574275, -4.16803071]), (StoufferCombinationTest, _z2, "undirected", [4.87462819, 4.16803071]), - (StoufferCombinationTest, _z2, "concordant", [4.55204117, -4.00717817]), + (StoufferCombinationTest, _z2, "concordant", [4.69574275, -4.16803071]), (FisherCombinationTest, _z1, "directed", [5.22413541]), (FisherCombinationTest, _z1, "undirected", [5.27449962]), - (FisherCombinationTest, _z1, "concordant", [5.09434911]), + (FisherCombinationTest, _z1, "concordant", [5.22413541]), (FisherCombinationTest, _z2, "directed", [5.22413541, -3.30626405]), (FisherCombinationTest, _z2, "undirected", [5.27449962, 4.27572965]), - (FisherCombinationTest, _z2, "concordant", [5.09434911, -4.11869468]), + (FisherCombinationTest, _z2, "concordant", [5.22413541, -4.27572965]), ] @@ -238,10 +238,14 @@ def test_combination_permutation_rejects_undirected_mode(combination_estimator): def test_combination_permutation_survives_saturated_p_values(): - """Concordant p caps at 1, so its z is -inf and cannot be compared. - - Ranking on z made every permutation tie at -inf, which read as "more - extreme than nothing" and returned the smallest achievable p-value. + """Concordant p caps at 1, and the permutation test must still rank on it. + + Ranking on z made every permutation tie, which read as "more extreme than + nothing" and returned the smallest achievable p-value. The tie is no longer + at -inf -- a capped p now reads as a z of exactly zero -- but z is still the + wrong thing to rank on, because it is signed: a strongly negative result and + a strongly positive one sit at opposite ends of it while carrying the same + concordant evidence. The log p-value orders them together. """ # Perfectly balanced z: neither tail wins, so both directed p-values are # 0.5 and the doubled minimum is capped at exactly 1. @@ -249,7 +253,8 @@ def test_combination_permutation_survives_saturated_p_values(): result = FisherCombinationTest(mode="concordant").fit_dataset(Dataset(y=z)).summary() assert np.ravel(result.p)[0] == 1.0 - assert not np.isfinite(result.z).all() # the condition that used to break it + # No evidence reads as zero, not as an infinite magnitude. + assert np.ravel(result.z)[0] == 0.0 # Ranking on z gave 1 / n_perm here. This data is as far from significant # as it gets, so the permutation p-value should sit at the other end. @@ -384,6 +389,71 @@ def test_stouffer_does_not_underflow_on_many_moderate_z(): assert np.allclose(fitted["logp"], ss.norm.logsf(60.0)) +#: The tail each mode's p-value describes, as a function mapping the reported z +#: back onto it. ``"concordant"`` doubles because its p-value carries the +#: correction for having tested both tails. +TAIL_OF_MODE = { + "directed": lambda z: ss.norm.sf(z), + "undirected": lambda z: ss.norm.sf(z), + "concordant": lambda z: 2 * ss.norm.sf(np.abs(z)), +} + + +@pytest.mark.parametrize("Cls", [StoufferCombinationTest, FisherCombinationTest]) +@pytest.mark.parametrize("mode", ["directed", "undirected", "concordant"]) +def test_reported_z_and_p_describe_the_same_tail(Cls, mode): + """A statistic and a p-value that disagree about their own tail area. + + This is the invariant the concordant mode broke: it reported ``norm.isf(p)`` + of an already-doubled p-value, so thresholding on z and thresholding on p + selected different results. It is asserted for all three modes because one + generic ``norm.isf(p)`` used to serve all three, and it happened to be right + for the two whose p-values are right-tailed -- which is why the mismatch went + unnoticed. + """ + fitted = Cls(mode).fit(_z2).params_ + + assert np.allclose(TAIL_OF_MODE[mode](fitted["z"]), fitted["p"], rtol=1e-9) + + +@pytest.mark.parametrize("Cls", [StoufferCombinationTest, FisherCombinationTest]) +def test_concordant_statistic_is_the_better_directed_combination(Cls): + """Winkler's T: ``max`` of the two directed combinations, not a shrunk copy. + + :footcite:t:`winkler2016non` define the concordant statistic as + ``T = max(-2 sum ln p_k, -2 sum ln (1 - p_k))`` and apply the correction for + two tests to its p-value. So the concordant statistic must equal whichever + directed statistic won, exactly -- inverting the corrected p-value directly + returns something smaller. + + References + ---------- + .. footbibliography:: + """ + directed = Cls("directed") + positive = np.ravel(directed.fit(_z2).params_["z"]) + negative = np.ravel(directed.fit(-_z2).params_["z"]) + winner = np.where(positive >= negative, positive, -negative) + + concordant = np.ravel(Cls("concordant").fit(_z2).params_["z"]) + + assert np.allclose(concordant, winner, rtol=1e-12, atol=1e-14) + + +def test_concordant_statistic_of_one_input_is_that_input(): + """The plainest case the old form got wrong: k = 1 must be the identity. + + A single z of 3.0 has a two-sided p of 0.0027, which is exactly what the + concordant test reports -- so the statistic beside it has to be 3.0. It was + 2.7822, the normal deviate of 0.0027 read as a one-tailed probability. + """ + for zin in (-3.0, -0.3, 0.3, 3.0): + fitted = StoufferCombinationTest("concordant").fit(np.array([[zin]])).params_ + + assert np.allclose(np.ravel(fitted["z"])[0], zin) + assert np.allclose(np.ravel(fitted["p"])[0], 2 * ss.norm.sf(abs(zin))) + + def test_public_p_value_survives_the_move_to_log_space(): """Subclasses now implement log_p_value, but p_value stays part of the API.""" est = StoufferCombinationTest() From d0ce4a3738414d6f4ef8bef7250b51ec8637653f Mon Sep 17 00:00:00 2001 From: James Kent Date: Thu, 20 Aug 2026 15:08:16 -0500 Subject: [PATCH 2/2] delete extra text Co-authored-by: James Kent --- pymare/estimators/combination.py | 8 -------- pymare/results.py | 15 --------------- 2 files changed, 23 deletions(-) diff --git a/pymare/estimators/combination.py b/pymare/estimators/combination.py index 78bc0fd..2f18731 100644 --- a/pymare/estimators/combination.py +++ b/pymare/estimators/combination.py @@ -99,14 +99,6 @@ def p_value(self, z, *args, **kwargs): def fit(self, z, *args, **kwargs): """Fit the estimator to z-values. - Notes - ----- - ``params_["z"]`` and ``params_["p"]`` describe the same tail, in every - mode: ``norm.sf(z) == p`` under ``"directed"`` and ``"undirected"``, - whose p-values are right-tailed, and ``2 * norm.sf(|z|) == p`` under - ``"concordant"``, whose p-value carries the correction for testing both - tails. - .. versionchanged:: 0.0.11 The ``"concordant"`` statistic was ``norm.isf(p)`` applied to that already-doubled p-value, which is a tail mismatch rather than a diff --git a/pymare/results.py b/pymare/results.py index b89baf1..b96767a 100644 --- a/pymare/results.py +++ b/pymare/results.py @@ -779,21 +779,6 @@ class CombinationTestResults: Array of right-tailed p-values. Default = None. logp : :obj:`numpy.ndarray`, optional Array of natural logarithms of the right-tailed p-values. Default = None. - - Notes - ----- - Any one of the three determines the other two, but not equally well, so the - container rebuilds from the most precise one it was given. ``logp`` is the - primitive: over a few hundred inputs a combined p-value falls below the - smallest positive double, where ``p`` flushes to zero and the ``z`` rebuilt - from that zero is ``+inf``. - - Rebuilding ``z`` assumes the p-value is right-tailed, as the parameter says. - A ``"concordant"`` p-value is not -- it carries the correction for having - tested both tails, and its statistic needs both ``norm.isf(p / 2)`` and the - direction of the tail that won, neither of which the p-value alone - determines. Pass ``z`` for that mode, as - :meth:`~pymare.estimators.combination.CombinationTest.summary` does. """ def __init__(self, estimator, dataset, z=None, p=None, logp=None):