[ENH] k h correction - #139
Merged
Merged
Conversation
A random-effects meta-regression estimates tau^2 and then treats it as known, so the model-based covariance (X'WX)^-1 understates how uncertain the coefficients are and the Wald test built on it rejects too often at small K. knapp_hartung_cov_and_dof scales that covariance by q = e'We / (K - P) and reports K - P degrees of freedom, which is the adjustment of Knapp & Hartung (2003) as metafor's rma.uni applies it. Written to reproduce metafor: the scale factor is RSS.knha / (k - p) under the inverse-variance weights, and conservative=True is metafor's test="adhoc", the modification of Rover, Knapp & Friede (2015). Nothing consumes it yet; the estimators are wired up separately. Notable details, all pinned by tests: - q is undefined when K - P < 1, since the residual sum of squares and its divisor are both zero. The inputs come back unadjusted with dof None, and a warning names the consequence. - A weighted residual sum of squares at or below machine epsilon is floored to a scale factor of exactly zero rather than to rounding noise. q is dimensionless -- scaling y by c and v by c^2 leaves every w_i e_i^2 unchanged -- so an absolute floor is meaningful here. Zero reaches get_fe_stats as a non-positive standard error, which it already reports as undefined rather than as maximally significant. - conservative=True floors q after that clamp, not before: the whole point of the modification is that the adjustment can only widen an interval, and a zeroed scale factor is the most extreme narrowing there is. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Three functions computed e'We under three spellings: np.square in correlated_effects_tau2, ** 2 in q_gen, and a named intermediate in knapp_hartung_cov_and_dof. They are the same quantity -- verified equal to 1.8e-16 across 200 random designs -- so a change to the formulation could reach one and leave the others behind. _weighted_rss is now that definition, and the three call sites read it. This also makes a claim the Knapp-Hartung docstring was only asserting true in the code: the scale factor is literally q_gen over K - P, which test_knapp_hartung_scale_factor_is_q_gen_over_its_expectation pins. q_gen cannot simply be called from there, because it refits beta and the residuals have to come from the fit that produced model_cov. Byte-identical across 144 estimator x weight-scheme x dataset-count configurations, including the rescale cases that exercise correlated_effects_tau2. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…efault Every meta-regression estimator gains small_sample_correction, and the four that estimate tau^2 default it to "knapp-hartung". This changes every standard error, p-value and interval reported without group labels; small_sample_correction="wald" returns the previous values, and is byte-identical to them across 144 configurations. WeightedLeastSquares defaults to "wald" instead. The adjustment accounts for uncertainty in an estimated tau^2 and that estimator takes tau^2 as given, which is why metafor warns when test="knha" is asked of a fixed-effects model. It is still accepted there, for a caller who plugged in a tau^2 obtained elsewhere. _robust_cov_and_dof becomes _inference_cov_and_dof and now resolves the whole choice in one place: dependent observations get the CR2 sandwich with Satterthwaite degrees of freedom, independent ones get Knapp-Hartung with K - P. They are alternatives rather than layers -- a sandwich already stops assuming the fitted weights are right, which is the entire content of the scale factor, so applying both would correct for the same thing twice. Group labels therefore take precedence and small_sample_correction is ignored, silently: the estimators do not share one default for it, so a warning keyed on "not the default" would fire on ordinary use of one estimator while staying quiet on the same request to another. Collapsing the decision into that one function let all five call sites drop both their "model_cov if robust_cov is None else robust_cov" expression and their "if dof is not None" guard, so estimators.py has the same statement count and eight fewer branches than before. The measurements behind the default are recorded separately, in validation/knapp_hartung. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Agreement with an independent implementation is the strongest evidence available for a formula, so validation/metafor regenerates 180 reference cases from a Docker image with pinned R and metafor versions, mirroring validation/robumeta. The adjustment agrees to machine precision in all 180: 1.8e-15 on the coefficients, 5.6e-16 on the standard errors, 3.0e-15 on the p-values, and 5.9e-11 absolute on the interval bounds -- the last bounded not by the adjustment but by scipy.stats.t.ppf and R's qt returning quantiles that differ in their final bits. The comparison is layered, because two questions are involved and mixing them would let one answer hide the other. Supplying PyMARE with metafor's own tau^2 isolates the adjustment from the tau^2 estimators and is what covers every metafor `method`; the closed-form estimators (FE, DL) are additionally compared end to end. Whether PyMARE's ML, REML and Hedges tau^2 match metafor's is a separate question that predates this work: ML and REML differ by up to the tau^2 search tolerance, one ML case lands on the opposite side of the tau^2 = 0 boundary, and Hedges differs by a definitional choice about whether the mean sampling variance is taken over raw or weighted rows. All three are visible with no correction applied, so none is caused by the adjustment. The README records the measurements; comparing them here would fold an optimizer's tolerance into a check on a closed-form scale factor. The pinned file keeps metafor's own test= spellings, because it records what metafor was asked. CORRECTIONS in the alignment module translates them, and a test asserts that map is total in both directions -- an unmapped metafor spelling would otherwise raise only for the cases using it, and an unmapped PyMARE value would simply never be compared. The four designs bracket the condition the literature says decides whether the adjustment behaves: sampling variances spanning a factor of 1, 250, 10,000 and 30, at 5 to 20 observations. The 10,000-fold cell is there because IntHout et al. (2014) and Rover et al. (2015) both report the adjustment exceeding its nominal level for few observations of very unequal precision. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Changing a default needs evidence, and validation/knapp_hartung is it: 144 cells at 10,000 replications each, crossing 5 to 40 observations, 1 to 3 predictors, weight ratios from 1 to 10,000, and no/moderate/large heterogeneity, for two tau^2 estimators and all three corrections. The result that decides it: "knapp-hartung" is at least as close to nominal as "wald" in every one of the 144 cells, for both estimators, to within Monte Carlo noise -- the worst excess measured is +0.0012 against a standard error of 0.0022. There is no cell where turning the correction on costs anything, including the cell where the Wald test is provably right (variances known, tau^2 truly zero), where "wald" measures 0.038 and "knapp-hartung" 0.049. Against that, the uncorrected test rejects at 0.098 on median at K = 5 and 0.310 in the worst cell. Near-zero cost when unnecessary, large gain when necessary. It does not fix the worst corner. At five observations with three predictors and a hundredfold spread of precisions -- two residual degrees of freedom, one observation carrying most of the weight -- it halves the error rather than removing it, and no reference distribution rescues that. The conservative variant helps there and overcorrects badly elsewhere (median 0.005 at K = 5, below 0.02 in a quarter of all cells), which is why it is offered rather than defaulted to. Self-contained on purpose: nothing in pymare imports from here and no measurement is pinned into the test suite. These are Monte Carlo estimates, so a pinned number would move on every run; --check enforces the three claims instead, each the reason for one decision. Not wired into CI either -- the grid takes about twenty minutes. Replications go on the parallel-dataset axis, so one fit covers a whole cell, which is what makes the grid affordable at all. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #139 +/- ##
==========================================
+ Coverage 94.76% 94.80% +0.04%
==========================================
Files 13 13
Lines 1929 1944 +15
==========================================
+ Hits 1828 1843 +15
Misses 101 101 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
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
1 participant
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.
very useful for small meta-analyses with only one image per study (no groups). The analogous correction to ch2 + satterwaite dof when there are groups.