An information-geometric tour of meta-analysis

The barycenter idea

gtmeta treats a meta-analysis as a problem in the geometry of distributions. Each study enters not as a point estimate plus a standard error, but as its whole estimated Gaussian sampling distribution N(θ_i, Σ_i) — a point on the manifold of Gaussians. The pooled result is the distribution that lies closest, on average, to all of them: a weighted average taken on the curved space of distributions, then read back as an effect and its uncertainty.

N_i = N(θ_i, Σ_i) · one study N · barycenter (pooled estimate)

N = arg min_N ∑ w_i d(N, N_i)²

The geometry — the distance d — is an explicit modelling choice. gtmeta implements three: Bures–Wasserstein (optimal transport; the closed-form default), Fisher–Rao (information geometry), and Wasserstein–Fisher–Rao (unbalanced transport, which yields a robust pool). The package is deliberately conservative: wherever a classical meta-analytic limit exists, the geometric estimator is proven — and tested — to reduce to it exactly. This vignette walks that tour: pooling, heterogeneity, inference, diagnostic accuracy, and robustness.

Univariate quickstart

A study with estimate yi and standard error sei is the Gaussian N(yi, sei^2). igmi_studies() builds the study objects, igmi_precision_weights() the usual inverse-variance weights, and bw_barycenter() the Bures–Wasserstein Frechet mean by fixed-point iteration.

yi  <- c(0.12, -0.30, 0.25, 0.05, 0.40, -0.10)
sei <- c(0.10,  0.25, 0.15, 0.30, 0.20,  0.12)

studies <- igmi_studies(yi = yi, sei = sei)
w       <- igmi_precision_weights(studies)

fit <- bw_barycenter(studies, weights = w)
c(theta = fit$theta, Sigma = fit$Sigma,
  iterations = fit$iterations, converged = fit$converged)
#>      theta      Sigma iterations  converged 
#> 0.08241437 0.01969837 2.00000000 1.00000000

In this scalar case the barycenter mean is the fixed-effect estimate — nothing familiar is lost. The pooled Sigma is the (weighted) barycenter of the study variances, not the variance of the pooled mean; inference comes later. The underlying distance is available directly:

bw_dist2(studies[[1]], studies[[2]])
#> [1] 0.1989

Forest-style plot: six blue study intervals and the green pooled barycenter estimate below them.

Heterogeneity: nothing familiar is lost

The Frechet variance — the minimised objective at the barycenter — decomposes exactly into a location part (spread of the means) and a scale part (spread of the covariances):

fv <- frechet_variance(studies, w)
c(V_F = fv$V_F, V_loc = fv$V_loc, V_scale = fv$V_scale)
#>         V_F       V_loc     V_scale 
#> 0.035088450 0.032230428 0.002858022

Calibrating the location part against its known distribution under homogeneity gives the IGMI heterogeneity index I_F2 and a trace method-of-moments estimator tr_T_hat of the between-study variance. In the scalar case these are not analogues of the classical quantities; they are identical to them:

het <- igmi_heterogeneity(studies, w)

ee <- metafor::rma(yi = yi, sei = sei, method = "EE")  # Q, Higgins-Thompson
dl <- metafor::rma(yi = yi, sei = sei, method = "DL")  # DerSimonian-Laird

knitr::kable(
  data.frame(
    quantity = c("heterogeneity index", "between-study variance"),
    gtmeta = c(het$I_F2, het$tr_T_hat),
    classical = c(max(0, (ee$QE - (length(yi) - 1)) / ee$QE), dl$tau2),
    row.names = c("I_F^2  =  I^2", "tr T_hat  =  tau^2_DL")
  ),
  digits = 10
)

	quantity	gtmeta	classical
I_F^2 = I^2	heterogeneity index	0.41679358	0.41679358
tr T_hat = tau^2_DL	between-study variance	0.01794957	0.01794957

The agreement is exact (to machine precision), not asymptotic: the scalar Bures–Wasserstein geometry re-derives Higgins–Thompson I^2 and the DerSimonian–Laird estimator.

Inference

Three routes to a confidence statement, in increasing order of assumptions dropped:

igmi_location_ci() — exact chi-square region treating the study covariances as known (no heterogeneity);
igmi_pivot_ci() — exact studentized t (Hotelling for m > 1) interval via the Frechet scatter, valid at every K under proportional total covariances; no tau^2 estimate is ever plugged in;
igmi_bootstrap() — nonparametric study resampling.

loc  <- igmi_location_ci(studies, w)
piv  <- igmi_pivot_ci(studies, w)
set.seed(1)
boot <- igmi_bootstrap(studies, w, B = 500)

rbind(location = loc$ci, pivot = piv$ci, bootstrap = as.numeric(boot$theta_ci))
#>                  [,1]      [,2]
#> location  -0.03775873 0.2025875
#> pivot     -0.12397127 0.2888000
#> bootstrap -0.06247987 0.2297302

Horizontal segments comparing four 95 percent confidence intervals: location, pivot, bootstrap, and the classical DL Wald interval.

The pivot interval is the one to reach for by default: it is exact at every study count and, unlike Hartung–Knapp–Sidik–Jonkman as usually presented, needs no between-study variance estimate at all.

