align metafor and PyMARE - #145
Merged
Merged
Conversation
PyMARE reproduces metafor's Knapp-Hartung adjustment to machine precision and has since #139, but that is one of four things rma.uni reports. The other three -- Cochran's Q and its p-value, I^2 and H, and the Q-profile interval for tau^2 -- were never compared against anything, and neither was the tau^2 point estimate outside the two closed-form estimators. Measuring them turned up one defect and located three divergences. q_profile computed both bounds by minimizing (Q(tau^2) - crit)**2 with scipy.optimize.minimize. Squaring turns a transversal crossing into a tangential minimum, so the gradient vanishes as the critical value is approached and the search stops early where Q is flattest. The upper bound was wrong by up to 3.6e-2 relative against metafor's confint.rma.uni -- 59.6127 where the root is 59.6160 on the test_stats design. Solving for the root with brentq instead brings both bounds to 1.3e-13 of metafor across the grid. The two tests that pinned the old upper bound asserted it to two decimals, which is why this never showed up as a failure. run_metafor.R now also records QE, QEp, I2, H2 and confint's tau^2 bounds for the existing 180-case grid. confint is asked for to convergence rather than at its default uniroot tolerance of eps^0.25 (~1.2e-4 relative), which would otherwise pin metafor's display precision rather than the bound it solves for. test_metafor_random_effects.py compares all of it. What agrees: Q, p(Q), logp(Q) 1.5e-13, all 60 design/model/method cells I^2, H for FE and DL 4.9e-14 tau^2 Q-profile interval 1.3e-13 Hedges tau^2, no moderators 2.2e-16 ML, REML tau^2 2.7e-5 (both profile numerically) Three divergences are pinned down by asserting their cause rather than their size, so the tests stay meaningful if a tolerance moves: - I^2 and H: PyMARE always reports the Q-based Higgins-Thompson pair. metafor reports that pair only for FE and DL, where it coincides with tau^2 / (tau^2 + v_t), and switches to the tau^2-based pair otherwise. - Hedges tau^2: metafor subtracts tr(PV) / (K - P), PyMARE subtracts sum(v) / K. Equal when the intercept is the only predictor, out by up to 0.14 relative with moderators. This refines what validation/metafor/README.md recorded: the divergence is specific to meta-regression, not general. - ML on extreme_k10 with one moderator: metafor stops at tau^2 = 0 and PyMARE at 0.0114. The README had this cell's direction backwards. The reference file is regenerated in full under the same pinned stack it already named (R 4.4.1, metafor 4.6-0) but a different BLAS, so every previously pinned number moved: at most 9.5e-14 relative on the inference path and 3.3e-9 on an ML tau^2, with dof unchanged. That is below every tolerance in the suite, but it does mean `git diff --exit-code` can no longer police this file; a numeric comparator follows. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
PyMARE's effect-size converters had no reference implementation behind them. escalc computes the same conversions, so this pins its output on an eight-design grid of summary statistics -- balanced and unbalanced groups, n from 5 to 400, zero to very large effects, correlations from 0 to 0.99 -- and compares measure by measure. Six comparisons agree and are now asserted: RM <-> MN estimate and variance, exact R <-> COR estimate, exact ZR <-> ZCOR estimate and variance, exact RMD <-> MD estimate, exact sdp <-> escalc's pooled SD, recovered as MD.yi / (SMD.yi / c(m)), exact SM, SMD estimates to the bias-correction bound below metafor has no single-group standardized mean, so SM's reference is escalc(measure="SMCC") with the second measurement set to zero and uncorrelated with the first, which reduces algebraically to m / sd with the exact correction on n - 1 degrees of freedom. PyMARE corrects for bias with 1 - 3/(4m - 1) where metafor uses the exact gamma-function c(m). Rather than pin a tolerance, the tests bound the error at 0.05 / m**2 -- the approximation's actual second order, measured at 0.043 / m**2 over the grid, worst 2.7e-3 at m = 4 and 2.0e-7 at m = 398. A first-order error would break that bound. Measuring the variances turned up three defects beyond the one PR #144 fixes, all of the same kind -- a missing pair of parentheses in expressions.json changing what the expression solves to: v_rmd solves to sd1**2/n1 - sd2**2/n2 instead of the sum, so the variance of a raw mean difference is negative whenever the second group is the more variable one, and exactly zero for two equally sized equally variable groups, which gives that study infinite weight. Four of the eight grid rows are affected. v_sm solves to A + d**2 where the noncentral-t variance is A - d**2, so the single-group Hedges' g variance is 9x to 110x too large and grows with the effect rather than being dominated by 1/n. v_d (one-sample) adds n * d**2 / j**2 where the variance subtracts d**2 / j**2, so the reported variance grows with the sample size. Each is recorded as xfail(strict=True) naming the expression and what it should be, together with the two-sample v_d that PR #144 fixes. strict is the point: correcting an expression turns the test green, pytest reports XPASS as a failure, and the marker has to go in the same change. All four were confirmed to flip to XPASS under the corresponding one-line fix, and under all four fixes together the rest of the suite -- 991 tests -- still passes, so nothing currently pins the wrong values. Two variances are divergences rather than defects and are recorded as such, with a test asserting which formula each side uses so the shape of the divergence cannot change silently: - the raw correlation variance: metafor's (1 - r**2)**2 / (n - 1) is the asymptotic sampling variance, PyMARE's (1 - r**2) / (n - 2) is the squared standard error under the null of no correlation. They are 51x apart at r = 0.99, so no tolerance relates them. - the single-group standardized-mean variances: metafor's are large-sample approximations, PyMARE's are the exact noncentral-t expressions and are the better quantity. They are 2.5x apart at n = 5, so only an order-of-magnitude bound can span them. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
PyMARE's permutation test had no reference implementation behind it. permutest does the same thing, and with exact = TRUE its p-value is a property of the data rather than of a generator, so it can be pinned: ten cases, being the designs small enough to enumerate -- 2^K sign flips for the intercept-only models up to K = 10, and K! orderings for the moderator ones at K = 5. The two disagree, for two reasons that together account for every counted permutation. The statistic. permutest counts |beta / se|; PyMARE counts |beta|. A permutation test needs a statistic whose null distribution does not move with what is being permuted away, and |beta| does move: refitting a permuted dataset re-estimates tau^2, which changes the weights and so the standard error. The two coincide only where the standard error happens to be invariant -- a fixed-effects intercept-only model under sign flipping, which is why four of the ten cases agree anyway. On unequal_k5 under DerSimonianLaird, PyMARE reports 0.5625 against permutest's 0.5. The tie. The observed estimate is computed by a different code path from the permuted ones and the two differ by one unit in the last place, so the inclusive comparison can drop the identity permutation -- the one that reproduces the observed data and must therefore count -- along with its sign-flipped mirror. That understates the p-value by 2/2^K whenever it bites: 0.033203125 against permutest's 0.03515625 on extreme_k10. metafor avoids this by comparing against |zval| - sqrt(eps). test_metafor_permutest_is_reproduced_by_the_z_statistic is the load-bearing test and it passes. It counts the same permutations of the same PyMARE fits, changing only the statistic to |beta / se| and reading the observed value out of the identity permutation in the same batch, and reproduces permutest exactly in all ten cases -- both coefficients of the moderator models included. So the permutation sets agree, the batched refits agree, and the inclusive comparison agrees; the disagreement is entirely in the statistic and the tie, which is also what the fix would be. What PyMARE currently reports is asserted by test_permutation_p_value_matches_metafor, a strict xfail naming both causes. It covers all ten cases in one test rather than parametrizing, because one of the two causes is a last-place rounding difference that need not reproduce on every platform in the test matrix while the other is structural: asserting them together keeps the xfail driven by the structural cause and out of reach of an unexpected pass on some runner. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
pymare.stats.cluster_robust_cov takes the CRn names from clubSandwich and
satterthwaite_dof implements the degrees of freedom that package pairs
with CR2, but neither had been checked against it. validation/robumeta
cannot: robumeta is the reference for the correlated-effects working model
-- how weight is spread across a study's rows -- and its model has
constant within-study weights by construction, which is precisely the
condition under which the two CR2 adjustments below coincide.
12 cases, the grid of variance column x model x closed-form estimator, on
the CSV the robumeta check already uses. That CSV carries two variance
columns for the same effects, one constant within each study and one
varying, and the pair turns out to be exactly what separates agreement
from divergence.
With the variances constant within a cluster, everything reported agrees
to machine precision over all six cases: coefficients 2.9e-15, the full
CR2 covariance including its off-diagonal entries 5.1e-15, standard errors
1.8e-15, Satterthwaite dof 2.4e-15, p-values 5.4e-15. Coefficients and
tau^2 agree in all twelve, which is what makes the divergence
attributable to the sandwich rather than to the fit.
With them varying, the covariance is out by up to 5.8e-2, the standard
errors 1.0e-2, the dof 3.9e-3 and the p-values 3.4e-2. The cause is exact.
Both implementations build the same Bell-McCaffrey matrix
B_j = W_j^-1 - X_j (X'WX)^-1 X_j'
and take its inverse square root, but in different metrics. With Psi = W^-1
the assumed target, clubSandwich forms
A_j = Psi_j^(1/2) (Psi_j^(1/2) B_j Psi_j^(1/2))^(-1/2) Psi_j^(1/2)
and _cr2_scores forms
A_j = W_j^(1/2) (W_j^(1/2) B_j W_j^(1/2))^(-1/2) W_j^(1/2)
A matrix square root does not commute with an asymmetric congruence, so
the two are equal if and only if W_j is a multiple of the identity -- when
the weights are constant within the cluster.
test_cr2_divergence_is_the_whitening_metric pins that rather than the size
of the gap: it writes both forms out from their definitions in one
function differing only in that congruence, then asserts the clubSandwich
form reproduces the pinned standard errors on every case, the PyMARE form
reproduces what PyMARE reports on every case, and the two agree when and
only when the within-cluster weights are constant. So the divergence is
not a different target covariance, a different cluster set, or a bug in
either sandwich.
PyMARE's form is Fisher & Tipton's A_j^C, which _cr2_scores' docstring
says, so it is answering a different question rather than answering this
one wrongly. But it is not clubSandwich's CR2 once the variances vary
inside a cluster, which is the common case in practice, and
cluster_robust_cov's docstring does say the CRn naming follows
clubSandwich. Which form method="CR2" should mean is a maintainer's
decision; these tests only make the choice visible and keep it from
changing by accident.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
The metafor reference values had no workflow re-verifying them against
metafor -- validation/metafor/README.md named adding one, and generalising
the robumeta comparator it would need, as a follow-up. This does that, and
brings the two new reference sets in with it.
validation/compare_reference.py replaces
validation/robumeta/compare_reference.py and serves all five reference
files. It walks any {"source": ..., "cases": [...]} document without being
told its shape, folding the keys that name a case into the label and
comparing everything else as numbers -- and raising on a value that is
neither, so a generator that stopped recording a quantity cannot pass the
check by having nothing left to compare.
The comparison has to be numeric, which is a change for metafor: its
Makefile target used `git diff --exit-code`, and that can no longer hold.
Regenerating metafor_reference.json under the same R 4.4.1 and metafor
4.6-0 the file already named, on a different BLAS, moves 1326 of its lines
-- by at most 9.5e-14 relative on the inference path and 3.3e-9 on an ML
tau^2, where an optimizer stops a step earlier or later. Hence the
per-quantity tolerance for tau^2, still three orders tighter than the
1e-4 the profiled-tau^2 test itself allows.
Verified both directions: the comparator accepts that real cross-BLAS
regeneration on all 1260 shared quantities, and rejects a single standard
error perturbed by 1e-9, a dropped field, and a non-numeric value.
.github/workflows/alignment.yml covers metafor and clubSandwich as two
matrix legs, each regenerating from its pinned image, comparing every
reference it owns, and then running its marker's tests against the values
just regenerated rather than only against the pin. robumeta keeps its own
workflow rather than becoming a third leg: its check is named "Check
robumeta alignment / PyMARE vs robumeta" and may be a required check on
the repository, which folding it in would silently rename. It now shares
the comparator, so the duplication is a workflow header and not logic.
Both R images take Rscript as their entrypoint instead of a single
generator, so one built image can run the three metafor scripts in turn.
Also here: a test_clubsandwich Makefile target and marker, a rewritten
validation/metafor/README.md recording what agrees to what tolerance and
what does not and why across all three metafor checks, and the correction
of two claims in it that measurement contradicted -- the Hedges tau^2
divergence is specific to meta-regression rather than general, and it is
metafor, not PyMARE, that stops at tau^2 = 0 on the ML boundary cell.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
The contributing guide described only the robumeta check; the metafor targets predate this branch and were never added to it, and this branch adds two more reference sets and a shared comparator. Replaces "Alignment with robumeta" with a section covering all six alignment modules, what each pins against, and where the divergences are written down -- with a pointer to read the validation README before changing an estimator, since several of the divergences are deliberate and the READMEs are the only place that says which. Also documents the strict-xfail convention the new modules use, because it has a trap in it: a contributor who fixes one of the recorded defects will see pytest report XPASS as a failure, and needs to know the marker is what to delete rather than something to work around. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #145 +/- ##
==========================================
+ Coverage 95.20% 95.26% +0.05%
==========================================
Files 13 13
Lines 2025 2050 +25
==========================================
+ Hits 1928 1953 +25
Misses 97 97 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Merging master brings in #144, which corrects the two-sample v_d expression. test_standardized_mean_difference_variance_matches_metafor was a strict xfail on exactly that, so the merge turned it green and pytest reported XPASS as a failure -- which is what strict is for. The marker comes out here, and the test becomes an ordinary passing one holding PyMARE's SMD variance to within 2 / (n1 + n2) of metafor's, the order at which the two approximations legitimately differ. Three markers remain, on v_rmd and on the two one-sample standardized-mean variances. Nothing else in the merge conflicts. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
All three are the defect PR #144 fixed for the two-sample Cohen's d, in three other expressions: a missing pair of parentheses, so the expression solves to something other than the variance it is named for. The escalc alignment check found them, and each has been recorded as a strict xfail naming the expression and its fix since that check was added. v_rmd read "v_rmd - (sd1**2 / n1) + (sd2**2 / n2)" and so solved to the *difference* of the two terms. The variance of a raw mean difference came out negative whenever the second group was the more variable one, and exactly zero for two equally sized equally variable groups -- which gives that study infinite weight in any inverse-variance meta-analysis. Non-positive on five of the eight rows of the reference grid. Now the sum, which matches escalc(measure="MD") exactly, to zero relative error on every row. v_sm read "... * j**2 * (1 / n + d**2) - d**2" -> A + d**2 v_d read "... * (1 / n + d**2) - d**2 / j**2 * n" -> A + n d**2/j**2 Both are the exact noncentral-t variances, and both had the final term added rather than subtracted; v_d additionally scaled it by n. With t noncentral-t on nu = n - 1 degrees of freedom and noncentrality lambda = delta sqrt(n), and d = t/sqrt(n): Var(d) = (n-1)/(n-3) (1/n + delta^2) - delta^2 / c^2 Var(g) = c^2 Var(d) = c^2 (n-1)/(n-3) (1/n + delta^2) - delta^2 The released expressions gave a single-group Hedges' g variance up to 93x too large and a single-group Cohen's d variance up to 6,399x too large, both growing with the effect, and v_d growing with the sample size -- backwards, and the same signature as the bug in #143. metafor has no exact counterpart to compare these two against: its single-group standardized-mean variance is the large-sample 1/n + y^2/(2n), which differs from the exact form by 2.5x at n = 5 and cannot settle an exact expression. So they are checked three ways instead -- against that approximation within a factor, against each other through the Var(g) = j^2 Var(d) identity, and for falling rather than rising with n. No existing test pinned the old values: the suite passes unchanged apart from the three markers, which come out here. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
The two misalignments with metafor that were defects rather than conventions. Both were recorded as measurements by the alignment suite before being fixed here; both change released numbers. permutation_test counted |beta| where metafor::permutest counts |beta / se|. A permutation test needs a statistic whose null distribution does not move with what is being permuted away, and |beta| does: refitting a permuted dataset re-estimates tau^2, which moves the weights and so the standard error. The two coincide only where that standard error is invariant under the permutation -- a fixed-effects intercept-only model under sign flipping -- so four of the ten pinned cases agreed anyway and the rest did not, unequal_k5 under DL coming out at 0.5625 against permutest's 0.5. Separately, the observed statistic is computed by a different code path from the permuted ones, which refit in one batched call, and the two could disagree by a unit in the last place. An exactly inclusive comparison then dropped the identity permutation -- the one that reproduces the observed data, and so must count -- and with sign flipping its mirror too, understating the p-value by 2/2**K whenever it bit: 0.033203125 against permutest's 0.03515625 on extreme_k10. The comparison now allows a relative slack of sqrt(eps), which is the same constant permutest uses for this, applied relatively rather than absolutely so it does not depend on the scale of the statistic. It sits seven orders of magnitude below the closest genuine near-tie on these designs. With both fixed, PyMARE reproduces permutest exactly in all ten pinned cases, both coefficients of the moderator models included, and test_permutation_p_value_matches_metafor stops being a strict xfail. Hedges tau^2 subtracted the mean sampling variance sum(v) / K. The term that belongs there is tr(PV) / (K - P), with P the OLS residual maker and V = diag(v); since only P's diagonal is needed that is sum_i (1 - h_i) v_i / (K - P) over the OLS leverages. The two are the same quantity when the intercept is the only predictor -- P is then I - J/K, every h_i is 1/K, and the sum collapses -- so PyMARE was applying a special case unconditionally, and tau^2 was out by up to 0.14 relative in a meta-regression. It now agrees with metafor's HE to 1.9e-15 across all twelve design-by-model cells, and intercept-only models are unchanged. test_hedges_estimator carried a comment saying metafor "always gives negligibly different values for tau2, likely due to algorithmic differences", and pinned PyMARE's 11.3881 rather than metafor's 11.3594. It was neither algorithmic nor negligible. tau^2 is now pinned at metafor's value there like every other quantity in that test. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
…rest stay v_r read (1 - r**2) / (n - 2), which is the squared standard error of r under the null hypothesis of no correlation rather than its sampling variance at the observed value. It is now (1 - r**2)**2 / (n - 1), which matches escalc(measure="COR") exactly, so R joins the measures compared on both halves rather than only the estimate. The tempting defence of the old expression was that it is conservative, being larger by (n - 1) / ((n - 2)(1 - r**2)). That does not survive contact with what a meta-analysis does with a variance: the inflation factor depends on the data, so the expression did not widen intervals uniformly, it reweighted studies against one another -- down-weighting those with strong correlations by up to 51x at r = 0.99 and pulling the pooled estimate toward the weak ones. ZR is unaffected, always agreed with metafor, and remains the measure to prefer for pooling correlations; the converter's docstring now says so. The four remaining divergences are deliberate, and the reasoning for each is now written down in the validation READMEs rather than living in a review thread. Each argument was checked rather than asserted: CR2 whitening metric. Both PyMARE's form and clubSandwich's satisfy the Bell-McCaffrey condition A_j B_j A_j' = Psi_j exactly -- 1.3e-15 and 1.7e-15 on the varying-variance grid -- and both give E[V_R] = (X'WX)^-1 under the working model, so both are exactly unbiased and PyMARE's is a CR2 in the defining sense. The condition does not determine A_j uniquely; clubSandwich takes it symmetric, PyMARE symmetric in the whitened metric (max |A - A'| of 2.8e-17 against 9.1e-2). What the whitened choice buys is that I - H_j is then identity-minus-rank-p, so its spectrum collapses to p non-unit eigenvalues at any group size and _cr2_low_rank_factors works in p x p; clubSandwich's is Psi^2 minus rank p, whose n_j distinct eigenvalues need the full decomposition -- 416x more work at n_j = 200 and 5,901x at 800. method="CR2" is kept, and the parameter and _cr2_scores now both say which CR2 it is and when it coincides with clubSandwich's. I^2 and H stay Q-based. It is the Higgins & Thompson definition the docstring cites, and estimator-independence is a feature: on unequal_k5 metafor reports I^2 of 80.31, 80.31, 98.92, 50.54 and 97.77 for FE, DL, HE, ML and REML on one dataset and one Q, where PyMARE reports 80.31 throughout. The ML/REML search tolerance stays. 2.7e-5 relative on a tau^2 whose own Q-profile interval spans a factor of 167 on that design, and bounded_scalar_min is built to fit 10^5 or more datasets in one vectorized search, so iterations multiply through the whole analysis. The bias correction stays approximate. Not only because the error is 0.043/m^2 -- a thousandth of SE(g) at m = 4 -- but because the converters are a symbolic system: substituting the exact gamma factor keeps the forward solve working and makes the reverse one raise NotImplementedError, sympy being unable to invert a ratio of gammas. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
Generalizing the HE correction term to tr(PV) / (K - P) introduced a second division by K - P, and at K == P both it and the residual sum of squares are zero in exact arithmetic. Dividing anyway let the sign of the rounding noise pick the answer: every leverage comes back as 1 +- 1e-16, so the subtracted term lands either side of zero and tau^2 came out +inf on this machine and NaN on the CI runners, where inf - inf replaced inf - (-inf). Only the NaN branch lost the Knapp-Hartung fallback warning that test_estimator_warns_and_falls_back_without_residual_dof asserts, so the suite passed locally and failed on all five CI platforms. The correction is now taken as a single division, which is also how metafor writes it -- (RSS - tr(PV)) / (K - P) -- and the saturated case is guarded explicitly. It reports no excess dispersion, which is what DerSimonianLaird and the likelihood estimators already report for these designs, and what Hedges itself already reported for K < P. Agreement with metafor's HE is unchanged at 1.9e-15 over the twelve cells, none of which is saturated. test_tau2_is_finite_without_residual_dof asserts finiteness rather than a value, across every variance estimator and both K < P and K == P. That is the property that was missing: a test pinning a number would have been just as platform-dependent as the bug. Verified to fail on the unguarded code here, where the bug presents as inf rather than as the NaN CI saw. Also installs arviz and cmdstanpy in the development environment used for this branch, so a local run collects the same tests as CI rather than skipping seventeen of them -- 1070 passed, 5 skipped, 5 deselected, which now matches the runners. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
test_tau2_is_finite_without_residual_dof, added in the previous commit to
stop rounding deciding the Hedges tau^2, immediately found the same defect
in DerSimonianLaird -- failing on the 3.10 runner while passing here,
which is precisely the platform-dependence it exists to catch. This one is
pre-existing rather than introduced by this branch.
With K == P the weighted fit is exact, so both halves of
tau^2 = max(0, (Q - (K - P)) / A)
are zero in exact arithmetic and the quotient is one rounding residue over
another. Measured here: Q = 2.5e-31 over A = -8.9e-16, giving -2.8e-16
which the floor turns into 0. On the runner the residues fall the other
way -- A underflows to exactly zero against a positive numerator -- and
tau^2 comes back +inf.
Guarded the same way as Hedges: a saturated design has no residual
dispersion to measure, so report none. That also covers K < P, where the
quotient was equally undefined and DL previously returned 2.6 on this
design.
VarianceBasedLikelihoodEstimator seeds its search scale from this function,
so ML and REML move on these designs as well -- from 0 to 1.4e-4 and from
4.0e9 to 7.4e5 at K < P. All of those are meaningless numbers for a design
with nothing left to fit; what changes is that they are now reached from a
deterministic starting point instead of a random one.
Nothing outside the degenerate case moves: DL is untouched for K > P, and
tau^2 against metafor is unchanged at 5.7e-14 (DL), 1.9e-15 (HE) and
2.7e-5 (REML) over the alignment grid, whose smallest design is K = 5 with
P at most 3.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
Codecov flagged 2 uncovered lines in the patch. Both were defensive branches, and looking at them found one of each kind worth having. pymare/results.py reshaped a 2-D permuted covariance to 3-D before taking its diagonal. inv_cov is (P, P, D) on every path permutation_test can drive -- checked against WeightedLeastSquares, DerSimonianLaird, Hedges, both likelihood estimators, the sample size-based one and the cluster-robust branch -- so the branch could not be taken. Removed, with a comment recording the invariant that makes it unnecessary. pymare/stats.py fell back to NaN when the bracket search for a Q-profile bound ran out of doublings. That one is reachable, and returns the right answer: a saturated design has K - P = 0, scipy.stats.chi2.ppf returns NaN for both critical values there, every comparison against NaN is False, so the loop runs out and both bounds come back NaN -- correct for a design with no residual left to profile. It was simply untested. test_q_profile_is_undefined_without_residual_dof pins it, and the comment above the loop now says what reaches the fallback rather than leaving it looking unreachable. That is the same saturated design that produced inf from a quotient of rounding residues in DerSimonianLaird and Hedges earlier on this branch. q_profile was already handling it correctly; this records that it does. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011o2j5ZNAdzPFL6LBqDqUsU
jdkent
commented
Sep 29, 2026
jdkent
left a comment
jdkent
left a comment
Member
Author
There was a problem hiding this comment.
The comparisons are verified, and the divergent solutions in PyMARE are justified!
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.