Meta-regression¶

Meta-regression relates study effect estimates to one or more study-level moderators. It is a natural extension of subgroup analysis, but it does not turn aggregate study data into individual-level evidence.

Fit a numeric moderator¶

Provide a generic effect plus exactly one of its sampling variance or standard error:

import pandas as pd
import meta_analyze as ma

studies = pd.DataFrame(
    {
        "citation": ["A", "B", "C", "D", "E", "F"],
        "effect": [0.12, 0.25, 0.41, 0.38, 0.62, 0.76],
        "se": [0.18, 0.20, 0.17, 0.23, 0.19, 0.21],
        "mean_age": [42.0, 48.0, 51.0, 55.0, 60.0, 64.0],
    }
)

result = ma.meta_regression(
    studies,
    effect="effect",
    standard_error="se",
    moderators=["mean_age"],
    study="citation",
    model="mixed",
    tau2_method="REML",
)

print(result.summary())
print(result.coefficients)
numeric_result = result

Numeric moderators are used exactly as supplied. PyMetaAnalysis does not center, scale, transform, or impute them automatically. The numerical solver uses an internal power-of-two column scaling and converts coefficients and covariance back to the supplied units; this does not change the model or the reported design matrix.

Encode categorical moderators explicitly¶

Every categorical moderator requires an ordered, complete list of levels. The first level is the treatment-coding reference:

studies = studies.assign(
    region=["Europe", "Asia", "North America", "Europe", "Asia", "North America"]
)
result = ma.meta_regression(
    studies,
    effect="effect",
    standard_error="se",
    moderators=["mean_age", "region"],
    categorical={
        "region": ["Europe", "Asia", "North America"],
    },
)

This produces terms such as region[Asia] and region[North America], each relative to Europe. The reference never depends on row order. Undeclared levels, levels absent after exclusions, and string moderators omitted from categorical= are errors rather than implicit recoding decisions. Categorical matching keeps booleans distinct from numeric levels. Python and NumPy integers are equivalent, and integer-valued floats such as 1.0 match the corresponding integer level 1. This accommodates pandas integer columns that become floating point after introducing a missing value. Non-integer floats such as 1.5 do not match integer levels.

Formula parsing, automatic interactions, splines, and polynomial terms are not implemented. Construct those columns explicitly before fitting when they are scientifically prespecified.

Array-like input¶

Use a mapping when moderators are arrays rather than DataFrame column names:

array_result = ma.meta_regression(
    effect=[0.10, 0.32, 0.45, 0.71],
    variance=[0.04, 0.05, 0.06, 0.08],
    moderators={"dose": [0.0, 1.0, 2.0, 3.0]},
    model="common",
)

A sequence such as moderators=["mean_age", "region"] is only meaningful with a DataFrame. A mapping may mix DataFrame column selectors and array-like values, provided every input has one value per row.

Models and inference¶

model="mixed" is the default. It estimates residual tau-squared after the moderators using REML, PM, or DL. model="common" fixes residual tau-squared at zero and supports normal inference only.

Setting	Coefficient tests and intervals	Joint moderator test
normal	z / normal	chi-squared
hartung_knapp	t with k-p df	F with k-p denominator df
hartung_knapp_adhoc	safeguarded t with k-p df	safeguarded F

The unmodified Hartung-Knapp result warns when its covariance is below the classic covariance. The adhoc choice explicitly applies a lower-bound scale of one. Common-effect models reject both Hartung-Knapp choices.

The design matrix must be full rank and leave positive residual degrees of freedom (k > p). PyMetaAnalysis does not silently drop collinear terms or use a pseudo-inverse.

Inspect the result¶

Meta-regression has no single pooled effect, so MetaRegressionResult does not provide a scalar estimate or ci. Inspect its coefficient table instead:

result.coefficients
result.coefficient_covariance
result.global_test
result.test_moderator("region")
result.heterogeneity
result.tau2
result.tau2_null
result.pseudo_r2
result.diagnostics

global_test tests all non-intercept terms. test_moderator(name) tests every encoded term belonging to that original moderator, so a multi-level category receives one joint test rather than separate interpretation through dummy-term p-values.

study_results retains every input row and contains fitted values, residuals, precision weights, normalized precision weights, and leverage. A regression precision weight is not a universal percentage contribution to every coefficient.

Leave-one-out sensitivity¶

Use exact deleted-model refits to inspect dependence on individual studies:

diagnostics = result.leave_one_out()

print(diagnostics.table)
print(diagnostics.coefficients)

Every successful deletion re-estimates residual tau-squared, coefficients, and inference with the same resolved model configuration. Deletions that make the design unidentifiable are retained explicitly instead of aborting the other refits. See sensitivity analysis for the result contract and interpretation limits.

Predict at moderator values¶

Prediction replays the fitted numeric and categorical encoding:

predictions = result.predict(
    pd.DataFrame(
        {
            "mean_age": [50.0, 65.0],
            "region": ["Europe", "Asia"],
        }
    )
)

Every model returns the fitted mean and its confidence interval. Mixed-effects models also return pi_low and pi_high for the distribution of true effects in a new study with those moderators. The interval does not include an additional, unknown sampling variance for a future observed estimate. Unknown categories and missing prediction inputs are rejected.

The default prediction rule follows coefficient inference: it uses a normal critical value with inference_method="normal" and t_(k-p) with either Hartung-Knapp method. Set prediction_interval_method="riley" while fitting to use a t_(k-p-1) prediction critical value instead. Riley therefore needs at least two residual degrees of freedom and is usually wider, especially in small samples. It is an explicit alternative approximation, not an automatic claim of superior coverage. The selected canonical method is retained in result.method.prediction_interval_method and reused by predict() and bubble().

