PaperlayerView exact revision on arXiv

Source-linked AI summary

Bayesian Weighted Mendelian Randomization for Causal Inference based on Summary Statistics

Jia Zhao, Jingsi Ming, Xianghong Hu, Gang Chen, Jin Liu, Can Yang

arXiv:1811.10223v2stat.ME

TL;DR

Mendelian randomization faces challenges from weak polygenic effects and pleiotropy when inferring causal relationships from GWAS summary statistics. The paper proposes BWMR, which uses Bayesian weighting and a variational expectation-maximization algorithm, and reports statistical efficiency and computational stability in simulations and real-data analysis.

  • Problem

    Mendelian randomization must infer causal effects despite many weak SNP-exposure effects arising from polygenic architectures and other unique challenges.

  • Method

    BWMR addresses these challenges through Bayesian weighting and a variational expectation-maximization algorithm for stable and efficient causal inference.

  • Results

    BWMR was shown through simulations and real-data analysis to be statistically efficient and computationally stable.

  • Takeaways & Limitations

    BWMR provides a summary-statistics approach for causal inference that was applied to relationships involving metabolites and complex human traits.

  • Takeaways & Limitations

    Existing competing methods can be computationally inefficient or have high estimation error when many weak effects are present.

Abstract

from arXiv · show

The results from Genome-Wide Association Studies (GWAS) on thousands of phenotypes provide an unprecedented opportunity to infer the causal effect of one phenotype (exposure) on another (outcome). Mendelian randomization (MR), an instrumental variable (IV) method, has been introduced for causal inference using GWAS data. Due to the polygenic architecture of complex traits/diseases and the ubiquity of pleiotropy, however, MR has many unique challenges compared to conventional IV methods. We propose a Bayesian weighted Mendelian randomization (BWMR) for causal inference to address these challenges. In our BWMR model, the uncertainty of weak effects owing to polygenicity has been taken into account and the violation of IV assumption due to pleiotropy has been addressed through outlier detection by Bayesian weighting. To make the causal inference based on BWMR computationally stable and efficient, we developed a variational expectation-maximization (VEM) algorithm. Moreover, we have also derived an exact closed-form formula to correct the posterior covariance which is often underestimated in variational inference. Through comprehensive simulation studies, we evaluated the performance of BWMR, demonstrating the advantage of BWMR over its competitors. Then we applied BWMR to make causal inference between 130 metabolites and 93 complex human traits, uncovering novel causal relationship between exposure and outcome traits. The BWMR software is available at https://github.com/jiazhao97/BWMR.

1 Introduction

GWAS summary statistics enable Mendelian randomization without individual-level data, but polygenicity, horizontal pleiotropy, and related biases challenge causal inference. BWMR addresses weak-effect uncertainty and pleiotropy in a unified, computationally stable framework and demonstrates efficient performance and novel causal relationships.

  • Challenges: MR faces weak SNP-exposure effects, ubiquitous horizontal pleiotropy, linkage disequilibrium, sample overlap, and other summary-statistics biases.Weak-effect uncertainty must be modeled, while unaccounted horizontal pleiotropy can invalidate classical instrumental-variable assumptions and produce false positives.
  • Related work: Existing methods have complementary limitations, including ignored SNP-exposure uncertainty, high error with weak pleiotropy, biased outlier sensitivity, ad-hoc detection, inflated type I error, and numerical instability.PRESSO, Egger, GSMR, and RAPS each address subsets of these challenges but retain important statistical or computational weaknesses.
  • Contributions: BWMR jointly models uncertainty in weak GWAS effects and horizontal pleiotropy, using Bayesian weighting plus a variational expectation-maximization algorithm for stable, efficient inference.An exact closed-form correction addresses posterior covariance underestimation in variational inference.
  • Results: Simulations showed BWMR was computationally stable and statistically efficient versus related methods, while real-data analysis revealed novel causal relationships between exposure and outcome traits.The method was applied to causal inference using GWAS summary statistics.

