Statistical methods¶

This page specifies the calculations implemented by PyMetaAnalysis. It is an implementation reference, not guidance for choosing an estimand or replacing a review protocol.

Notation¶

For study i:

• y_i is the study effect on the model scale;
• v_i is its sampling variance;
• w_i is an analysis weight;
• k is the number of included studies;
• tau^2 is the estimated between-study variance.

Treatment and control are ordered exactly as their argument names imply. A positive difference measure means treatment minus control is positive.

Inverse-variance pooling¶

Common-effect weights and the pooled estimate are:

w_i = 1 / v_i
y_hat = sum(w_i y_i) / sum(w_i)
Var(y_hat) = 1 / sum(w_i)

Random-effects weights replace v_i with v_i + tau^2:

w_i* = 1 / (v_i + tau^2)
mu_hat = sum(w_i* y_i) / sum(w_i*)
Var(mu_hat) = 1 / sum(w_i*)

Reported normalized_weight values equal the selected raw weights divided by their sum. Pooled estimates, Q statistics, and tau-squared estimating equations use relative weights internally when necessary to avoid overflow without changing their ratios or equation roots. A sampling variance too small to produce a finite float64 inverse weight is rejected explicitly.

Between-study variance¶

Random-effects inverse-variance models require at least two included studies. The default estimator is REML.

DerSimonian-Laird¶

With common-effect weights, Cochran's Q, df = k - 1, and:

C = sum(w_i) - sum(w_i^2) / sum(w_i)
tau^2_DL = max(0, (Q - df) / C)

DL is closed form and records zero iterations.

Internally, C is evaluated as 2 * sum(i < j, w_i * w_j) / sum(w_i) using relative weights and cumulative sums. This avoids subtracting nearly equal totals when one study dominates. Residual calculations first subtract the highest-weight study's effect, then center those offsets. The REML score also divides by the positive weight trace before subtraction. These are algebraic rearrangements, not a rule that forces extreme inputs to zero heterogeneity. Unrepresentable relative precisions raise InvalidStudyDataError; an unrepresentable iterative bracket raises ConvergenceError.

Paule-Mandel¶

PM solves:

Q(tau^2) = sum((y_i - mu_hat(tau^2))^2 / (v_i + tau^2)) = k - 1

If the equation is already non-positive at zero, the estimate is the zero boundary. Otherwise PyMetaAnalysis brackets the non-negative root and solves it with a bounded scalar root finder.

The numerical root equation is multiplied by a positive variance scale before evaluation. This leaves its root unchanged while avoiding overflow from raw inverse weights.

Restricted maximum likelihood¶

REML solves the restricted-likelihood score equation over non-negative tau^2. As with PM, a non-positive score at zero produces a boundary estimate; otherwise a bracketed root is found iteratively.

The implemented score is evaluated in an algebraically equivalent scaled form; the positive scale factor does not change its sign or root.

result.diagnostics records convergence, iteration count, and whether the solution reached the zero boundary. Boundary means the constrained estimate is exactly zero; a positive numerical root is not relabeled as a boundary merely because it is smaller than the configured absolute root-finding tolerance. Failure to bracket or converge raises ConvergenceError; no fallback estimator is selected silently.

Confidence intervals¶

Normal interval¶

The default interval is:

estimate +/- z_(1 - alpha/2) * SE

This is available for common- and random-effects inverse-variance models and for common-effect Mantel-Haenszel and Peto pooling.

Hartung-Knapp intervals¶

For random-effects inverse-variance models:

q_star = sum(w_i* (y_i - mu_hat)^2) / (k - 1)
Var_HK(mu_hat) = q_star / sum(w_i*)
CI = mu_hat +/- t_(k - 1, 1 - alpha/2) * sqrt(Var_HK)

hartung_knapp uses this variance without modification. If it is below the classic variance, the result records a warning. hartung_knapp_adhoc instead uses:

max(Var_classic, Var_HK)

The variants are separate public choices. PyMetaAnalysis does not silently apply the safeguard to an unmodified Hartung-Knapp request.

Prediction interval¶

Random-effects inverse-variance models with at least three included studies use a k - 2 prediction rule:

