Small-study effects and asymmetry tests

PyMetaAnalysis provides the Begg-Mazumdar rank-correlation test, classical Egger regression, and Harbord and Peters tests for two-group binary odds ratios as companions to the descriptive funnel plot. They diagnose funnel-plot asymmetry or small-study effects; none is a direct test for publication bias.

Run trim-and-fill sensitivity analysis

Use result.trim_and_fill() to estimate and mirror potentially missing studies around an iteratively estimated center. The default L0 estimator and automatic side selection follow metafor; estimator="R0" and an explicit side are also available:

import pandas as pd
import meta_analyze as ma

studies = pd.DataFrame(
    {
        "effect": [0.10, 0.15, 0.18, 0.22, 0.24, 0.29, 0.33, 0.37, 0.95, 1.15],
        "variance": [0.04, 0.035, 0.03, 0.028, 0.025, 0.022, 0.02, 0.018, 0.012, 0.01],
    }
)
studies["standard_error"] = studies["variance"] ** 0.5
result = ma.meta_analysis(studies, effect="effect", variance="variance")
filled = result.trim_and_fill(estimator="L0")
print(filled)
filled.augmented_studies
filled.adjusted_result.summary()

The augmented table marks synthetic rows with imputed=True and identifies the observed study mirrored by each fill. The adjusted fit retains the source model, interval method, and tau-squared estimator. Treat this as sensitivity analysis: funnel asymmetry can have causes other than publication bias, and trim-and-fill relies on symmetry assumptions.

The minimum input is three included studies, but this does not guarantee that every trimmed intermediate model is estimable. A random-effects refit must retain at least two studies; otherwise ConvergenceError explains the failed trimming step. Invalid estimator/iteration argument types raise InvalidStudyDataError (a MetaAnalysisError). Identical effects produce no fill and an explicit note that rank-based uncertainty is unavailable.

Automatic direction selection records the standard-error meta-regression slope in side_selection_statistic; explicit selection records no statistic. Report the original and adjusted estimates together. Do not choose among Begg, Egger, Harbord, and Peters tests after inspecting whichever produces the smallest p-value: their responses, predictors, null hypotheses, and applicability differ.

Run the Begg-Mazumdar rank test

Start from any fitted MetaAnalysisResult:

begg = result.begg_test()
print(begg.tau, begg.pvalue)

Following metafor::ranktest(), the test first centers included model-scale effects at their common-effect inverse-variance mean and standardizes each residual by sqrt(v_i - 1/sum(1/v_i)). It then computes Kendall's tau-b between those standardized residuals and sampling variances. The diagnostic is independent of the source common/random model, pooling weights, tau-squared estimator, and mean-effect inference.

By default, PyMetaAnalysis uses an exact two-sided permutation distribution when there are no tied effects or variances and fewer than 50 studies. It uses the tie-adjusted asymptotic normal approximation otherwise. Control this explicitly with:

exact = result.begg_test(exact=True)
asymptotic = result.begg_test(exact=False)
corrected = result.begg_test(
    exact=False,
    continuity_correction=True,
)

Exact inference rejects tied data rather than silently switching methods. Requesting a continuity correction with exact=None selects asymptotic inference; exact=True and continuity_correction=True are incompatible. The result records the response and predictor definitions, inference path, and the numbers of standardized-response, variance, and joint tied pairs.

This contract follows metafor::ranktest(effect, variance). R meta also describes its Begg option in terms of standardized treatment estimates, but cross-package results should still be checked against the exact versions and inputs used.

Run the classical Egger test

Start from any fitted MetaAnalysisResult:

import meta_analyze as ma

result = ma.meta_analysis(
    studies,
    effect="effect",
    standard_error="standard_error",
    model="random",
)

egger = result.egger_test()
print(egger)

Only rows with included=True enter the regression. The original pooling model, tau-squared estimate, confidence-interval method, weights, and result object are not changed. The Egger calculation is the same whether the source result was fitted as common or random effects because this first API implements only the classical regression form.

Pair the numerical result with the plot:

ax = result.funnel(contour_levels=(0.90, 0.95, 0.99))
print(egger.statistic, egger.df, egger.pvalue)

The optional contours mark two-sided significance regions around the null effect. They can make it easier to see whether an apparently missing part of the funnel lies mainly in a non-significant region, but they do not establish that studies are missing.

Harbord test for binary odds ratios

For an analysis fitted from retained two-group counts with meta_binary(..., measure="OR"), Harbord regression uses the efficient score and its variance:

trials = pd.DataFrame(
    {
        "event_treat": [1, 2, 4, 3, 5],
        "n_treat": [80, 90, 100, 110, 120],
        "event_control": [3, 5, 2, 4, 7],
        "n_control": [85, 100, 110, 100, 125],
    }
)
odds_ratios = ma.meta_binary(
    trials,
    event_treat="event_treat",
    n_treat="n_treat",
    event_control="event_control",
    n_control="n_control",
    measure="OR",
    method="MH",
)

harbord = odds_ratios.harbord_test()
print(harbord.statistic, harbord.df, harbord.pvalue)

harbord.intercept is the tested asymmetry coefficient in the standardized score regression. harbord.limit_estimate is its efficient-score limit coefficient; display_limit_estimate and display_limit_ci exponentiate that coefficient. Treat these as regression outputs, not replacements for the pooled estimate.

Harbord derives null efficient scores directly from raw treatment/control counts. It does not use continuity-corrected study log odds ratios, so changing the source study-level correction does not change the diagnostic. It is also independent of whether the source analysis used MH, inverse variance, random-effects inverse variance, or Peto pooling. Single-arm zero-event studies remain usable when their total event and non-event margins are positive.