2 Methods

BWMR estimates causal effects from GWAS summary statistics by modeling weak pleiotropic effects and identifying strong pleiotropic outliers through Bayesian weighting. A VEM algorithm enables efficient posterior approximation, complemented by a closed-form correction for underestimated posterior variance.

  • Data preparation: The model uses GWAS summary-statistic effect estimates and standard errors after selecting exposure-associated SNPs at p_Xj ≤5 × 10^-8 and applying LD clumping for independence.The resulting input is D = {γ̂, Γ̂, σ_X, σ_Y}.
  • BWMR model: BWMR models direct SNP effects as weak pleiotropic noise while using Bayesian weights to identify and downweight outliers caused by strong horizontal pleiotropy.This relaxes the exclusion restriction and supports robust causal-effect estimation when some direct effects are unusually large.
  • BWMR model: BWMR provides posterior mean and variance for the causal effect β while estimating model parameters τ^2 and σ^2 from summary statistics.It uses latent variables z = {β, π_1, w, γ} and fixed hyperparameters σ_0 = 1 × 10^6 and a_0 = 100.
  • VEM algorithm: Because exact posterior evaluation is intractable, BWMR uses mean-field variational expectation-maximization with closed-form updates for efficient and stable inference.The E-step updates the variational distribution, while the M-step optimizes model parameters to increase the ELBO.
  • Posterior variance correction: A closed-form correction inspired by linear-response methods addresses the posterior variance underestimation produced by variational inference.The correction uses sensitivity of the posterior mean under a perturbation to approximate the true posterior variance.

3 Results

Simulations showed that BWMR maintained accurate estimation, type I error control, and statistical power under horizontal pleiotropy, while real-data analyses identified established and novel metabolite–trait causal relationships. BWMR results were also stable across independent datasets.

  • Simulation study: Using independent exposure datasets avoided selection bias, whereas analyses using selected effect estimates produced biased estimates of β for all MR methods.Selection on pXj ≤ 1 × 10^-5 makes E[γ̂j|γj; pXj ≤ 1 × 10^-5] ≠ γj; independent datasets were therefore used for selection and effect estimation.
  • Simulation study: BWMR showed satisfactory estimation accuracy, type I error control, and statistical power when 20%, 50%, or 80% of IVs had horizontal pleiotropy.The authors attribute this performance to adaptive modeling of weak pleiotropic effects, Bayesian weighting against strong pleiotropy, and a stable, efficient estimation algorithm.
  • Simulation study: GSMR had the highest statistical power but failed nominal type I error control under greater horizontal pleiotropy, while Egger was most conservative and BWMR and RAPS performed similarly.GSMR’s inflated error was attributed to ignored weak pleiotropic effects and underestimated standard errors; Egger’s intercept could not model weak pleiotropy.
  • Real-data analysis: Real-data results replicated established associations, including LDL-C increasing CAD risk (β̂ = 0.13, p-value = 1.90 × 10^-25) and serum urate increasing gout risk (β̂ = 0.27, p-value = 1.49 × 10^-36).The LDL-C association was confirmed by RCTs (Ference et al., 2017), and the serum urate association supported previous individual participant data analysis.

4 Conclusion

BWMR is a GWAS summary-statistics method that models uncertainty in weak effects and weak horizontal pleiotropy while adaptively detecting outliers from a few large pleiotropic effects. Simulations and real-data analysis show that BWMR is statistically efficient and computationally stable, with an R package available online.

  • BWMR accounts for uncertainty in estimated weak effects and weak horizontal pleiotropic effects while adaptively detecting outliers caused by a few large horizontal pleiotropic effects.
  • BWMR is shown through comprehensive simulations and real-data analysis to be statistically efficient and computationally stable.
  • The BWMR R package, including the functions used in BWMR and an example application to real data, is available at https://github.com/jiazhao97/BWMR.

