diff --git a/pymare/estimators/combination.py b/pymare/estimators/combination.py index 93f4af1..2f18731 100644 --- a/pymare/estimators/combination.py +++ b/pymare/estimators/combination.py @@ -97,7 +97,28 @@ 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. + + .. 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 +133,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/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()