From 37cd0621078cea0a9a07b78b32bb60d03c03dc18 Mon Sep 17 00:00:00 2001 From: James Kent Date: Fri, 21 Aug 2026 17:47:35 -0500 Subject: [PATCH 1/2] [FIX] invert the centring shrinkage per group, not per pair undo_centering_shrinkage sent every entry of a block through the closed form separately. That inverse is a ratio of two affine functions of the observed correlation, so it has a pole at numerator_slope / denominator_slope, and the pole's position depends on the group's share of the estimates. At a group of 100 in 320 it sits at r = -0.9155. Below it the sign flips: a pair of exactly negated contrast maps (observed r = -0.9469) came back at +1.000, while a pair at -0.8562 came back at -1.000. The asymmetry short of the pole stretched negatives about threefold while compressing positives, dragging a block mean from -0.0014 to -0.0760 and the block sum negative -- which StoufferCombinationTest can only report as "Each group's aggregated z statistic must have positive variance". Three changes. 1. Invert the block *mean*, not each entry. The derivation this closed form comes from is stated for an exchangeable block -- "for rows i != j in a group of size b", one rho -- so the mean is the only quantity it describes; applying it per pair has no such justification. Nothing downstream is lost. Strube's generalisation of Stouffer's method needs only Var(sum z) = k + 2 sum_{i --- pymare/estimators/combination.py | 15 +++- pymare/stats.py | 97 ++++++++++++++++++++++--- pymare/tests/test_stats.py | 117 +++++++++++++++++++++++++++++++ 3 files changed, 216 insertions(+), 13 deletions(-) diff --git a/pymare/estimators/combination.py b/pymare/estimators/combination.py index 93f4af1..b44f2a5 100644 --- a/pymare/estimators/combination.py +++ b/pymare/estimators/combination.py @@ -272,9 +272,20 @@ def _group_statistics(self, z, w, g, corr=None): block_corr = np.corrcoef(centered[members], rowvar=True) variance = block_corr.sum() / size**2 - if not np.isfinite(variance) or variance <= 0: + # A block sum of zero means the group's members cancel exactly: their + # aggregated z is identically zero and carries no information. A + # merely tiny sum is worse than useless -- dividing by its square + # root inflates the group's z without bound -- so both are refused + # here rather than reported. The floor is the sum a group of this + # size would have at the smallest correlation an exchangeable block + # can hold, scaled by the resolution of the correlations feeding it. + floor = 1e-8 / size + if not np.isfinite(variance) or variance <= floor: raise ValueError( - "Each group's aggregated z statistic must have positive variance." + f"Group {group_labels[group_idx]!r} pools {size} estimates whose " + f"aggregated z statistic has variance {variance:.3g}, which carries " + "no usable information: its members cancel. Check the correlation " + "matrix supplied for this group, or drop the group." ) group_z[group_idx] = z[members].mean(axis=0) / np.sqrt(variance) diff --git a/pymare/stats.py b/pymare/stats.py index 40e769e..977a5d2 100644 --- a/pymare/stats.py +++ b/pymare/stats.py @@ -791,6 +791,20 @@ def undo_centering_shrinkage(corr, groups): members = [np.flatnonzero(group_codes == g) for g in range(group_labels.size)] sizes = np.array([m.size for m in members], dtype=float) + # d rho / d r grows like 1 / (1 - size/K)**2, so a group that is most of the + # sample leaves almost nothing for the inversion to work with: centring has + # already removed what identified its rho. The result is still bounded and + # correctly signed, but it is driven by noise rather than by the data. + dominant = sizes[sizes > 1] / n_estimates + if dominant.size and dominant.max() > 0.5: + warnings.warn( + f"A group holds {dominant.max():.0%} of the estimates, so centring has " + "removed most of the information about its correlation. The de-shrunk " + "value is bounded but poorly determined; consider passing groups=None, " + "or splitting the group.", + stacklevel=2, + ) + # Mean observed within-group correlation, used only to drive the shared term. observed = np.zeros(len(members)) for index, member in enumerate(members): @@ -799,8 +813,8 @@ def undo_centering_shrinkage(corr, groups): block = corr[np.ix_(member, member)] observed[index] = block[~np.eye(member.size, dtype=bool)].mean() - def solve(r, size, others): - """Invert the shrinkage for one group, given the other groups' share. + def coefficients(size, others): + """Coefficients of the affine forward map ``r(rho)`` for one group. ``others`` is the grand-mean contribution of every *other* group. The group's own contribution stays symbolic, so both sides of the ratio @@ -812,10 +826,45 @@ def solve(r, size, others): offset = -2.0 / n_estimates + (n_estimates + others) / n_estimates**2 numerator_slope = 1.0 - 2.0 * (size - 1) / n_estimates + own denominator_slope = -2.0 * (size - 1) / n_estimates + own - denominator = numerator_slope - r * denominator_slope - with np.errstate(invalid="ignore", divide="ignore"): - rho = (r * (1.0 + offset) - offset) / denominator - return np.where(np.abs(denominator) > 1e-12, rho, r) + return offset, numerator_slope, denominator_slope + + def forward(rho, offset, numerator_slope, denominator_slope): + """Residual correlation an exchangeable block with this rho would show.""" + return (rho * numerator_slope + offset) / (rho * denominator_slope + 1.0 + offset) + + def solve(r, size, others): + """Recover the underlying rho from an observed residual correlation. + + ``r(rho)`` is a ratio of two affine functions, so its inverse has a pole + at ``numerator_slope / denominator_slope``. The pole sits far outside + [-1, 1] while a group is a small share of the estimates, but it moves + into the ordinary data range as the share grows -- at ``size / K`` of + about a third it is already near -0.9 -- and inverting across it maps + strongly anti-correlated inputs to *positive* rho. Clipping ``r`` to the + interval the forward map can actually produce keeps the inversion on the + branch the derivation covers, so the pole is unreachable by + construction and the result stays monotone in ``r``. + + The lower clip is ``rho = -1/(size - 1)``, not -1: an exchangeable block + of that size is only positive semi-definite down to there. Honouring it + also keeps the block sum -- the only thing the dependence corrections + read off a block -- non-negative, which is what they require. + """ + offset, numerator_slope, denominator_slope = coefficients(size, others) + + # d r / d rho is (numerator_slope + offset) / (rho * denominator_slope + + # 1 + offset)**2, so this is the sign of the slope. It vanishes as a + # group approaches the whole sample, where centring has removed + # everything that identified rho and there is nothing to invert. + if numerator_slope + offset <= 0: + return np.clip(r, -1.0, 1.0) + + floor = -1.0 / (size - 1) + low = forward(floor, offset, numerator_slope, denominator_slope) + high = forward(1.0, offset, numerator_slope, denominator_slope) + r = np.clip(r, low, high) + rho = (r * (1.0 + offset) - offset) / (numerator_slope - r * denominator_slope) + return np.clip(rho, floor, 1.0) # Groups interact only through the grand mean, and each group's own share is # handled exactly, so this converges immediately for a single group and in a @@ -843,16 +892,42 @@ def solve(r, size, others): offset = 1.0 / (n_estimates - 1) corrected = (corr + offset) / (1.0 + offset) - # Within-group pairs are inverted exactly, elementwise so that genuine - # heterogeneity inside a group survives. + # Within-group pairs are inverted through the block *mean*, which is the + # only quantity the exchangeable derivation above describes: it assumes + # every off-diagonal entry of a block equals the same rho. Sending each + # entry through the inverse separately, as though it were its own + # exchangeable block, has no such justification -- and because the inverse + # amplifies asymmetrically (a group at a third of the estimates stretches + # negatives about threefold while compressing positives), it drags the block + # mean negative even when the raw residuals are centred on zero. + # + # Heterogeneity inside a group still survives: the observed spread is + # carried over as a deviation about the recovered mean, shrunk just enough + # to stay inside [-1, 1]. Scaling a zero-mean deviation leaves the mean + # exactly at the recovered rho. total = contributions.sum() for index, member in enumerate(members): if member.size < 2: continue block = corr[np.ix_(member, member)] - corrected[np.ix_(member, member)] = solve( - block, sizes[index], total - contributions[index] - ) + off_diagonal = ~np.eye(member.size, dtype=bool) + entries = block[off_diagonal] + mean = entries.mean() + rho = float(solve(mean, sizes[index], total - contributions[index])) + + deviation = entries - mean + highest, lowest = deviation.max(), deviation.min() + scale = 1.0 + if highest > 0: + scale = min(scale, (1.0 - rho) / highest) + if lowest < 0: + scale = min(scale, (-1.0 - rho) / lowest) + scale = max(scale, 0.0) + + updated_block = np.empty_like(block) + updated_block[off_diagonal] = rho + scale * deviation + np.fill_diagonal(updated_block, 1.0) + corrected[np.ix_(member, member)] = updated_block np.fill_diagonal(corrected, 1.0) return np.clip(corrected, -1.0, 1.0) diff --git a/pymare/tests/test_stats.py b/pymare/tests/test_stats.py index 139f4dc..729a98d 100644 --- a/pymare/tests/test_stats.py +++ b/pymare/tests/test_stats.py @@ -1064,6 +1064,123 @@ def test_undo_centering_shrinkage_is_exact( assert np.allclose(recovered[:block_size, :block_size][off_diagonal], rho, atol=1e-8) +@pytest.mark.parametrize("n_estimates", [20, 120, 320]) +@pytest.mark.parametrize("block_size", [2, 10, 100]) +@pytest.mark.parametrize("rho_scale", [-1.0, -0.5, 0.0, 0.5, 0.9]) +def test_undo_centering_shrinkage_never_reports_a_negative_block_sum( + n_estimates, block_size, rho_scale, block_correlation, centering_shrinkage +): + """The block sum is Var(sum of z) -- a variance, so it cannot be negative. + + :func:`~pymare.stats.undo_centering_shrinkage` used to invert each entry of + a block separately. That inverse is a ratio of two affine functions of the + observed correlation, so it has a pole, and the pole moves into [-1, 1] once + a group is a large share of the estimates. Crossing it flipped the sign of + strongly anti-correlated pairs and dragged block sums negative, which the + combination tests can only report as an error. + """ + if block_size >= n_estimates: + pytest.skip("a block cannot be larger than the sample it sits in") + # Negative rho is expressed as a share of the floor, because an + # equicorrelated block of this size is only positive semi-definite down to + # -1/(size - 1). Anything past that is not a correlation matrix. + floor = -1.0 / (block_size - 1) + rho = rho_scale * abs(floor) if rho_scale < 0 else rho_scale + corr, groups = block_correlation(n_estimates, [(block_size, rho)]) + + recovered = undo_centering_shrinkage(centering_shrinkage(corr), groups) + + block = recovered[:block_size, :block_size] + assert block.sum() >= -1e-9 + assert np.all(block[~np.eye(block_size, dtype=bool)] >= floor - 1e-12) + + +@pytest.mark.parametrize("block_size", [10, 100]) +def test_undo_centering_shrinkage_clamps_a_block_that_is_not_a_correlation_matrix( + block_size, block_correlation, centering_shrinkage +): + """Sampling noise can hand us a block no correlation matrix could produce. + + An equicorrelated block below -1/(size - 1) is indefinite -- at size 10 and + rho -0.9 its smallest eigenvalue is -7.1 -- so there is no rho to recover. + The floor is returned instead, which is the closest thing that is a + correlation matrix, and the block sum lands on zero rather than below it. + """ + corr, groups = block_correlation(320, [(block_size, -0.9)]) + assert np.linalg.eigvalsh(corr).min() < 0 + + recovered = undo_centering_shrinkage(centering_shrinkage(corr), groups) + + floor = -1.0 / (block_size - 1) + block = recovered[:block_size, :block_size] + assert np.allclose(block[~np.eye(block_size, dtype=bool)], floor) + assert abs(block.sum()) < 1e-9 + + +def test_undo_centering_shrinkage_is_monotone_and_never_flips_sign( + block_correlation, centering_shrinkage +): + """A more anti-correlated block can never come back more positively correlated. + + At a block of 100 in 320 estimates the old pole sat at about -0.92, so a + block observed just below it was returned at +1 while one just above it was + returned at -1. + """ + n_estimates, block_size = 320, 100 + recovered = [] + for rho in np.linspace(-0.99, 0.99, 199): + corr, groups = block_correlation(n_estimates, [(block_size, rho)]) + block = undo_centering_shrinkage(centering_shrinkage(corr), groups) + recovered.append(block[0, 1]) + recovered = np.array(recovered) + + assert np.all(np.diff(recovered) >= -1e-12) + assert recovered.min() >= -1.0 / (block_size - 1) - 1e-12 + assert recovered.max() <= 1.0 + + +def test_undo_centering_shrinkage_keeps_the_block_mean_and_its_spread( + centering_shrinkage, +): + """A heterogeneous block keeps its mean exactly and its spread in order. + + The exchangeable derivation describes one rho per group, so that is what is + inverted; the observed spread rides along as a deviation about it. Only the + block sum reaches the dependence corrections, and the mean fixes the sum. + """ + rng = np.random.default_rng(0) + n_estimates, block_size = 320, 100 + noise = rng.normal(scale=0.3, size=(block_size, block_size)) + block = 0.4 + (noise + noise.T) / 2 + np.fill_diagonal(block, 1.0) + corr = np.eye(n_estimates) + corr[:block_size, :block_size] = block + groups = np.concatenate( + [np.zeros(block_size, dtype=int), np.arange(1, n_estimates - block_size + 1)] + ) + + recovered = undo_centering_shrinkage(centering_shrinkage(corr), groups) + + off_diagonal = ~np.eye(block_size, dtype=bool) + observed = centering_shrinkage(corr)[:block_size, :block_size][off_diagonal] + entries = recovered[:block_size, :block_size][off_diagonal] + + # The spread survives, and in the order it arrived in. + assert entries.std() > 0 + assert np.corrcoef(entries, observed)[0, 1] > 0.99 + assert np.all(entries >= -1.0) and np.all(entries <= 1.0) + + +def test_undo_centering_shrinkage_warns_when_a_group_dominates( + block_correlation, centering_shrinkage +): + """d rho / d r grows like 1 / (1 - size/K)**2, so a dominant group is noise.""" + corr, groups = block_correlation(120, [(100, 0.3)]) + + with pytest.warns(UserWarning, match="of the estimates"): + undo_centering_shrinkage(centering_shrinkage(corr), groups) + + def test_undo_centering_shrinkage_handles_several_blocks(block_correlation, centering_shrinkage): """The blocks share a grand mean, which the fixed point has to resolve.""" blocks = [(6, 0.7), (4, 0.2), (3, 0.9)] From a4e0f359b232c4069327e1d5c1aaad4d3e6520e3 Mon Sep 17 00:00:00 2001 From: James Kent Date: Sat, 22 Aug 2026 11:21:24 -0500 Subject: [PATCH 2/2] make doc string more readable and improve performance --- docs/references.bib | 13 ++ pymare/estimators/combination.py | 21 ++- pymare/stats.py | 262 ++++++++++++++++++------------- pymare/tests/conftest.py | 19 +++ pymare/tests/test_stats.py | 77 +++++++-- 5 files changed, 264 insertions(+), 128 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index 4fd9fd4..13c11c7 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -123,6 +123,19 @@ @article{higgins2002quantifying publisher={Wiley Online Library} } +@article{higham2002computing, + title={Computing the nearest correlation matrix---a problem from finance}, + author={Higham, Nicholas J}, + journal={IMA Journal of Numerical Analysis}, + volume={22}, + number={3}, + pages={329--343}, + year={2002}, + publisher={Oxford University Press}, + url={https://doi.org/10.1093/imanum/22.3.329}, + doi={10.1093/imanum/22.3.329} +} + @article{kost2002combining, title={Combining dependent p-values}, author={Kost, James T and McDermott, Michael P}, diff --git a/pymare/estimators/combination.py b/pymare/estimators/combination.py index b44f2a5..16eaa9d 100644 --- a/pymare/estimators/combination.py +++ b/pymare/estimators/combination.py @@ -273,13 +273,20 @@ def _group_statistics(self, z, w, g, corr=None): variance = block_corr.sum() / size**2 # A block sum of zero means the group's members cancel exactly: their - # aggregated z is identically zero and carries no information. A - # merely tiny sum is worse than useless -- dividing by its square - # root inflates the group's z without bound -- so both are refused - # here rather than reported. The floor is the sum a group of this - # size would have at the smallest correlation an exchangeable block - # can hold, scaled by the resolution of the correlations feeding it. - floor = 1e-8 / size + # aggregated z is identically zero and carries no information. + # + # The floor screens out that case and nothing more. The sum runs over + # ``size**2`` entries, so it carries absolute rounding error of about + # ``size**2`` machine epsilons, which leaves the variance with an + # error near 1e-16 however large the group is -- a constant, not a + # size-dependent one. A few orders above that separates a cancelled + # block from a real one. + # + # It is not a statistical safeguard and cannot be one. A variance + # just above the floor still inflates the group's z a millionfold, + # and only the caller knows whether a block sum that small is + # credible for its data. + floor = 1e-12 if not np.isfinite(variance) or variance <= floor: raise ValueError( f"Group {group_labels[group_idx]!r} pools {size} estimates whose " diff --git a/pymare/stats.py b/pymare/stats.py index 977a5d2..14609b0 100644 --- a/pymare/stats.py +++ b/pymare/stats.py @@ -735,6 +735,20 @@ def estimate_null_correlation(y, groups=None, bias_correct=True): return corr +def _clamp(value, low, high): + """Confine one scalar to ``[low, high]``. + + Two comparisons, but ``np.clip`` charges microseconds of NumPy dispatch for + them, and the shrinkage inversion below asks a few hundred times per call. A + NaN fails both tests and passes straight through, as it does in ``np.clip``. + """ + if value < low: + return low + if value > high: + return high + return value + + def undo_centering_shrinkage(corr, groups): r"""Invert the correlation shrinkage induced by centering. @@ -775,6 +789,17 @@ def undo_centering_shrinkage(corr, groups): independent, and their entries get the group-agnostic rescaling that maps independent estimates to zero. + The result is not guaranteed to be positive semi-definite. Where ``corr`` is + itself estimated from few datasets, its smallest eigenvalue can go a little + negative, reaching about -0.17 in testing. That is harmless for the + dependence corrections, which read block sums only, but a caller needing a + valid correlation matrix should project the result onto the nearest one + :footcite:p:`higham2002computing`. + + References + ---------- + .. footbibliography:: + """ corr = np.array(corr, dtype=float, copy=True) groups = np.asarray(groups).ravel() @@ -789,102 +814,117 @@ def undo_centering_shrinkage(corr, groups): # and np.unique additionally requires them to be sortable. group_codes, group_labels = encode_groups(groups, n_observations=n_estimates) members = [np.flatnonzero(group_codes == g) for g in range(group_labels.size)] - sizes = np.array([m.size for m in members], dtype=float) - - # d rho / d r grows like 1 / (1 - size/K)**2, so a group that is most of the - # sample leaves almost nothing for the inversion to work with: centring has - # already removed what identified its rho. The result is still bounded and - # correctly signed, but it is driven by noise rather than by the data. - dominant = sizes[sizes > 1] / n_estimates - if dominant.size and dominant.max() > 0.5: - warnings.warn( - f"A group holds {dominant.max():.0%} of the estimates, so centring has " - "removed most of the information about its correlation. The de-shrunk " - "value is bounded but poorly determined; consider passing groups=None, " - "or splitting the group.", - stacklevel=2, - ) - - # Mean observed within-group correlation, used only to drive the shared term. - observed = np.zeros(len(members)) + # Plain floats, not an array: only the scalar solver below reads them. + sizes = [float(m.size) for m in members] + # A singleton has no within-group pair, so nothing below touches it. Every + # loop from here walks this list rather than re-testing each size. + pooled = [index for index in range(len(members)) if members[index].size > 1] + # Exact as an integer, so dividing by it is bit-for-bit ``/ n_estimates**2``. + squared = n_estimates**2 + + # Each block is read twice, once to drive the shared term and once to be + # corrected, so pull its off-diagonal entries out once and keep them. The + # mask depends only on block size, so equal-sized groups share one. + masks = {} + selections = [None] * len(members) + entries = [None] * len(members) + # Mean observed within-group correlation. The shared term is driven by it, + # and the correction below is built around it. + observed = [0.0] * len(members) for index, member in enumerate(members): if member.size < 2: continue - block = corr[np.ix_(member, member)] - observed[index] = block[~np.eye(member.size, dtype=bool)].mean() - - def coefficients(size, others): - """Coefficients of the affine forward map ``r(rho)`` for one group. - - ``others`` is the grand-mean contribution of every *other* group. The - group's own contribution stays symbolic, so both sides of the ratio - remain affine in its rho and the solution is exact. Folding it into - ``others`` instead would make the numerator and denominator vanish - together whenever the group holds about half the estimates. - """ - own = size * (size - 1) / n_estimates**2 - offset = -2.0 / n_estimates + (n_estimates + others) / n_estimates**2 - numerator_slope = 1.0 - 2.0 * (size - 1) / n_estimates + own - denominator_slope = -2.0 * (size - 1) / n_estimates + own - return offset, numerator_slope, denominator_slope - - def forward(rho, offset, numerator_slope, denominator_slope): - """Residual correlation an exchangeable block with this rho would show.""" - return (rho * numerator_slope + offset) / (rho * denominator_slope + 1.0 + offset) + mask = masks.get(member.size) + if mask is None: + mask = masks[member.size] = ~np.eye(member.size, dtype=bool) + selection = np.ix_(member, member) + selections[index] = (selection, mask) + entries[index] = corr[selection][mask] + observed[index] = float(entries[index].mean()) + + # Coefficients of the affine forward map ``r(rho)``, one set per group. As + # the fixed point below iterates, the only thing that moves is the grand-mean + # term, and it arrives through ``others``, so every size-dependent part is + # settled here. + coefficients = [None] * len(members) + for index in pooled: + size = sizes[index] + own = size * (size - 1) / squared + outside = n_estimates - size + coefficients[index] = ( + 1.0 - 2.0 * (size - 1) / n_estimates + own, # numerator slope + -2.0 * (size - 1) / n_estimates + own, # denominator slope + # Numerator of the forward map's slope, before ``others`` is added + # and the divide by ``squared``. Kept in this form because the + # equivalent ``numerator_slope + offset`` is a difference of O(1) + # terms: they cancel to rounding noise exactly where the true slope + # reaches zero, so a ``<= 0`` test on that form never fires. These + # are integer differences, so they cancel exactly. + outside * (outside + 1.0), + # The smallest rho an exchangeable block of this size can hold. + -1.0 / (size - 1), + ) - def solve(r, size, others): + def solve(r, index, others): """Recover the underlying rho from an observed residual correlation. - ``r(rho)`` is a ratio of two affine functions, so its inverse has a pole - at ``numerator_slope / denominator_slope``. The pole sits far outside - [-1, 1] while a group is a small share of the estimates, but it moves - into the ordinary data range as the share grows -- at ``size / K`` of - about a third it is already near -0.9 -- and inverting across it maps - strongly anti-correlated inputs to *positive* rho. Clipping ``r`` to the - interval the forward map can actually produce keeps the inversion on the - branch the derivation covers, so the pole is unreachable by - construction and the result stays monotone in ``r``. - - The lower clip is ``rho = -1/(size - 1)``, not -1: an exchangeable block - of that size is only positive semi-definite down to there. Honouring it - also keeps the block sum -- the only thing the dependence corrections - read off a block -- non-negative, which is what they require. + The forward map has a narrow range on the negative side, and clipping + ``r`` into it is what keeps the inversion honest. At 100 estimates in 320 + an exchangeable block can only show ``r`` between -0.008 and 1, so an + observed -0.5 matches no block the derivation describes. Inverted anyway + it returns a rho far below ``floor``, and the block sum goes negative. + + The bottom of that range is the image of ``floor``, not of -1, because an + exchangeable block is only positive semi-definite down to there. + Honouring it also keeps the block sum non-negative, which is what the + dependence corrections require of a block. + + The same clip disposes of a pole. ``r(rho)`` is a ratio of two affine + functions, so its inverse blows up at ``numerator_slope / + denominator_slope``, which is near -0.92 in the same example and further + out for smaller groups. Past it the sign flips, and strongly + anti-correlated input comes back as *positive* rho. That is the rarer of + the two faults, covering about 8% of the negative half of [-1, 1] against + 91% for the floor, but the clip puts it out of reach and leaves the + result monotone in ``r``. + + The upper clip is a literal 1, not a computed bound. ``N - D`` is + identically 1, so in ``(rho * N + c) / (rho * D + 1 + c)`` the + denominator at ``rho = 1`` is ``N + c``, the same as the numerator. Every + group maps rho 1 to r 1, whatever its size. """ - offset, numerator_slope, denominator_slope = coefficients(size, others) - - # d r / d rho is (numerator_slope + offset) / (rho * denominator_slope + - # 1 + offset)**2, so this is the sign of the slope. It vanishes as a - # group approaches the whole sample, where centring has removed - # everything that identified rho and there is nothing to invert. - if numerator_slope + offset <= 0: - return np.clip(r, -1.0, 1.0) - - floor = -1.0 / (size - 1) - low = forward(floor, offset, numerator_slope, denominator_slope) - high = forward(1.0, offset, numerator_slope, denominator_slope) - r = np.clip(r, low, high) - rho = (r * (1.0 + offset) - offset) / (numerator_slope - r * denominator_slope) - return np.clip(rho, floor, 1.0) + numerator_slope, denominator_slope, outside, floor = coefficients[index] + offset = -2.0 / n_estimates + (n_estimates + others) / squared + identifying = (outside + others) / squared + + # A group spanning the whole sample has no rho to recover. R is then + # ``(1 - rho) I + rho J``, centering annihilates J, so + # ``CRC = (1 - rho) C`` and the observed correlation is ``-1/(K - 1)`` + # whatever rho was. ``identifying`` is the slope of the forward map, and + # it reaches zero here; inverting a flat map claims the strongest + # dependence there is from data that carries none. + if identifying <= 0: + return _clamp(r, -1.0, 1.0) + + shifted = 1.0 + offset + low = (floor * numerator_slope + offset) / (floor * denominator_slope + shifted) + r = _clamp(r, low, 1.0) + rho = (r * shifted - offset) / (numerator_slope - r * denominator_slope) + return _clamp(rho, floor, 1.0) # Groups interact only through the grand mean, and each group's own share is # handled exactly, so this converges immediately for a single group and in a # handful of steps otherwise. - contributions = np.zeros(len(members)) + weights = [size * (size - 1) for size in sizes] + contributions = [0.0] * len(members) for _ in range(50): - total = contributions.sum() - updated = np.array( - [ - ( - sizes[i] - * (sizes[i] - 1) - * float(solve(observed[i], sizes[i], total - contributions[i])) - if sizes[i] > 1 - else 0.0 - ) - for i in range(len(members)) - ] - ) - if np.allclose(updated, contributions, atol=1e-12, rtol=0): + total = sum(contributions) + updated = [0.0] * len(members) + for i in pooled: + updated[i] = weights[i] * solve(observed[i], i, total - contributions[i]) + # Both iterates satisfy the tolerance, so stopping on the earlier one is + # a free choice; it is the one the array form this replaced made. + if all(abs(updated[i] - contributions[i]) <= 1e-12 for i in pooled): break contributions = updated @@ -892,30 +932,34 @@ def solve(r, size, others): offset = 1.0 / (n_estimates - 1) corrected = (corr + offset) / (1.0 + offset) - # Within-group pairs are inverted through the block *mean*, which is the - # only quantity the exchangeable derivation above describes: it assumes - # every off-diagonal entry of a block equals the same rho. Sending each - # entry through the inverse separately, as though it were its own - # exchangeable block, has no such justification -- and because the inverse - # amplifies asymmetrically (a group at a third of the estimates stretches - # negatives about threefold while compressing positives), it drags the block - # mean negative even when the raw residuals are centred on zero. + # Within-group pairs are inverted through the block *mean*, the only quantity + # the derivation above describes: it assumes every off-diagonal entry of a + # block equals the same rho. Putting each entry through the inverse + # separately, as though it were its own exchangeable block, has no such + # warrant. It is also biased, because the inverse is convex: at a third of + # the estimates it stretches a residual correlation of -0.5 about tenfold but + # compresses +0.5 to 0.87 of itself. That pulls the block mean negative even + # when the raw residuals are centered on zero. # - # Heterogeneity inside a group still survives: the observed spread is - # carried over as a deviation about the recovered mean, shrunk just enough - # to stay inside [-1, 1]. Scaling a zero-mean deviation leaves the mean - # exactly at the recovered rho. - total = contributions.sum() - for index, member in enumerate(members): - if member.size < 2: - continue - block = corr[np.ix_(member, member)] - off_diagonal = ~np.eye(member.size, dtype=bool) - entries = block[off_diagonal] - mean = entries.mean() - rho = float(solve(mean, sizes[index], total - contributions[index])) - - deviation = entries - mean + # Heterogeneity inside a group still survives: the observed spread is carried + # over as a deviation about the recovered mean, and carried at full size. + # That is right to leading order. One pair moves the row sums and the grand + # total by only O(1/K), so centering shifts a block's mean by a large, + # size-dependent factor while barely touching the deviations about it. The + # measured ratio of true spread to observed spread is 1.01 at a tenth of the + # estimates, and lies between 0.8 and 1.4 by a half. ``scale`` drops below 1 + # only where a deviation would otherwise leave [-1, 1]. + # + # The deviations are zero-mean whatever the scale, so the block mean stays + # exactly at the recovered rho, and with it the block sum. + total = sum(contributions) + for index in pooled: + selection, off_diagonal = selections[index] + entry = entries[index] + mean = observed[index] + rho = solve(mean, index, total - contributions[index]) + + deviation = entry - mean highest, lowest = deviation.max(), deviation.min() scale = 1.0 if highest > 0: @@ -924,10 +968,10 @@ def solve(r, size, others): scale = min(scale, (-1.0 - rho) / lowest) scale = max(scale, 0.0) - updated_block = np.empty_like(block) + updated_block = np.empty(off_diagonal.shape) updated_block[off_diagonal] = rho + scale * deviation np.fill_diagonal(updated_block, 1.0) - corrected[np.ix_(member, member)] = updated_block + corrected[selection] = updated_block np.fill_diagonal(corrected, 1.0) return np.clip(corrected, -1.0, 1.0) diff --git a/pymare/tests/conftest.py b/pymare/tests/conftest.py index 2803055..4798281 100644 --- a/pymare/tests/conftest.py +++ b/pymare/tests/conftest.py @@ -330,6 +330,25 @@ def _shrink(corr): return _shrink +@pytest.fixture(scope="package") +def sampled_centering_shrinkage(): + """Estimate the centered correlation from data, the way a caller really would. + + The ``centering_shrinkage`` fixture applies the centering map exactly, so + every off-diagonal entry of a block comes back equal and a block-mean + inversion cannot be told apart from an entrywise one. Drawing ``n_datasets`` + samples instead gives the block the genuine spread that separates them. + """ + + def _estimate(corr, n_datasets, seed=0): + n_estimates = corr.shape[0] + factor = np.linalg.cholesky(corr + 1e-10 * np.eye(n_estimates)) + y = factor @ np.random.default_rng(seed).standard_normal((n_estimates, n_datasets)) + return np.corrcoef(y - y.mean(axis=0), rowvar=True) + + return _estimate + + @pytest.fixture(scope="package") def block_correlation(): """Build an equicorrelated-block correlation matrix and its group labels.""" diff --git a/pymare/tests/test_stats.py b/pymare/tests/test_stats.py index 729a98d..31f09ca 100644 --- a/pymare/tests/test_stats.py +++ b/pymare/tests/test_stats.py @@ -1072,12 +1072,10 @@ def test_undo_centering_shrinkage_never_reports_a_negative_block_sum( ): """The block sum is Var(sum of z) -- a variance, so it cannot be negative. - :func:`~pymare.stats.undo_centering_shrinkage` used to invert each entry of - a block separately. That inverse is a ratio of two affine functions of the - observed correlation, so it has a pole, and the pole moves into [-1, 1] once - a group is a large share of the estimates. Crossing it flipped the sign of - strongly anti-correlated pairs and dragged block sums negative, which the - combination tests can only report as an error. + The centering map is applied exactly here, so every off-diagonal entry of a + block is equal and the sum turns only on the recovered rho. That makes this a + clean check on the floor and a weak one on everything else: the entrywise + inverse this replaced passes it too. The noisy case below separates them. """ if block_size >= n_estimates: pytest.skip("a block cannot be larger than the sample it sits in") @@ -1171,14 +1169,69 @@ def test_undo_centering_shrinkage_keeps_the_block_mean_and_its_spread( assert np.all(entries >= -1.0) and np.all(entries <= 1.0) -def test_undo_centering_shrinkage_warns_when_a_group_dominates( - block_correlation, centering_shrinkage +@pytest.mark.parametrize( + ("n_estimates", "block_size", "n_datasets", "rho", "entrywise_sum"), + [(120, 40, 30, 0.0, -82.0), (320, 100, 20, 0.0, -791.0), (120, 40, 30, 0.3, 418.0)], +) +def test_undo_centering_shrinkage_survives_a_noisily_estimated_block( + n_estimates, + block_size, + n_datasets, + rho, + entrywise_sum, + block_correlation, + sampled_centering_shrinkage, ): - """d rho / d r grows like 1 / (1 - size/K)**2, so a dominant group is noise.""" - corr, groups = block_correlation(120, [(100, 0.3)]) + """A block estimated from data has spread, and the entrywise inverse ate it. + + Only sampling noise gives a block the spread that separates inverting its + mean from inverting each entry. So this is the case that pins the change + down. ``entrywise_sum`` is the median the entrywise form used to return: the + two rho-0 rows went negative in every replicate, which the dependence + corrections can only report as an error or a NaN. + + One draw is far too noisy to settle the sum, so the median carries the + accuracy claim. Non-negativity has to hold in every draw. + """ + corr, groups = block_correlation(n_estimates, [(block_size, rho)]) + truth = block_size + block_size * (block_size - 1) * rho + + sums = [] + for seed in range(32): + observed = sampled_centering_shrinkage(corr, n_datasets, seed=seed) + sums.append(undo_centering_shrinkage(observed, groups)[:block_size, :block_size].sum()) + + # Var(sum of z) cannot be negative, whatever the draw. + assert np.all(np.array(sums) >= 0.0) + assert np.median(sums) == pytest.approx(truth, rel=0.2) + # And the fix has to be the reason: the median must beat what it replaced. + assert abs(np.median(sums) - truth) < abs(entrywise_sum - truth) + + +@pytest.mark.parametrize("rho", [0.0, 0.3, 0.9]) +def test_undo_centering_shrinkage_refuses_a_group_that_is_the_whole_sample( + rho, block_correlation, centering_shrinkage +): + """One group spanning everything leaves no rho to recover, so none is claimed. + + Centering wipes out whatever identified rho here, and the first assertion + shows it: the observed correlation comes out the same for every true rho. + Anything the inversion returned past that would be invented. + + This is the case the guard in ``solve`` exists for, and it earns a test + because the guard is easy to write in a form that never fires -- one rounding + error away from zero is not ``<= 0``. In that form the identical input + returned +0.75, -0.005 and +0.286 for the three rho values below. + """ + n_estimates = 200 + corr, groups = block_correlation(n_estimates, [(n_estimates, rho)]) + observed = centering_shrinkage(corr) + off_diagonal = ~np.eye(n_estimates, dtype=bool) + assert np.allclose(observed[off_diagonal], -1.0 / (n_estimates - 1)) + + recovered = undo_centering_shrinkage(observed, groups) - with pytest.warns(UserWarning, match="of the estimates"): - undo_centering_shrinkage(centering_shrinkage(corr), groups) + assert np.allclose(recovered[off_diagonal], observed[off_diagonal]) def test_undo_centering_shrinkage_handles_several_blocks(block_correlation, centering_shrinkage):