Supplementary Document … B.3 Inference of posterior variance of one specific latent variable

The supplementary document develops BWMR inference through VEM, derives mean-field updates, and corrects MFVB’s underestimated posterior variance using an exact closed-form covariance adjustment. It then specializes this correction to the posterior variance of the latent variable β.

  • Supplementary Document: BWMR uses VEM because exact posterior evaluation is intractable, maximizing an ELBO through E- and M-steps to approximate the posterior and increase marginal likelihood.The E-step updates the variational distribution, while the M-step optimizes model parameters.
  • A.1 E-step: Mean-field variational Bayes factorizes q(z), yielding Gaussian variational distributions for β and γj, a Beta distribution for π1, and a Bernoulli distribution for wj.These coordinate-ascent updates form the variational E-step.
  • A.2 M-step: The M-step updates σ2 directly, while τ2 requires convexity-based bounds and a bounded ELBO optimization to obtain an updating equation.Direct differentiation with respect to τ2 does not yield a closed-form update.
  • B Inference: MFVB can provide accurate posterior means but often underestimates posterior variance; BWMR therefore proposes an exact closed-form correction for accurate inference of β.The correction is inspired by linear response methods.
  • B.1 An example showing the properties of MFVB: In a multivariate normal example, MFVB gives accurate mean estimates but provides no covariance information between variables and often underestimates variance.This example motivates correcting MFVB covariance estimates.
  • B.2 An example showing the intuition for linear response variational Bayes (LRVB): LRVB uses the sensitivity of the posterior mean under perturbation at t = 0 to estimate posterior covariance, and the Gaussian example shows it can correct MFVB covariance estimation.The derivation assumes accurate MFVB posterior means for perturbations and an accurate first-order approximation at t = 0.
  • B.3 Inference of posterior variance of one specific latent variable: For BWMR, the LRVB-style inference targets only the posterior variance of β, approximating Var(z1) by perturbing the corresponding component and applying MFVB solutions.The derivation sets t = (t1, 0) and obtains the variance approximation from the perturbed variational posterior.

B.4 Detailed derivations for inference

This section derives an estimator for the posterior variance of latent β by combining cumulant-generating-function identities with MFVB properties. It also details computation of the required matrices and vectors and introduces a reparameterization to improve numerical stability.

  • Posterior variance correction: Cumulant-generating-function properties and MFVB solution properties provide the theoretical basis for correcting posterior variance estimates.The derivation perturbs the posterior and approximates it with an optimally parameterized MFVB distribution under stated conditions.
  • Posterior variance correction: The derivation combines Eqs. (38, 39, 40) to obtain an estimator of the posterior variance of β.The approach uses accurate MFVB posterior-mean estimation while avoiding repeated posterior perturbations in practice.
  • Computational details: The implementation computes the ELBO-based matrix H and vector g needed for the variance correction.The section gives the first and nonzero second derivatives of the ELBO and identifies the associated calculation steps.
  • Numerical stability: When variational parameter π_wj approaches 0 or 1, H becomes numerically noninvertible, so a reparameterization replaces the unstable boundary term with a stable matrix expression.The reparameterized formulation uses exponential-family sufficient statistics and enables stable calculation of the inverse matrix.
  • Computational details: The standard error of β is obtained from the first-row, first-column element of the derived matrix expression.The derivation separately specifies the corresponding H1 and g_m calculations.

C Comparison of different MR Methods · D More simulation results