Plot a single numeric moderator¶

After installing the plot extra, an intercept-containing model with exactly one numeric moderator can be displayed as a weighted bubble plot:

ax = numeric_result.bubble(
    moderator_label="Mean age (years)",
    effect_label="Treatment effect",
    show_confidence_interval=True,
    show_prediction_interval=True,
)

Bubble area is proportional to normalized fitted precision weight. The fitted line and interval bands reuse result.predict(), including the selected coefficient covariance and prediction-interval rule. A prediction band is available only for mixed-effects models.

PyMetaAnalysis rejects bubble plots for categorical, multiple-moderator, or no-intercept fits. Drawing a marginal or partial-effect line for those models requires explicit choices for the other moderator values, and the library does not silently choose them.

Residual heterogeneity and pseudo-R²¶

The result reports residual QE, I-squared, H-squared, and tau-squared. For an intercept-containing mixed model, PyMetaAnalysis also refits the same included rows without moderators and reports:

pseudo-R² = max(0, 1 - tau²_model / tau²_null)

The raw value is retained as pseudo_r2_raw. If the null-model tau-squared is zero, pseudo-R² is undefined; if moderators increase estimated tau-squared, the public value is truncated to zero and the negative raw value is retained with a warning. Pseudo-R² is not the proportion of outcome variance explained in ordinary individual-level regression.

Inspect deleted-study influence¶

Use exact deletion diagnostics when a study appears unusual or the fitted association may depend strongly on one row:

influence = result.influence()

influence.table
influence.dfbetas
influence.flagged

Every included study is omitted once, with tau-squared and coefficient inference re-estimated under the original settings. The result reports the externally standardized deleted residual, Cook's distance, and term-specific DFBETAS values. Failed reduced models remain visible rather than being silently dropped.

The result also exposes its numerical screening thresholds. These thresholds identify rows for review; they do not prove that a study is erroneous, authorize automatic exclusion, or correct for trying multiple model specifications. See sensitivity analysis for the full output and interpretation contract.

Inspect moderator collinearity¶

Use collinearity() to inspect coefficient inflation and the geometry of the fitted weighted design:

collinearity = result.collinearity()

collinearity.term_vif
collinearity.moderator_gvif
collinearity.condition_indices
collinearity.variance_proportions
collinearity.concerning_dimensions

term_vif contains one VIF and its square-root standard-error inflation factor (sif) for every encoded non-intercept term. moderator_gvif keeps all dummy terms for one categorical moderator together and reports its GVIF plus the dimension-adjusted gsif = gvif ** (1 / (2 * term_count)). This grouping makes a multi-level moderator comparable without pretending that each treatment- coded term is a separate scientific variable.

The condition table uses the fitted inverse-variance weights and normalizes each weighted design column before singular-value decomposition. It is therefore invariant to positive changes of moderator units. The result retains the original unweighted, scale-dependent design condition number separately as raw_condition_number.

The fixed references of condition index greater than 30 and variance proportion greater than 0.5 are conventional screening heuristics. A concerning dimension must meet the condition-index reference and concentrate more than half of the coefficient variance for at least two terms. The library does not attach automatic cutoffs to VIF/GVIF, remove variables, refit models, or turn these diagnostics into evidence for a preferred specification.

Test explicit linear contrasts¶

Inspect result.design_info.term_names, then define a scientific comparison by term name rather than coefficient position:

north_america_vs_asia = result.contrast(
    {
        "region[North America]": 1.0,
        "region[Asia]": -1.0,
    },
    name="North America - Asia",
)

north_america_vs_asia.table
north_america_vs_asia.contrast_matrix
north_america_vs_asia.joint_test

The default null hypothesis is C beta = 0. Use rhs= for a nonzero null:

age_threshold = result.contrast(
    {"mean_age": 1.0},
    name="Age slope equals 0.02",
    rhs=0.02,
)

For several prespecified hypotheses, pass a named mapping of mappings:

contrasts = result.contrast(
    {
        "North America - Asia": {
            "region[North America]": 1.0,
            "region[Asia]": -1.0,
        },
        "Age slope": {"mean_age": 1.0},
    },
    rhs={
        "North America - Asia": 0.0,
        "Age slope": 0.02,
    },
)

A DataFrame is also accepted: rows are named contrasts and columns are fitted term names. Unspecified fitted terms receive weight zero. Unknown terms, nonfinite weights, zero rows, duplicate names, and rank-deficient multi-row contrast matrices are errors rather than implicit repairs.

Each row reports the contrast estimate, standard error, null value, estimate-minus-null, z or t statistic, unadjusted p-value, and confidence interval for C beta. A full-row-rank set also receives a joint chi-squared or F test in joint_test. Individual p-values are deliberately not adjusted for multiple testing; prespecify the contrast family and apply an external multiplicity procedure when the scientific protocol requires one.

Missing values, provenance, and reports¶

missing="raise" identifies every missing effect, uncertainty, study label, or moderator field by row. missing="drop" excludes the entire row from every model calculation while retaining it in study_results with all applicable reasons.

result.provenance
result.method_details()
result.report().to_json()

Provenance records moderator input roles and categorical treatment coding. Reports include coefficients, their covariance, residual heterogeneity, moderator tests, encoding, diagnostics, row decisions, and the study-level interpretation warning.

See statistical methods for the equations and scope and limitations before using the model for consequential work.