mu_hat +/- t_(k - 2, 1 - alpha/2)
          * sqrt(tau^2 + Var_selected(mu_hat))

For ci_method="normal", Var_selected is the classic pooled-mean variance and the method is recorded as HTS. For either Hartung-Knapp choice, it is the same unmodified or safeguarded HK variance used by the mean interval, and the prediction method is recorded as HK-PR. This matches metafor::predict(..., predtype="Riley"). No interval is returned below three studies. An interval calculated with three or four studies carries an explicit uncertainty warning.

Heterogeneity and inconsistency¶

For inverse-variance fits, Cochran's Q uses common-effect inverse-variance weights:

Q = sum(w_i (y_i - y_hat_common)^2)
df = k - 1
p = P(chi-square_df >= Q)

For Mantel-Haenszel analyses, residuals are centered on the MH pooled estimate while the individual-effect inverse variances supply the Q weights. Peto uses the separate observed-minus-expected definition documented under Peto pooling.

Common-effect IV, MH, and Peto results use Q-based inconsistency:

I^2 = max(0, (Q - df) / Q)
H^2 = Q / df

This Q-based H^2 is untruncated and can be below one when Q < df, while I^2 is truncated at zero (also when Q is zero). Consequently, H^2 = 1 / (1 - I^2) does not hold after that truncation. This existing convention differs from the tau-squared-based definition below.

Random-effects results use the fitted tau^2 and a typical within-study variance:

C = sum(w_i) - sum(w_i^2) / sum(w_i)
v_typical = (k - 1) / C
I^2 = tau^2 / (tau^2 + v_typical)
H^2 = 1 + tau^2 / v_typical

Q-profile confidence interval¶

For a candidate non-negative value t, define generalized Q using weights w_i(t) = 1 / (v_i + t):

mu_hat(t) = sum(w_i(t) y_i) / sum(w_i(t))
Q(t) = sum(w_i(t) (y_i - mu_hat(t))^2)

Under the conventional normal random-effects model with known sampling variances, Q(tau^2) follows a chi-squared distribution with k-1 degrees of freedom. For confidence level 1-alpha, PyMetaAnalysis numerically solves:

Q(tau2_low)  = chi2_(k-1, 1-alpha/2)
Q(tau2_high) = chi2_(k-1, alpha/2)

and constrains both bounds to the non-negative parameter space. If the lower root is negative but the upper root is non-negative, the reported lower bound is zero. If both roots are negative, the formal confidence set is empty; Tau2ConfidenceInterval records is_empty=True and uses [0, 0] as the explicit constrained display, following metafor's output convention.

The tau, I-squared, and H-squared bounds are monotonic transformations of the same tau-squared bounds. The interval is exact only under its model assumptions; in applications with estimated sampling variances or non-normal effects it should be interpreted as an approximation. See Viechtbauer (2007).

result.i2_method records q_based or tau2_typical_variance. I-squared is a proportion internally and is formatted as a percentage in human-readable output. With one study, Q is zero and its p-value, I-squared, and H-squared are unavailable.

Begg-Mazumdar rank-correlation test¶

For included model-scale effects y_i and sampling variances v_i, first compute the common-effect center and standardized residual used by metafor::ranktest():

w_i = 1 / v_i
theta = sum(w_i * y_i) / sum(w_i)
v_bar = 1 / sum(w_i)
y_i* = (y_i - theta) / sqrt(v_i - v_bar)

Then define Kendall's concordance statistic between y_i* and v_i:

S = sum_(i<j) sign((y_i*-y_j*) * (v_i-v_j))
tau_b = S / sqrt((n0-n1) * (n0-n2))
n0 = k * (k-1) / 2

Here n1 and n2 are the numbers of tied pairs in standardized responses and variances. Without ties and for k < 50, the default two-sided p-value uses the exact permutation distribution of S. Otherwise the test uses the standard tie-adjusted null variance of S and a normal approximation. An optional continuity correction moves nonzero S one unit toward zero before standardization.

The implementation computes tau-b with an O(k log k) rank algorithm and calculates tie moments without pairwise matrices. The result is independent of the fitted pooling model and matches metafor::ranktest(). It requires at least three studies and variation in both effects and variances. See small-study effects for inference controls and interpretation limits.

Classical Egger regression test¶