This section compares IVW, Egger, GSMR, and RAPS as approaches to Mendelian randomization with horizontal pleiotropy, emphasizing their modeling choices and limitations. The supplied passages do not include substantive content from the “D More simulation results” subsection.

  • C Comparison of different MR Methods: IVW is presented as the standard Mendelian randomization approach and can be viewed as a meta-analysis of single causal estimates.
  • C Comparison of different MR Methods: Egger extends IVW with an intercept term to address horizontal pleiotropy, but may perform poorly with weak non-constant pleiotropy or strong pleiotropic outliers.
  • C Comparison of different MR Methods: GSMR uses HEIDI-outlier detection to remove strong horizontal-pleiotropy effects while accounting for SNP-exposure measurement error and weak linkage disequilibrium between SNPs.
  • C Comparison of different MR Methods: HEIDI-outlier tests each variant against a target variant and removes significantly different variants, but selecting an inappropriate target can bias estimates or inflate type I error.
  • C Comparison of different MR Methods: RAPS models weak horizontal pleiotropy with a random-effects variance component and can reduce outlier influence using robust losses such as Huber or Tukey’s biweight loss.
  • C Comparison of different MR Methods: Compared with RAPS, GSMR ignores weak pleiotropic effects by setting τ^2 = 0, which could lead to inflated type I errors.
  • D More simulation results: No substantive passage from the supplied input reports the “D More simulation results” subsection.

D.1 Summary-level simulations

Summary-level simulations tested BWMR and related methods across five data-generation cases, varying pleiotropy and simulation parameters. Estimation accuracy, type I error control, and statistical power were evaluated against Egger, GSMR, and RAPS.

  • D.1 Summary-level simulations: Summary-level data were simulated to test the robustness of BWMR and related methods across five cases.The simulations generated summary statistics from latent γ_j and Γ_j values, with σ_Xj and σ_Yj sampled independently from U[c, d].
  • D.1 Summary-level simulations: Case-3 modeled Γ_j = βγ_j + α_j with independent normal γ_j and pleiotropic effects, while Case-5 used Laplace-distributed α_j with rate r.Case-3 and Case-5 were among the explicitly specified data-generation settings.

D.2 Individul-level simulations · E Real data analysis

Individual-level simulations evaluated BWMR, Egger, GSMR, and RAPS under settings with and without selection bias, including varying sample sizes. Results were summarized across 100 replications.

  • D.2 Individul-level simulations: Individual-level simulations compared BWMR, Egger, GSMR, and RAPS with and without selection bias.The comparison is shown in Figure 24.
  • D.2 Individul-level simulations: The simulations varied β and used π10 = 0.08, π01 = 0.08, π00 = 0.82, SNR1 = SNR2 = 1 : 1, and a p-value threshold of 1 × 10−5.The reported β setting includes 0.02.
  • D.2 Individul-level simulations: Simulations without selection bias were denoted “.u” or “unbias,” whereas those with selection bias were denoted “.b.”These labels distinguish the two simulation conditions.
  • D.2 Individul-level simulations: The individual-level simulation results were summarized from 100 replications.This replication count applies to the reported comparison.
  • D.2 Individul-level simulations: A separate selection-bias analysis compared individual-level simulations across different sample sizes.The comparison is shown in Figure 25.
  • D.2 Individul-level simulations: The sample-size simulations used β ∈ {0.0, 0.1, 0.2, 0.3, 0.4, 0.5} and N0 = 10, 000.They also used n1 = n2, SNR2 = 1 : 1, and a p-value threshold of 1 × 10−5.

E.1 Data sources · E.2 MR implementation · E.3 Examination of selection bias

The analyses draw on GWAS sources for metabolites and complex human traits and diseases, using standardized MR preprocessing and implementations across several methods. Selection-bias analyses compare causal-effect estimates and standard errors with and without selection bias for RAPS, GSMR, and Egger.

  • E.1 Data sources: GWAS sources for metabolites are compiled in a supplementary table.
  • E.1 Data sources: GWAS sources for complex human traits and diseases are compiled in a second supplementary table.
  • E.2 MR implementation: Exposure and outcome alleles were harmonized with TwoSampleMR, and SNP effects were standardized with gsmr (Zhu et al., 2018).
  • E.2 MR implementation: Egger and RAPS used TwoSampleMR, GSMR used gsmr, and BWMR used the authors’ BWMR R package.
  • E.3 Examination of selection bias: Selection-bias comparisons for RAPS show estimated causal effects and standard errors with bias on the x-axis and without bias on the y-axis, with a dashed diagonal reference.
  • E.3 Examination of selection bias: The same selection-bias comparison design was used for GSMR and Egger, plotting estimated causal effects with standard-error bars against a dashed diagonal.

