An R Package for Causal Meta-Analysis by Integrating Multiple Observational Studies
Causal comparisons across multiple observational studies are obstructed by two simultaneous imbalances: covariates differ between the groups being compared, and the studies themselves differ from the larger natural population the investigator cares about. WMAP (Weighted Meta-Analysis with Pseudo-Populations) implements the unified weighting framework of Guha and Li (2024)1, which extends inverse probability weighting to settings with J studies and K groups. The package estimates multiple propensity scores for study-group membership, constructs balancing weights that transform the observed population into a covariate-balanced pseudo-population, and infers arbitrary features of counterfactual group outcomes (means, standard deviations, medians, group differences) with bootstrap variability estimates. Among the three implemented methods, FLEXOR is distinguished by explicitly maximizing effective sample size over a family of admissible pseudo-populations rather than optimizing for a single prespecified estimand.
Weighting methods for observational studies are overwhelmingly built for one study and two groups. Inverse probability weights balance covariates across a treatment contrast2 but become unstable when subjects carry extreme propensity values. Overlap weights address this by minimizing the asymptotic variance of the weighted average treatment effect3, and generalized overlap weights extend the idea to multiple groups within a single study4.
Three limitations follow. These methods are optimal only for restricted estimand types under theoretical conditions that applied data rarely satisfy, and study goals often involve estimands quite different from group mean contrasts, such as medians, percentiles, or correlations among multivariate responses. They do not account for inferential efficiency when many admissible pseudo-populations exist, including those constrained to match known population characteristics such as disease prevalences. And they do not extend to meta-analysis across multiple studies each containing more than two groups.
That last structure is common in practice. A multi-site cancer consortium contributes seven cohorts, each recording patients across two or more disease subtypes, with site-level referral patterns confounded with subtype composition. Pooling naively conflates site effects with subtype effects; analyzing sites separately discards power.
Guha and Li (2024) resolve this by generalizing the propensity score to the multiple propensity score over study-group pairs, then defining a family of pseudo-populations in which a patient's study, group, and covariates are mutually independent by design1. Within that family, existing methods reappear as special cases: integrative combined (IC) weights extend the combined pseudo-population of Li et al. (2018)3, and integrative generalized overlap (IGO) weights extend the generalized overlap pseudo-population of Li and Li (2019)4. A new member, FLEXOR, is obtained by optimizing effective sample size directly. WMAP is the first implementation of this framework5.
The package ships a demo dataset mirroring the structure of a multi-site breast cancer study from The Cancer Genome Atlas: 450 patients across 7 medical centers, split into 2 groups by subtype (infiltrating ductal versus infiltrating lobular carcinoma), with 30 demographic and clinicopathological covariates and 8 outcomes recording mRNA expression for genes implicated in breast cancer research. The distributed data are simulated to reproduce the layout and variable types of the original, which is available from the GDC Data Portal on registration; this lets users explore the methods without controlled-access credentials.
User data require five inputs: a study vector S over {1,…,J}, a group vector Z over {1,…,K}, an N × p covariate matrix X of continuous or binary measurements with factors expanded to dummies, an N × L outcome matrix Y, and, for FLEXOR only, a probability vector naturalGroupProp giving group prevalences in the target population. The framework assumes each subject belongs to one study and that all K groups are observed in every study.
The analysis is two-stage, exposed through two functions.
balancing.weights() estimates the multiple propensity score and returns normalized subject weights with the sample effective sample size. Propensity scores come from a joint (S,Z) ranger probability forest over all K × J cells. The min.node.size hyperparameter is tuned over {10, 25, 50, 100} using 5-fold cross-fitting with folds stratified by study-group cell, and the smallest candidate satisfying three weight-stability criteria is selected: fewer than 5% of subjects with own-cell propensity below 0.01, a maximum-to-median weight ratio under 10, and ESS/N above 20%. Out-of-fold predictions from the winning candidate become the final scores, truncated at a floor of 0.005 and renormalized.
Balancing weights are then formed as the ratio of pseudo-population to observed-population densities. For IC and IGO the weights follow in closed form. For FLEXOR, a two-step iterative procedure alternates between the analytically optimal tilting function at fixed study masses, given in closed form by Theorem 1 of Guha and Li (2024)1, and numerical maximization of sample ESS over those masses, run from multiple random starting points to guard against local optima.
causal.estimate() calls Stage 1, then computes weighted estimates of group means, standard deviations, medians, and mean group differences. These estimators are consistent and asymptotically normal by Theorem 2 of Guha and Li (2024)1; because N is often too small to justify the asymptotic approximation, the package reports percentile confidence intervals from B bootstrap replicates instead. Bootstrapping parallelizes through future with automatic backend selection: forking on macOS and Linux, separate sessions on Windows and RStudio. An optional outcome-modeling extension substitutes random-forest predictions for observed outcomes, with K-fold cross-fitting available to limit overfitting.
Both functions run non-blocking diagnostics: warnings fire when ESS falls below 5% of N, when the maximum-to-mean weight ratio exceeds 20, when any study-group cell holds fewer than 5 subjects, or when more than 5% of subjects have a pre-truncation own-cell propensity below the floor.
On the demo dataset, all three methods agree that counterfactual mean COL9A3 expression does not differ between subtypes, while FLEXOR additionally detects significantly greater expression variability in the ductal group, consistent with that subtype's known biological heterogeneity. CXCL12 and IGF1 both show lower counterfactual mean expression in ductal relative to lobular carcinoma, aligning with the genes' established roles in cell migration and growth signaling. In nearly every estimand examined, FLEXOR produced narrower confidence intervals than IC or IGO.
A simulation study of 250 independent synthetic datasets (N = 500 each, J = 7, K = 2) generated with known ground truth confirms the pattern. FLEXOR achieved substantially higher effective sample size than IC and IGO across all datasets, and lower absolute bias and standard deviation of the weighted average treatment effect in all 250. Notably, IGO weights are theoretically optimal for ATE estimation under homoscedastic outcomes and correct model specification4, conditions the simulation deliberately violates, and FLEXOR's advantage under that violation is the practical case for the method.
Like all weighting-based methods, the framework is sensitive to limited covariate overlap and to small study-group cell sizes; both destabilize propensity estimation and inflate weight variance, and the built-in diagnostics are designed to surface these conditions rather than correct them. The framework assumes every study observes every group, which excludes designs where a site contributes only one arm.
Confidence intervals are percentile bootstrap intervals; their finite-sample coverage has not yet been characterized in simulation, and no bias correction is applied. Bootstrap replicates reuse the node size selected on the original data and fit a single forest each, rather than re-tuning and cross-fitting per replicate, so the intervals do not reflect variability from the model-selection step. Reproducibility under parallel execution depends on the core count as well as the seed.
High-dimensional biomarker data remain a challenge. Covariates are currently passed through a truncated singular value decomposition before propensity estimation, applied to the covariates as supplied without centering or scaling, a step whose sensitivity to covariate measurement scale has not been systematically assessed.
WMAP translates a general theoretical framework for integrative causal inference into working software, making estimand-agnostic meta-analysis across multiple studies and multiple groups available to applied researchers for the first time. The design keeps the common case simple, two function calls from raw vectors to confidence intervals, while exposing the tuning parameters and diagnostics that heterogeneous real-world data demand.
Planned work includes a simulation-based characterization of bootstrap interval coverage across all three weighting methods, an automated test suite, support for hybrid designs integrating randomized trials with retrospective cohorts, and improved handling of high-dimensional biomarker covariates.