The result records the standardized-score response, square-root score-variance predictor, equivalent weighting convention, residual dispersion, scaled-design condition number, and the fact that no continuity correction is used.

Peters test for binary odds ratios

When the analysis was fitted from retained two-group counts with meta_binary(..., measure="OR"), use the outcome-specific Peters regression:

odds_ratios = ma.meta_binary(
    trials,
    event_treat="event_treat",
    n_treat="n_treat",
    event_control="event_control",
    n_control="n_control",
    measure="OR",
    method="MH",
)

peters = odds_ratios.peters_test()
print(peters.statistic, peters.df, peters.pvalue)

The test reconstructs conventional study log odds ratios from the retained four-cell counts and the analysis's recorded study-level continuity correction. Its result is therefore independent of whether the source pooled estimate used MH, inverse variance, random-effects inverse variance, or Peto. If the source used Peto, a note makes the different study-effect construction explicit.

The tested coefficient is peters.slope, the slope of log OR on inverse total sample size. peters.limit_estimate is the extrapolated log OR as total sample size tends to infinity; display_limit_estimate and display_limit_ci are on the OR scale. These are regression extrapolations, not replacements for the pooled estimate. The result also records residual_dispersion, weight_method, the continuity-correction contract, corrected-study count, and the scaled-design condition number.

Peters regression is unavailable for generic effects whose original four-cell counts are no longer known, and for RR, RD, continuous, correlation, or diagnostic-accuracy analyses.

Inspect the Egger result

The tested coefficient is the intercept in the standardized-normal-deviate form of the Egger regression:

egger.intercept
egger.intercept_standard_error
egger.intercept_ci
egger.statistic
egger.df
egger.pvalue

The result also exposes the extrapolated effect as the standard error tends to zero:

egger.limit_estimate
egger.limit_standard_error
egger.limit_ci
egger.display_limit_estimate
egger.display_limit_ci

limit_estimate remains on the analysis model scale. Its display counterpart is exponentiated for OR/RR and back-transformed with tanh for Fisher's z correlations. It is an extrapolated regression intercept, not an automatically bias-corrected replacement for the fitted pooled estimate.

to_dict() returns a detached mapping containing the coefficients, intervals, test, method identifiers, condition number, scales, and warnings:

payload = egger.to_dict()

Pass an explicit confidence level for the two coefficient intervals, or omit it to reuse the fitted analysis level:

egger_90 = result.egger_test(confidence_level=0.90)

Changing the confidence level does not change the coefficient estimates, test statistic, or p-value.

Statistical form

The classical equation is:

y_i / s_i = alpha + beta * (1 / s_i) + error_i

where y_i is the effect and s_i is its standard error. The two-sided test is H0: alpha = 0 and uses a t distribution with k-2 degrees of freedom. The fitted beta is the limit estimate.

PyMetaAnalysis evaluates the algebraically equivalent weighted regression of y_i on s_i, using inverse sampling-variance weights and a multiplicative residual-dispersion estimate. The implementation scales the weights and design columns for numerical stability without changing the coefficients or their covariance. This contract corresponds to:

metafor::regtest(
  effect,
  variance,
  model = "lm",
  predictor = "sei"
)

The random/mixed-effects regression version of metafor::regtest() is a different model and is not silently substituted.

Applicability checks

At least three included studies are mathematically required. PyMetaAnalysis returns a result for 3 <= k < 10 but records a warning because funnel- asymmetry tests generally have low power with fewer than ten studies. The Cochrane Handbook uses ten studies as a rule of thumb and also advises against testing when study standard errors are all similar.

Exactly or numerically non-identifiable standard errors produce an error instead of an unstable coefficient. For less extreme cases, inspect the study size distribution and egger.condition_number; the library does not invent a universal cutoff for “enough” variation.

The classical Egger test is particularly problematic for some effect measures because the effect and its standard error can be inherently associated. PyMetaAnalysis therefore records an additional warning for:

  • odds ratios, for which binary-outcome alternatives such as Harbord or Peters may be preferable;
  • standardized mean differences, for which the same association can produce distorted funnel plots.

Harbord and Peters regression are available for OR analyses created by meta_binary(). A generic effect labeled GENERIC cannot reveal the original two-group counts, so neither binary-specific method accepts it.

For Harbord regression, the efficient-score variances must vary enough to identify the asymmetry intercept. Every included study must have positive total events and non-events. The method uses raw counts and never applies the source continuity correction. At least three studies are mathematically required, and the same fewer-than-ten warning applies.

For Peters regression, total sample sizes must vary enough to identify the slope. The method uses S*F/N weights, where S and F are the raw total events and non-events. Its continuity correction affects only the reconstructed study log OR, not these marginal-count weights. At least three studies are mathematically required, and the same fewer-than-ten warning applies.

Interpretation

A small p-value indicates evidence of the association defined by the selected diagnostic: standardized centered effect with sampling variance for Begg, effect with standard error for Egger, standardized efficient score with score precision for Harbord, or log OR with inverse total sample size for Peters. Possible explanations include genuine heterogeneity, design or population differences, selective outcome reporting, other non-reporting mechanisms, artefactual associations, and chance.

A large p-value does not demonstrate symmetry or exclude missing evidence, especially with few studies. A small p-value does not establish publication bias. Interpret the test alongside the funnel plot, heterogeneity, study characteristics, protocol information, and sensitivity analyses. Contour- enhanced funnels add useful significance context, but neither locate missing studies nor identify the mechanism behind asymmetry. The original Egger paper and metafor::regtest provide the methodological and software references for the Egger implementation. The Begg implementation follows the documented metafor::ranktest contract and Begg and Mazumdar (1994). The Harbord and Peters implementations follow the documented meta::metabias contract and the corresponding Harbord et al. (2006) and Peters et al. (2006) methods cited there.