Diagnostic test accuracy and the SROC curve

A 2x2 diagnostic study becomes a bivariate Gaussian on the (logit sensitivity, logit specificity) plane. Because sensitivity and specificity come from disjoint patient groups, the within-study covariance is exactly diagonal — the within-study correlation that primary reports never provide is 0 by design, not a plug-in guess.

The five toy studies below show the classic threshold effect: read across them and sensitivity is traded for specificity, so the between-study correlation of the two logits is negative — which is exactly what gives a summary ROC curve its meaning.

tp <- c(96, 88, 85, 76, 70); fn <- c( 4, 12, 15, 24, 30)
fp <- c(28, 20, 14,  9,  5); tn <- c(72, 80, 86, 91, 95)

dta <- igmi_dta_studies(tp = tp, fn = fn, fp = fp, tn = tn)
dta_studies <- lapply(dta, function(s) igmi_gaussian(s$theta, s$Sigma))

dta_fit <- bw_barycenter(dta_studies)
dta_fit$theta
#> [1] 1.781012 1.880824
plogis(dta_fit$theta)          # summary sensitivity and specificity
#> [1] 0.8558218 0.8677057

The between-study covariance comes from the matrix method-of-moments estimator, and the summary ROC curve is read directly off the pooled Gaussian — by construction it coincides with the Reitsma bivariate-model SROC at the same parameters:

T_hat <- igmi_mm_cov(dta_studies)$T_hat
sr <- igmi_sroc(dta_fit, T_hat)
attr(sr, "summary_point")
#>      sens      spec 
#> 0.8558218 0.8677057

$Left: blue per-study covariance ellipses with the green pooled barycenter ellipse in logit space. Right: ROC plane with study points, the green SROC curve, and the summary operating point.$

The Wasserstein–Fisher–Rao robustness dial

Unbalanced transport may destroy mass rather than move it. For the WFR pool this becomes a single robustness tuning constant: a length scale delta such that a study farther than pi * delta from the consensus stops counting entirely. The one-atom WFR Frechet mean is exactly the Andrews-sine redescending M-estimator, and as delta -> Inf it recovers the fixed-effect mean — the dial interpolates from classical pooling to full outlier rejection. (This section uses the pure-R solver; the optional Python bridge is never needed.)

set.seed(21)
dat <- sim_contaminated(k = 12, mu = 0.3, eps = 0.25, shift = 3)
dat[dat$is_outlier, c("yi", "sei")]
#>         yi       sei
#> 2 3.378984 0.2009782
#> 3 3.313957 0.3797009
#> 6 2.817528 0.4674734

robust <- wfr_pool(dat$yi, dat$sei, delta = 0.5)
c(estimate = robust$estimate, se = robust$se,
  mass_ratio = robust$mass_ratio, n_active = robust$n_active)
#>   estimate         se mass_ratio   n_active 
#> 0.28700147 0.05369795 0.59934956 9.00000000

naive <- wfr_pool(dat$yi, dat$sei, delta = Inf)   # = fixed effect
c(robust = robust$estimate, naive = naive$estimate, truth = 0.3)
#>    robust     naive     truth 
#> 0.2870015 0.8694082 0.3000000

The retained-mass fraction mass_ratio is a built-in contamination diagnostic: mass is reported as destroyed, never silently renormalised.

Left: the WFR pooled estimate as a function of delta in violet, approaching the amber fixed-effect limit. Right: the Andrews-sine influence function with the rejection point at pi delta marked.

Classical benchmarks

The classical estimators ship alongside the geometric ones, so every comparison in this vignette can be reproduced without leaving the package:

bm <- benchmark_models(yi, sei)
knitr::kable(bm[, c("model", "estimate", "se", "ci.lb", "ci.ub",
                    "scale_par", "scale", "AIC", "BIC")],
             digits = 4, row.names = FALSE)

model	estimate	se	ci.lb	ci.ub	scale_par	scale	AIC	BIC
FE	0.0824	0.0613	-0.0378	0.2026	phi	1.0000	0.5612	0.3530
RE(ML)	0.0840	0.0736	-0.0603	0.2282	tau2	0.0071	2.3785	1.9620
UWLS	0.0824	0.0803	-0.1240	0.2888	phi	1.7147	2.1293	1.7128

And the headline reduction, stated as code:

all.equal(unname(fit$theta), bm$estimate[bm$model == "FE"])
#> [1] TRUE

Where to go next

Bures–Wasserstein (bw_barycenter()) is the closed-form default; its barycenter mean is the fixed-effect estimate.
Fisher–Rao (fr_mean(), fr_karcher_check()) is reported within its proven uniqueness region.
Wasserstein–Fisher–Rao (wfr_pool()) turns unbalanced transport into a principled robustness dial.

Multivariate pooling (igmi_studies(theta = , Sigma = )), prediction-model meta-analysis (igmi_pred_studies()), simulation generators (sim_meta(), sim_multivariate(), sim_network()), and Cochrane-style data adapters (cochrane_tidy(), igmi_study_list()) follow the same pattern: build Gaussian study objects, choose a geometry, pool.

Every estimator here is anchored to a numbered theorem and tested against it in tests/testthat/ — the test suite is the package’s second vignette.