For included study effect y_i, sampling variance v_i, and standard error s_i = sqrt(v_i), the classical Egger regression is:

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

This is fitted as the algebraically equivalent weighted regression:

y_i = beta + alpha * s_i + epsilon_i
weight_i = 1 / v_i

Let X contain an intercept and s_i, let W = diag(1/v_i), and let df = k-2:

b = (X' W X)^(-1) X' W y
phi = sum(w_i * (y_i - X_i b)^2) / df
Cov(b) = phi * (X' W X)^(-1)

The implementation uses a scaled singular-value decomposition rather than an explicit inverse. The second coefficient is alpha, reported as the EggerTestResult.intercept, and is tested against zero with a two-sided t test using k-2 degrees of freedom. The first coefficient is beta, reported as the model-scale limit_estimate; its interval uses the same t distribution.

This is the weighted-regression-with-multiplicative-dispersion version produced by metafor::regtest(model="lm", predictor="sei"). It is distinct from an additive random/mixed-effects meta-regression. The test requires at least three studies and identifiable variation in standard errors. Interpretation and effect-measure limitations are documented under small-study effects.

Harbord efficient-score test for binary odds ratios¶

For an included two-group binary study, let a_i be treatment events, n_ti and n_ci the group sizes, N_i=n_ti+n_ci, S_i the total events, and F_i=N_i-S_i. The null efficient score and score variance are:

Z_i = a_i - S_i * n_ti / N_i
V_i = n_ti * n_ci * S_i * F_i / (N_i^2 * (N_i - 1))

Harbord regression is fitted in the standardized form:

Z_i / sqrt(V_i) = alpha + beta * sqrt(V_i) + epsilon_i
phi = sum(epsilon_i^2) / (k - 2)
b = [alpha, beta]'
Cov(b) = phi * (X' X)^(-1)

This is algebraically equivalent to fitting Z_i/V_i on 1/sqrt(V_i) with weights V_i. alpha, reported as HarbordTestResult.intercept, is tested against zero with a two-sided t_(k-2) test. beta is reported as the efficient-score limit_estimate and is also available with an exponential display transformation; it is not an automatic replacement for the source pooled estimate.

The implementation calculates V_i through bounded group fractions and raw event margins, then uses a column-scaled SVD rather than an explicit inverse. It does not apply or depend on study-level continuity corrections, source study-effect estimators, or source pooling weights. This matches meta::metabias(method.bias="Harbord") for two-group odds ratios. Applicability and interpretation limits are documented under small-study effects.

Peters regression test for binary odds ratios¶

For each included two-group binary study, let y_i be its conventional log odds ratio, N_i its raw total sample size, S_i its raw total event count, and F_i=N_i-S_i its raw total non-event count. Peters regression fits:

y_i = beta_0 + beta_1 * (1 / N_i) + epsilon_i
w_i = S_i * F_i / N_i
df = k - 2

With X_i = [1, 1/N_i] and W=diag(w_i), the coefficients and multiplicative dispersion covariance are:

b = (X' W X)^(-1) X' W y
phi = sum(w_i * (y_i - X_i b)^2) / df
Cov(b) = phi * (X' W X)^(-1)

The implementation uses scaled square-root weights and a column-scaled SVD; it does not explicitly invert the information matrix. The tested asymmetry coefficient is beta_1, reported as PetersTestResult.slope and compared with zero using a two-sided t_(k-2) test. beta_0 is the log-OR limit estimate as total sample size tends to infinity and is also exposed on the OR display scale.

Study log odds ratios are reconstructed from retained raw four-cell counts using the source result's recorded study-level continuity correction. The S*F/N weights and 1/N predictor always use uncorrected raw margins. This matches meta::metabias(method.bias="Peters") for two-group odds ratios and does not depend on the source pooling model. Applicability and interpretation limits are documented under small-study effects.

Duval-Tweedie trim-and-fill¶

For a chosen missing side, effects are oriented so the algorithm trims the upper tail, sorted, and centered on the current inverse-variance pooled fit. Absolute centered effects receive first-occurrence ranks with their original signs. L0 uses the positive signed-rank sum; R0 uses the most extreme negative signed rank. The non-negative integer missing count is iterated until stable, then the selected extreme effects are mirrored around the final center while retaining their sampling variances. Observed and imputed effects are refitted with the source common/random model and tau-squared estimator.

With side=None, the sign of the sampling-standard-error coefficient from a source-model-compatible meta-regression selects the side, matching metafor::trimfill(). R0 reports its conventional test of k0=0; L0 does not manufacture a p-value. This is an exploratory sensitivity method whose symmetry assumption can fail under genuine heterogeneity or design differences.

First-occurrence ranks are assigned after a stable sort of the oriented effects. Tied inputs have dedicated R comparisons, but those cases do not establish equivalence for every possible tie pattern. If all included effects are identical, no studies are imputed and missing-count uncertainty is unavailable; this is an explicit boundary policy, whereas the pinned metafor R0 implementation errors on that input. A negative L0 variance approximation likewise yields an unavailable standard error, not zero. Random-effects intermediate refits must retain at least two studies; otherwise the operation raises ConvergenceError without an adjusted fit.

Meta-regression¶

For k study effects, let X be the full-rank k-by-p design matrix. It contains an intercept by default plus numeric moderator columns and treatment- coded categorical terms:

y = X beta + u + epsilon
epsilon ~ N(0, V),  V = diag(v_i)
u ~ N(0, tau^2 I)

Given residual tau-squared:

W = diag(1 / (v_i + tau^2))
B = (X' W X)^(-1)
beta_hat = B X' W y
P = W - W X B X' W
residual df = k - p

The implementation power-of-two scales design columns, uses a weighted QR solve, and returns coefficients and covariance in the supplied moderator units. This makes rank checks and fitted results invariant to representable positive changes of units without silently changing the user-visible design. For a design containing a constant column, residual calculations first anchor effects on a canonical observed value so a large common location does not erase residual geometry. A rank-deficient design or k <= p is rejected; terms are not silently removed and a pseudo-inverse is not used.

Residual tau-squared¶

For generalized DL, define W0 and P0 at tau-squared zero:

QE0 = y' P0 y
C = trace(P0)
tau^2_DL = max(0, (QE0 - (k-p)) / C)

PM solves:

y' P(tau^2) y = k - p

REML solves the restricted score:

y' P(tau^2)^2 y - trace(P(tau^2)) = 0

All estimators are constrained to non-negative tau-squared. PM and REML use a bracketed scalar root; a non-positive equation at zero returns the boundary. Residual tau-squared requires k-p > 0. A positive root remains an interior solution regardless of its size relative to the numerical tolerance. Failure raises ConvergenceError without estimator fallback.

Generalized DL and REML need trace(P(t)), which can be much smaller than the two matrix traces in its usual subtractive expression. Internally this is evaluated as the positive sum sum(w_i * (1-h_i)); high-leverage complements use determinant ratios so precision imbalances do not cancel the residual information. The REML equation is divided by this positive trace before its final subtraction. These are algebraic rescalings and do not change an estimating-equation root.

Coefficient and moderator inference¶

Normal inference uses covariance B, z coefficient tests, and a chi-squared Wald test for one or more moderator terms.

For both Hartung-Knapp choices:

q = y' P(tau^2) y / (k-p)
Cov_HK(beta_hat) = q B

Coefficient tests use a t distribution with k-p degrees of freedom. Joint tests divide their Wald statistic by the number of tested terms and use an F distribution with that numerator df and k-p denominator df. hartung_knapp_adhoc replaces q with max(1, q); unmodified hartung_knapp retains q and warns when it is below one.

With exactly zero residual dispersion, unmodified HK covariance is zero and coefficient intervals collapse to points. Scalar tests report an undefined statistic for 0 / 0 and signed infinity for a nonzero coefficient divided by zero. Global, moderator, and joint contrast tests are unavailable because their covariance cannot be inverted. The result records a warning and reports serialize unavailable statistics as JSON null. No automatic safeguard is applied; hartung_knapp_adhoc remains a separate explicit choice.

The global moderator test covers every non-intercept term. test_moderator(name) covers every encoded term for one original moderator, so a multi-level category receives a joint test.

Residual heterogeneity and pseudo-R-squared¶

Residual heterogeneity uses:

QE = y' P0 y
QE ~ chi-square_(k-p) under the null

Common models derive residual I-squared and H-squared from QE. Mixed models use:

v_typical = (k-p) / trace(P0)
I^2_residual = tau^2 / (tau^2 + v_typical)
H^2_residual = 1 + tau^2 / v_typical

For mixed models with an intercept, the same rows are refitted with an intercept-only design and the same tau-squared estimator:

R^2_raw = 1 - tau^2_model / tau^2_null
R^2 = max(0, R^2_raw)

The raw value is retained. Pseudo-R-squared is unavailable when the null tau-squared is zero and is not equivalent to ordinary regression R-squared.

Meta-regression prediction¶

For a new design vector x0:

fitted = x0' beta_hat
Var_mean = x0' Cov(beta_hat) x0

Mean-effect intervals use the selected normal or t_(k-p) critical value. A mixed model additionally reports a true-effect prediction interval:

fitted +/- critical * sqrt(tau^2 + Var_mean)

The default critical value is standard normal for normal inference and t_(k-p) for either Hartung-Knapp mode. The optional Riley rule uses t_(k-p-1) and therefore requires k-p >= 2. Both are approximations because tau-squared is estimated; the Riley option is not selected automatically. Neither interval includes a sampling variance for a future observed effect. See ADR 0003 for the compatibility and validation decision.

Meta-regression deletion diagnostics¶

For every included study i, influence() deletes that study and refits the model with the same configuration, including a fresh residual tau-squared estimate. Let beta_(-i), tau^2_(-i), and q_(-i) denote the deleted-model coefficient estimate, residual variance estimate, and inference scale. The deleted prediction and externally standardized residual are:

pred_(-i) = x_i' beta_(-i)
e_(-i) = y_i - pred_(-i)
SE(e_(-i)) = sqrt(
    q_(-i) * (v_i + tau^2_(-i))
    + x_i' Cov(beta_(-i)) x_i
)
rstudent_i = e_(-i) / SE(e_(-i))

The inference scale is one for normal inference, the deleted-model Hartung-Knapp scale for unmodified Hartung-Knapp, and its lower-bounded value for the safeguarded variant.

Let delta_i = beta_hat - beta_(-i). Cook's distance is:

Cook_i = delta_i' Cov(beta_hat)^(-1) delta_i

For DFBETAS, the deleted tau-squared is inserted into the full design matrix:

B_i = (X' W(tau^2_(-i)) X)^(-1)
DFBETAS_ij = delta_ij / sqrt(q_(-i) * B_i[j,j])

This matches the deletion definitions used by R metafor for the supported independent-effect models. The implementation uses linear solves and does not substitute a pseudo-inverse.

The output provides three transparent screening references:

potential outlier: abs(rstudent_i) > z(0.975)
Cook's-distance flag: Cook_i > median(chi-square_p)
DFBETAS flag: any abs(DFBETAS_ij) > 1

They are heuristic review aids. The residual reference is pointwise rather than multiplicity-adjusted, and no flag automatically removes a study.

Meta-regression collinearity diagnostics¶

VIF/GVIF uses the correlation matrix R of the classic coefficient covariance for non-intercept terms. A scalar Hartung-Knapp adjustment does not change this correlation matrix. For a term set S with m encoded coefficients:

GVIF(S) = det(R_SS) det(R_notS,notS) / det(R)
GSIF(S) = GVIF(S)^(1 / (2m))

For one coefficient, GVIF is the ordinary VIF and GSIF is its square-root standard-error inflation factor. Categorical terms are grouped by their original moderator. The library reports values without applying an automatic 5-or-10 cutoff.

Condition diagnostics use the actual fitted precision weights. Let D contain the Euclidean column norms of W^(1/2) X, and define:

Z = W^(1/2) X D^(-1)
Z = U diag(s_1, ..., s_p) V'
condition_index_d = s_1 / s_d

For coefficient j and singular dimension d, the variance-decomposition proportion is:

component_jd = V[j,d]^2 / s_d^2
proportion_jd = component_jd / sum_l(component_jl)

Column normalization makes these condition indices invariant to positive changes of moderator units. The original unweighted design condition number is also returned separately and remains scale-dependent.

A dimension is marked concerning only when its condition index exceeds 30 and more than 50% of the variance of at least two coefficients is concentrated in that dimension. These conventional references are descriptive screening heuristics; no term is removed or model refitted automatically.

Meta-regression linear contrasts¶

Let C be a user-supplied, full-row-rank matrix with one column per fitted coefficient, and let d be the supplied rhs vector. With fitted coefficient vector beta_hat and its selected normal or Hartung-Knapp covariance Sigma:

theta_hat = C beta_hat
Cov(theta_hat) = C Sigma C'
SE(theta_hat_j) = sqrt(Cov(theta_hat)[j,j])
test_j = (theta_hat_j - d_j) / SE(theta_hat_j)

Individual tests use the same distribution as coefficient inference: standard normal for inference_method="normal" and t_(k-p) for either Hartung-Knapp choice. Confidence intervals are centered on theta_hat, not on theta_hat - d.

For q contrast rows, the joint Wald quantity is:

W = (C beta_hat - d)' (C Sigma C')^(-1) (C beta_hat - d)

Normal inference compares W with chi-squared on q degrees of freedom. Hartung-Knapp inference reports W/q against F_(q,k-p). The implementation uses a linear solve and rejects a rank-deficient contrast matrix instead of using a pseudo-inverse.

Term-name mappings and labeled DataFrames are expanded to the complete fitted term order with zero weight for unspecified terms. Individual p-values are reported without multiplicity adjustment. The joint test evaluates the whole prespecified set but is not a substitute for a protocol-appropriate multiplicity strategy.

Binary study effects¶

For a treatment/control 2-by-2 table:

	Event	No event
Treatment	a	b
Control	c	d

Odds ratio¶

y_i = log(a d / (b c))
v_i = 1/a + 1/b + 1/c + 1/d

Risk ratio¶

y_i = log((a / (a + b)) / (c / (c + d)))
v_i = 1/a - 1/(a + b) + 1/c - 1/(c + d)

OR and RR results remain on the log model scale. Display properties exponentiate them. Continuity-correction and double-zero/double-all behavior is specified in zero-event studies.

Risk difference¶

With p_t = a / (a + b) and p_c = c / (c + d):

y_i = p_t - p_c
v_i = p_t (1 - p_t) / (a + b) + p_c (1 - p_c) / (c + d)

The effect always uses raw counts. The variance uses corrected counts for every study selected by correction_scope; with the default only_zero_studies scope, this includes every table containing a zero cell. The separate rd_zero_variance policy determines whether tables whose raw variance is exactly zero are corrected or excluded.

Mantel-Haenszel pooling¶

MH is implemented only for common-effect OR, RR, and RD. It uses raw tables by default and has a continuity-correction option separate from the study-effect correction. Every included stratum must have a positive total count; empty inputs and zero-total strata are rejected before ratio or variance arithmetic.

For OR, with n_i = a_i + b_i + c_i + d_i:

OR_MH = sum(a_i d_i / n_i) / sum(b_i c_i / n_i)

For RR, with treatment total n1_i = a_i + b_i and control total n0_i = c_i + d_i:

RR_MH = sum(a_i n0_i / n_i) / sum(c_i n1_i / n_i)

The OR and RR model-scale estimates are the logarithms of the pooled ratios. Their standard errors use the Greenland-Robins equations implemented in the estimator.

For RD, define w_i = n1_i n0_i / n_i, W = sum(w_i), and RD_i = a_i / n1_i - c_i / n0_i. The MH point estimate is:

RD_MH = sum(w_i RD_i) / W

Its Sato-Greenland-Robins sampling variance is calculated from:

A = sum(c_i (n1_i / n_i)^2
        - a_i (n0_i / n_i)^2
        + (n1_i / n_i) (n0_i / n_i) (n0_i - n1_i) / 2)

B = sum(a_i (n0_i - c_i) / n_i
        + c_i (n1_i - a_i) / n_i)

Var(RD_MH) = (RD_MH A + B / 2) / W^2

RD remains on the identity scale and its normal confidence interval is not clipped to [-1, 1]. The result records mh_rd_variance="Sato-Greenland-Robins" in its method options.

An undefined raw OR/RR estimator or a non-positive RD sampling variance raises InvalidStudyDataError. If the review protocol permits it, an explicit positive mh_continuity_correction can make a sparse fit estimable; the decision is recorded separately from the correction used for study effects.

Peto one-step odds ratio¶

Peto is implemented only for common-effect OR with a normal confidence interval. For study i, let O_i = a_i, m_i = a_i + c_i, treatment and control totals be n1_i and n0_i, and N_i = n1_i + n0_i. The expected treatment events and their hypergeometric information are:

E_i = m_i n1_i / N_i
V_i = m_i (N_i - m_i) n1_i n0_i / (N_i^2 (N_i - 1))

The one-step study estimates and their variances are:

y_i = (O_i - E_i) / V_i
v_i = 1 / V_i

The pooled log odds ratio, its variance, and heterogeneity statistic are:

y_Peto = sum(O_i - E_i) / sum(V_i)
Var(y_Peto) = 1 / sum(V_i)
Q_Peto = sum(((O_i - E_i) - y_Peto V_i)^2 / V_i)

Normalized study weights are V_i / sum(V_i). Pooling and Q always use raw tables; continuity_correction affects only the displayed one-step study estimates and variances. Double-zero and double-all rows are excluded from all synthesis calculations.

The approximation is intended for rare outcomes, similar treatment/control group sizes within studies, and effects that are not large. Every result records this caution. Peto is not an automatic replacement for MH or inverse-variance pooling when those conditions are doubtful. See ADR 0006.

Continuous study effects¶

Mean difference¶

y_i = mean_treat - mean_control
v_i = sd_treat^2 / n_treat + sd_control^2 / n_control

The variance is unpooled and does not require equal group variances.

Standardized mean difference¶

The pooled within-study variance and Cohen's d are:

df = n_treat + n_control - 2
s_pooled^2 = ((n_treat - 1) sd_treat^2
              + (n_control - 1) sd_control^2) / df
d = (mean_treat - mean_control) / s_pooled

PyMetaAnalysis applies the exact gamma-function correction:

J(df) = Gamma(df / 2) / (sqrt(df / 2) Gamma((df - 1) / 2))
g = J(df) d

The implemented LS sampling variance is:

v_i = 1/n_treat + 1/n_control + g^2 / (2 (n_treat + n_control))

Correlation study effects¶

meta_correlation(measure="ZCOR") transforms each independent Pearson correlation with Fisher's r-to-z transformation:

y_i = atanh(r_i) = 0.5 log((1 + r_i) / (1 - r_i))
v_i = 1 / (n_i - 3)

The common- or random-effects model, heterogeneity statistics, confidence interval, and prediction interval are calculated on the Fisher's z scale. Display values use the inverse transformation:

r = tanh(z)

The accepted domain is -1 < r_i < 1 and integer n_i >= 4. The strict correlation boundary prevents infinite transformed effects; the sample-size boundary ensures a finite, strictly positive sampling variance. The pooled display correlation is the inverse transform of the pooled z value, not a direct arithmetic average of raw correlations.

This implementation assumes one independent correlation per study. It does not construct covariance matrices for multiple dependent correlations from the same sample. See ADR 0007.

Subgroup differences¶

Subgroup fits use the same selected model independently within each group. For random effects, tau-squared is estimated separately for every subgroup and the overall analysis whenever a subgroup contains at least two included studies. A single-study subgroup cannot identify tau-squared or Hartung-Knapp inference, so it is retained as that study's common-effect estimate and normal interval with an explicit warning. The overall fit remains random effects.

The formal subgroup-differences test weights subgroup pooled estimates by the inverse square of their classic model standard errors:

W_g = 1 / SE_g^2
theta_bar = sum(W_g theta_g) / sum(W_g)
Q_between = sum(W_g (theta_g - theta_bar)^2)
df_between = number_of_groups - 1

Q_between is compared with a chi-squared distribution. This is not a test obtained by comparing separate subgroup p-values. Hartung-Knapp adjustments affect subgroup confidence intervals but are not substituted for the classic model variances in this Wald test.

Numerical and software references¶

The implementation is independent Python code. Regression fixtures generated by R metafor, hand calculations, invariants, and numerical edge cases are described in validation. Exact supported combinations and defaults are listed in the API reference.