E.4 Consistency of analysis results

The section compares causal-effect analyses using Data-A versus Data-B as exposure data for RAPS, GSMR, and Egger. Each comparison displays estimated causal effects with standard-error bars against a diagonal reference.

  • RAPS: RAPS results are compared between Data-A and Data-B exposure data using estimated causal effects and standard-error bars.The comparison is shown in Figure 29, with Data-A on the x-axis and Data-B on the y-axis.
  • GSMR: GSMR results are compared between Data-A and Data-B exposure data using estimated causal effects and standard-error bars.The comparison is shown in Figure 30, with Data-A on the x-axis and Data-B on the y-axis.
  • Egger: Egger results are compared between Data-A and Data-B exposure data using estimated causal effects and standard-error bars.The comparison is shown in Figure 31, with Data-A on the x-axis and Data-B on the y-axis.

E.5 Some illustrative examples in real data analysis

Real-data examples illustrate how weak or correlated horizontal pleiotropy and selection bias can undermine several MR methods. Specifically, GSMR may produce false discoveries, Egger may yield biased estimates, and RAPS may lack robustness under correlated pleiotropy.

  • GSMR examples: Across LDL cholesterol–height and Gp–Crohn’s disease examples, GSMR can have inflated type I error under weak horizontal pleiotropy, leading to false discoveries.The LDL cholesterol–height example additionally shows that GSMR’s HEIDI outlier method may select an inappropriate pleiotropic SNP as the top variant.
  • Egger examples: In the waist circumference–L.LDL.L example, Egger is sensitive to outliers caused by a few large pleiotropic effects.The corresponding MR results and causal-inference comparisons are reported in Table 5 and Figure 34.
  • Egger examples: For age at menarche–S.VLDL.C, Egger may mistakenly detect pleiotropic IVs or outliers, producing biased estimates due to selection bias in the summary statistics.The example’s MR results and method comparisons are presented in Table 6 and Figure 35.
  • RAPS example: In the CAD–M.LDL.C example, RAPS may not be robust when horizontal pleiotropic effects correlate with SNP-exposure effects.The corresponding comparisons are shown in Figure 36 and Table 7.

E.6 MR-based causal inference between metabolites and complex traits

BWMR was applied to estimate causal effects in both directions between 130 metabolites and 92 complex human traits, with significance assessed after Bonferroni correction. Comparisons included RAPS, GSMR, and Egger, while repeated analyses exposed numerical instability in RAPS.

  • BWMR: BWMR estimated causal effects from 130 metabolites to 92 complex human traits and in the reverse direction, with significant effects identified after Bonferroni correction at level 0.05.SNP instruments were selected using the p-value threshold 5 × 10−8.
  • RAPS: RAPS produced numerically unstable estimates across repeated analyses, including cases where maximum value −minimum value exceeded 50.Each exposure–outcome pair was analyzed 300 times.
  • RAPS: RAPS estimates also differed across repetitions with 2 < maximum value −minimum value < 50 and with 0.01 < maximum value −minimum value < 2.These stability analyses were conducted for metabolite–human-trait causal-effect estimation.
  • Comparative methods: GSMR and Egger provided additional bidirectional analyses of metabolite–trait effects, using SNPs at p-value threshold 5 × 10−8 and Bonferroni correction at level 0.05.The reported analyses covered effects of metabolites on complex traits and effects of complex traits on metabolites.
Loading 1811.10223v2…