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 93f4af1..16eaa9d 100644 --- a/pymare/estimators/combination.py +++ b/pymare/estimators/combination.py @@ -272,9 +272,27 @@ 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. + # + # 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( - "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..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,53 +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) - - # 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 solve(r, size, others): - """Invert the shrinkage for one group, given the other groups' share. + 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), + ) - ``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. + def solve(r, index, others): + """Recover the underlying rho from an observed residual correlation. + + 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. """ - 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 - 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) + 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 @@ -843,16 +932,46 @@ 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. - 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] - ) + # 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, 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: + 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(off_diagonal.shape) + updated_block[off_diagonal] = rho + scale * deviation + np.fill_diagonal(updated_block, 1.0) + 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 139f4dc..31f09ca 100644 --- a/pymare/tests/test_stats.py +++ b/pymare/tests/test_stats.py @@ -1064,6 +1064,176 @@ 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. + + 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") + # 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) + + +@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, +): + """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) + + assert np.allclose(recovered[off_diagonal], observed[off_diagonal]) + + 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)]