Comparison of statistical properties of MR methods

We assessed the statistical performances of several methods designed to either detect (MR-PRESSO global test, the Q test^17,18, Q (modified) test¹⁹, Q’ test ^17,18 and Q’ (modified) test¹⁹) or correct for horizontal pleiotropy (correction of the global average horizontal pleiotropic effect using MR-Egger regression^20,21 or Multi-variable MR^22,23; outlier detection and removal approaches using MR-PRESSO outlier test, Cook’s distance^24,25, Studentized residuals²⁴, Q (modified) outlier test or Q’ (modified) outlier test) in multi-instrument summary-level MR. For these comparisons, we performed 10,000 simulations under different scenarios using the model described in Supplementary Figure 1b. We varied several parameters in the simulations including the main causal effect (β_causal = 0, 0.1, 0.2, 0.5 with β_pleiotropic = 0.1), the percentage of horizontal pleiotropic outlier variants (0, 2, 4, 10, 50, 90%), the type of horizontal pleiotropy (positive or balanced) and whether or not we assume the Instrument Strength Independent of Direct Effect (InSIDE) condition²⁰ (Online Methods; Supplementary Note).

False positive rate/power to detect horizontal pleiotropy

Using simulations, we first assessed the false positive rate (type 1 error) and power (1 – type 2 error) of the above MR methods that detect horizontal pleiotropy, including the MR-PRESSO global test, the Q test^17,18, Q (modified) test¹⁹, Q’ test^17,18 and Q’ (modified) test¹⁹ (Online Methods; Supplementary Note).

Under the null hypothesis of no horizontal pleiotropy (Supplementary Table 1), the MR-PRESSO global test, the Q (modified) test and Q’ (modified) test had controlled false positive rates (~5%) whereas the rates of the original Q and Q’ tests were inflated (between 5 and 25%). Since the Q and Q’ tests were found to be inflated, we proceeded with the MR-PRESSO global test, the Q (modified) test and the Q’ (modified) test in power analyses. The power to detect horizontal pleiotropy was similar for the MR-PRESSO global test and the Q (modified) test (Table 1). We observed acceptable power to detect horizontal pleiotropy across all three tests under simulations when the percentage of horizontal pleiotropic variants was ≥ 10% (which corresponds to 5 horizontal pleiotropic variants out of 50). Similar results were observed across all magnitudes of causal estimates, including instances in which there was no causal relationship (β_causal = 0) (Table 1). In addition, we conducted several sensitivity analyses (violation of the InSIDE condition, percentage of horizontal pleiotropic variants ≥ 50% and perfectly overlapping samples to estimate the effect sizes on the outcome and exposure) to further evaluate the robustness of these methods. All three tests had high power to detect horizontal pleiotropy under a wide range of these parameters; however, simulations showed a reduction in power in those with perfectly overlapping samples (Supplementary Results and Supplementary Tables 2–4).

Evaluation of bias in the causal estimates

We investigated how the causal effect estimate (bias) and corresponding standard deviation (precision) of the MR-PRESSO outlier test were affected by horizontal pleiotropy. We then compared MR-PRESSO to other established MR methods that can correct for horizontal pleiotropy, including those that correct for a global average horizontal pleiotropic effect across all variants (e.g. inclusion of the intercept in MR-Egger regression^20,21, adjustment for multiple exposures in Multi-variable MR (MMR)^22,23), horizontal pleiotropic outlier detection methods (Cook’s distance^24,25, Studentized residuals²⁴, Q (modified) outlier test, Q’ (modified) outlier test) and outlier-robust methods (weighted median²⁶, mode-based estimate²⁷) (Online Methods).

The MR-PRESSO outlier test and Cook’s distance had similar power to identify the correct horizontal pleiotropic outliers whereas Studentized residuals had low power (Supplementary Table 5; Supplementary Note). However, Cook’s distance had lower specificity to identify the correct horizontal pleiotropic outliers compared to the MR-PRESSO outlier test (inflated false positive rate (family-wise error rate), > 97% when there was no horizontal pleiotropy). The false positive rate (family-wise error rate) and power of the Q (modified) outlier test and the Q’ (modified) outlier test were very similar to the MR-PRESSO outlier test (Supplementary Table 5).

MR-Egger regression generally had lower precision than the MR-PRESSO outlier test and MMR. The $I_{GX}^2$ index, which informs on when MR-Egger regression should be employed (i.e. $I_{GX}^2$>90%), was lower than 90% on average in most settings (Supplementary Table 6; Supplementary Results).

Cook’s distance and Studentized residuals had similar bias and precision in the causal estimate compared to the MR-PRESSO outlier test (Supplementary Table 7). The weighted median had less bias but also less precision in the causal estimate compared to the MR-PRESSO outlier test, particularly when the percentage of horizontal pleiotropic variants was < 50%. The mode-based estimate generally had very low precision compared to the other methods.

All four methods (mode-based estimate, weighted median, Cook’s distance and Studentized residuals) had similar limitations as the MR-PRESSO outlier test (e.g. inflated causal estimates and low precision) when there was a very high percentage of horizontal pleiotropic variants (≥ 50%) (Supplementary Table 7; Supplementary Note). Our simulations showed that when the InSIDE assumption was invalid, all methods exhibited bias in the causal estimate (Supplementary Tables 8 and 9; Supplementary Note).

The standard IVW approach showed an expected bias in the causal estimate due to horizontal pleiotropy. When the InSIDE assumption was valid and the percentage of horizontal pleiotropic variants was small (≤ 10%), the causal estimate of the MR-PRESSO outlier adjustment was less biased and had better precision (smaller standard deviation) than IVW, MR-Egger and MMR. However, when the percentage of horizontal pleiotropic variants was high (≥ 50%), the opposite was found. These trends were expected since outlier detection methods such as the MR-PRESSO outlier test require, as a baseline assumption, that at least 50% of the genetic variants to be valid instruments (no horizontal pleiotropy), have balanced pleiotropy and for the InSIDE assumption to be valid. Conversely, methods that correct for a global average horizontal pleiotropic effect amongst all variants (either by including the regression intercept in MR-Egger or covariate adjustment in MMR) are best suited when there is a large percentage of horizontal pleiotropic variants (> 50%).

Horizontal pleiotropy detection in MR for complex traits

We applied the MR-PRESSO global test, along with two methods (Q (modified) test, Q’ (modified) test) that can detect horizontal pleiotropy in MR to all possible pairs of 82 complex traits and diseases retrieved from publicly available GWA datasets (Table 2; Online Methods). In total, we conducted 4,250 tests for each of the three MR approaches. We accounted for multiple testing of the 4,250 tests using the Bonferroni correction. We note that this correction is overly stringent since many of the traits and diseases are correlated. Using a Bonferroni-corrected threshold of P < 1.17 × 10⁻⁵, the MR-PRESSO global test was statistically significant in 21.69% (n = 922) of the 4,250 tests. When restricting to statistically significant causal estimates in the IVW meta-analysis, we detected significance of the MR-PRESSO global test at a higher rate of 48.69% (n = 93 of 191 tests). Both the Q (modified) test and the Q’ (modified) test provided similar estimates. As a sensitivity analysis, we restricted to a subset of traits and diseases that were less correlated (e.g. Pearson r < 0.30). We observed that 24.17% of MR tests were statistically significant for the MR-PRESSO global test amongst significant causal relationships (Supplementary Table 10).

Horizontal pleiotropy correction in MR for complex traits

We evaluated five methods to correct for horizontal pleiotropy in MR (Table 3). This included: 1) removing outliers detected by the MR-PRESSO outlier test, Q (modified) outlier test and Q’ (modified) outlier test (at Bonferroni-corrected threshold); and 2) adjusting for significant covariates individually or all together from the main MR test (e.g. significant causal effect in IVW meta-analysis at Bonferroni-corrected threshold). As shown in the simulations, the MR-PRESSO global test can be used to determine if there is any remaining horizontal pleiotropy after applying correction strategies to minimize initial horizontal pleiotropy in the MR test (Supplementary Note; Supplementary Table 11). In Table 3, we observed that the outlier removal approach using the MR-PRESSO outlier test was effective in eliminating statistical significance in the MR-PRESSO global test in 46% of the 922 tests. The Q (modified) outlier test and Q’ (modified) outlier test provided similar estimates. Furthermore, the covariate adjustment approach – defined by accounting for traits that were shown to have a significant causal effect on the same outcome – eliminated significance in the MR-PRESSO global test in 22% (n = 20) of the 93 tests when adjusting for a single covariate. When adjusting for all significant covariates in the same model, the covariate adjustment approach eliminated significance in 34% (n = 22) of the 42 tests. Taken together, these two correction strategies (MR-PRESSO outlier removal and MMR) were successful in 47% (n = 438) of the 922 tests total. Furthermore, we note that the covariate adjustment approach is limited in that it requires a priori knowledge of the trait responsible for the horizontal pleiotropic effect.

Outliers effect on the distortion of MR causal estimates

We evaluated the extent to which outliers cause distortion in the causal estimates resulting from MR. Using the MR-PRESSO distortion test, we compared the causal estimates from the IVW meta-analysis before and after removal of outlier variants detected by the MR-PRESSO outlier test (Online Methods). Using a Bonferroni-corrected threshold, we observed a significant distortion (of −93% and 35%) in 2.5% (n = 2) of significant causal estimates (n = 81 total). Since the Bonferroni correction is overly stringent, we considered the commonly-used nominal threshold of P < 0.05 that the majority of MR studies to date have used for statistical significance. A significant distortion was observed in almost 10% (n = 22) of the causal relationships (n = 229 total) with a distortion between −131% and 201% on average (Figure 2).

Open in a separate window

Figure 2:

Distribution of the distortion of causal estimates before and after correction for horizontal pleiotropy using the MR-PRESSO distortion test.

The distortion coefficient was calculated for all exposure-outcome pairs with a significant causal estimate (P < 0.05) using an inverse variance weighted meta-analysis (n = 229). The distortion coefficients were then tested using MR-PRESSO (Mendelian Randomization Pleiotropy RESidual Sum and Outlier) distortion test which provides an empirical P-value. The distortion coefficients are colored according to whether the distortion is statistically significant (blue) or not (red) at a Bonferroni-corrected threshold of P < 0.05 / 229 in the MR-PRESSO distortion test. The distortion estimate represents the change in the causal estimate as a result of horizontal pleiotropic outlier variants (Online Methods). A positive distortion represents a decrease in the outlier-corrected causal estimate.

Below, we provide one example to highlight the role of pleiotropy in MR (Supplementary Figure 2). We observed that the causal effect of body mass index (BMI) on C-reactive protein was estimated to be 0.39 (P = 7.02 × 10⁻⁸) by IVW. The MR-PRESSO global test showed statistical significance (P < 10⁻⁶) and the MR-PRESSO outlier test identified one significant outlier variant (P < 10⁻⁶, rs2075650 in the APOE locus). Examining this further, we observed that this variant was highly pleiotropic, having associations with several traits and diseases including Alzheimer’s Disease, body mass index, C-reactive protein, high-density lipoprotein cholesterol, low-density lipoprotein cholesterol, plasma triglycerides, waist circumference, hip circumference and waist-hip ratio at P = 6 × 10⁻⁴ (Bonferroni-corrected P = 0.05 / 82; Supplementary Table 12). Furthermore, this variant was associated with several other traits and diseases in the public NHGRI-EBI GWAS catalog²⁸ (P < 5 × 10⁻⁸; Supplementary Table 13). We believe some of these genetic associations with other traits and diseases will be due to horizontal pleiotropy while others will be in the same causal pathway (e.g. vertical pleiotropy). After removing this outlier variant, we observed a lower estimation of the causal estimate for BMI on C-reactive protein (β_causal = 0.35, P = 3.45 × 10⁻¹⁶) with this single variant alone causing a 12% distortion in the causal estimate (Supplementary Figure 2).

Outliers effect on false positive causal relations in MR

We evaluated the extent to which outliers (as detected by the MR-PRESSO outlier test) can induce false positive causal relationships. A false positive causal relationship was defined as an exposure-outcome pair in which the causal estimate was no longer statistically significant in the outlier-corrected IVW model but was previously significant in the naïve IVW model. According to this definition, we identified false positive relationships in 10% of putatively causal relationships (n = 24 out of 229 total tests) using the common nominal P < 0.05 threshold and 1.2% (n = 1 out of 81 total tests) using the stringent Bonferroni-corrected threshold (P < 1.17 × 10⁻⁵). We note that the outlier removal approach decreases the power of the outlier-corrected IVW model test because a smaller number of variants is included after removal of horizontal pleiotropic outlier variants. Therefore, we expect that a small number of false positive causal relationships may be due to reduced statistical power.

Existing methods to detect and correct for horizontal pleiotropy in MR

Several existing methods have been developed to detect horizontal pleiotropy in MR. Methods to detect horizontal pleiotropy include the Q test and the Q’ test, as well as the recently proposed modified versions of these tests¹⁹, which are traditionally used to identify over dispersion and have been applied in the context of MR¹⁷.

Several existing MR methods have also been designed to correct for horizontal pleiotropy. Methods that we evaluated in the current study are described below.

MR-PRESSO: MR Pleiotropy RESidual Sum and Outlier

We developed the MR Pleiotropy RESidual Sum and Outlier (MR-PRESSO) test to detect and correct for horizontal pleiotropic outliers. MR-PRESSO is comprised of three components (Figure 1):

• a)

  detection of horizontal pleiotropy (and violation of the exclusion restriction assumption, Supplementary Figure 1) in MR (global test);

• b)

  correction by removal of offending IVs that are due to horizontal pleiotropy (outlier test);

• c)

  testing of significant differences in the causal estimates before and after outlier removal (distortion test).

MR-PRESSO extends the framework of IVW meta-analysis and is based on the following rationale. IVW meta-analysis consists of fitting a weighted linear regression between the set of effect sizes on the outcome and the effect sizes on the exposure. The slope of the regression provides an estimate of the causal effect of the exposure on the outcome. Under the null hypothesis where horizontal pleiotropy does not exist, all variants are expected to be close to the regression line (i.e. to have small residuals in the regression). However, when a variant is subject to horizontal pleiotropy, the effect size on the outcome can be larger or smaller than the effect size mediated by the exposure in question and therefore the variant can deviate from the true slope of the regression line. MR-PRESSO is designed to detect whether a subset of variants is significantly deviating from the regression line. The MR-PRESSO outlier test requires that at least 50% of the variants are valid instruments and relies on the InSIDE (Instrument Strength Independent of Direct Effect) condition²⁰, which states that the effect sizes of the variants on the exposure should not depend on the horizontal pleiotropic effects on the outcome.

MR-PRESSO global test

The MR-PRESSO global test evaluates for the presence of horizontal pleiotropy and is comprised of four steps (Figure 1a):

• 1)

  for each variant j, we removed the variant in question and refit a IVW regression. This allows us to calculate the slope of the regression line on the remaining variants, denoted ${\widehat{\beta }}_{-j}$ which represents the causal estimate without variant j;

• 2)

  the estimated causal effect (slope) without variant j is used to predict the expected effect size on the outcome as the product of${\widehat{\beta }}_{-j}$ and the effect size of the same variant on the exposure ${\widehat{\gamma }}_j$. Then, we calculated the observed residual sum of squares (RSS) as the difference between the observed effect size of the variant on the outcome (${\widehat{\Gamma }}_j$) and the predicted effect size of the same variant on the outcome $RS{S}_{obs}\left(j\right)={\left({\widehat{\Gamma }}_j-{\widehat{\beta }}_{-j}{\widehat{\gamma }}_j\right)}^2$. The global observed RSS is then obtained by summing over the J $RS{S}_{obs}\left(j\right)$: $RS{S}_{obs}={\displaystyle \sum }_{j}RS{S}_{obs}\left(j\right)={\displaystyle \sum }_j{\left({\widehat{\Gamma }}_j-{\widehat{\beta }}_{-j}{\widehat{\gamma }}_j\right)}^2$;

• 3)

  the observed RSS is compared to a simulated expected distribution of RSSs. The expected distribution is simulated under the null hypothesis (0% of variants are outliers). First, we simulated a distribution of effect sizes on the exposure ${\widehat{\gamma }}_j^{random}$ from a Gaussian distribution $\mathcal{N}\left({\widehat{\gamma }}_j,\mathbb{V}\left({\widehat{\gamma }}_j\right)\right)$. Second, we simulated a distribution of effect sizes on the outcome ${\widehat{\Gamma }}_j^{}$_random around the predicted effect size on the outcome (${\widehat{\beta }}_{-j}{\widehat{\gamma }}_j$) by drawing in a Gaussian distribution $\mathcal{N}\left({\widehat{\beta }}_{-j}{\widehat{\gamma }}_j,\mathbb{V}\left({\widehat{\Gamma }}_j\right)\right)$. The expected RSS is then obtained as ${\displaystyle \sum }_{j}RS{S}_{exp}\left(j\right)={\left({\widehat{\Gamma }}_j^{random}-{\widehat{\beta }}_{-j}{\widehat{\gamma }}_j^{random}\right)}^2$. The procedure was repeated multiple times (K) to obtain a null distribution of the K expected RSSs, $RS{S}_{exp}^k={\displaystyle \sum }_{j}RS{S}_{exp}^k\left(j\right)={\displaystyle \sum }_j{\left({\widehat{\Gamma }}_{jk}^{random}-{\widehat{\beta }}_{-j}{\widehat{\gamma }}_{jk}^{random}\right)}^2$;

• 4)

  an empirical P-value is computed as the number of expected RSSs greater than the observed RSS divided by the total number of times the procedure is repeated:$P=\frac{{{\displaystyle \sum }}_k{1}_{>RS{S}_{obs}}\left(RS{S}_{exp}^k\right)}{K}$.

K corresponds to the number of random drawings used to generate the expected distribution. The magnitude of K depends on the desired precision of the P-value which is 1/K. We recommend to use at least K = 1,000, which will produce a precision of at least 10⁻³.

MR-PRESSO outlier test

The MR-PRESSO outlier test allows for the detection of specific horizontal pleiotropic outlier variants. For a given variant j, we compared the observed j^th RSS $RS{S}_{obs}\left(j\right)$ (obtained in step 2 of the global test) with the distribution of K expected j^th RSSs $RS{S}_{exp}^k\left(j\right)$ (obtained in step 3 of the global test). Finally, an empirical P-value is computed as $P_j=\frac{{{\displaystyle \sum }}_k{1}_{>RS{S}_{obs}\left(j\right)}\left(RS{S}_{exp}^k\left(j\right)\right)}{K}$ (Figure 1b) which is then multiplied by the number of variants J to account for multiple testing using the Bonferroni correction.

Given that J outlier tests are performed, we recommend to use a Bonferroni correction of the J P-values.

MR-PRESSO distortion test

The MR-PRESSO distortion test quantifies the distortion in the causal estimate due to significant horizontal pleiotropic outlier variants (Figure 1c). Distortion (D) is defined as the percentage of the causal estimate that is due to significant horizontal pleiotropic outlier variants. It is calculated as$D=100\times \frac{{\widehat{\beta }}_{causal,o}-{\widehat{\beta }}_{causal}}{\left|{\widehat{\beta }}_{causal}\right|}$, with ${\widehat{\beta }}_{causal,o}$the original causal effect estimated using all variants and $\widehat{\beta }$_causal the corrected causal estimate obtained after removing outliers identified by MR-PRESSO. Normalizing by the absolute value of the corrected causal estimate provides a direction for the magnitude of D. To test for statistical significance of D, we calculated an empirical P-value by generating a null distribution (the null hypothesis corresponds to the expected distortion due to a random set of variants). We defined $n_O$ as the number of variants detected as outliers by the MR-PRESSO outlier test and $n_E$ as the total number of variants robustly associated with the exposure. The null distribution is generated by substituting $n_O$ variants detected as outliers by the MR-PRESSO outlier test with $n_E-2n_O$ non-outliers, which are drawn with replacement from the entire set of non-outlier variants. This results in the total number of variants being fixed at $n_E-n_O$. We repeated this procedure K times to generate the null distribution. An empirical P-value is then calculated as the number of times that the observed distortion is greater than the expected distortion under the null hypothesis divided by K.

Simulation framework

We performed simulations to evaluate the statistical properties (false positive rate and power) to detect and correct for horizontal pleiotropy. We simulated the standard MR framework shown in Supplementary Figure 1b with an outcome as well as two exposures E₁ and E₂. A total of 50 variants was simulated per case. Horizontal pleiotropy was induced by a certain percentage of the 50 variants G_j having a significant effect on both E₁ (through ${\text{γ}}_{1,j}$) and E₂ (through ${\text{γ}}_{2,j}$). We varied the following parameters:

• the main causal effect (of exposure 1 on the outcome), β_causal> = 0 (no main causal effect) or β_causal = 0.1, 0.2, 0.5 whereas the causal effect of exposure 2 was always set to β_pleiotropic= 0.1;

• percentage of horizontal pleiotropy (0, 2, 4, 10, 50, 90);

• type of horizontal pleiotropy: positive (all ${\text{γ}}_{2,j}$ positive) or balanced horizontal pleiotropy (approximately half of the ${\text{γ}}_{2,j}$ positive and half negative);

• verification of the InSIDE²⁰ (Instrument Strength Independent of Direct Effect) assumption (Supplementary Note).

For each scenario, 10,000 simulations were performed.

Collection of genome-wide association (GWA) summary statistics

We retrieved publicly available genome-wide association (GWA) summary statistics data for 82 complex traits and diseases (Supplementary Note). We performed the following steps to ensure that all datasets were uniform and standardized. For each, we retrieved the appropriate variant annotation (build, rsid, chromosome, position, reference and alternate alleles) and summary statistics (effect size, standard errors, P-values and sample size of the study). All variant coordinates (chr, pos) were lifted over to hg19 using the UCSC Genome Browser LiftOver Tool. We imputed Z-scores of variants using ImpG³⁹ using 1000 Genomes Phase 3 European panel⁴⁰ (n = 503) as a reference panel. Effect sizes, standard errors and P-values were then calculated using the variance of the trait estimated from genotyped variants and allele frequencies calculated on the same subset of individuals from the 1000 Genomes reference panel. Sets of GWA-significant variants were manually retrieved from the corresponding GWA manuscripts. In total, we retrieved GWA summary statistics for 82 traits and diseases (Supplementary Table 14).

We thank the various genome-wide association consortia for generously sharing the genome-wide association summary statistics. We thank Amanda Dobbyn and Daniel Jordan for critical comments on the manuscript. R.D is supported by R35GM124836 from the National Institute Of General Medical Sciences of the National Institutes of Health, R01HL139865 from the National Heart, Lung, Blood Institute of the National Institutes of Health, research grants from AstraZeneca and Goldfinch Bio and previously an American Heart Association Cardiovascular Genome-Phenome Discovery grant (15CVGPSD27130014). B.N. is supported by 1R01MH094469, 1R01MH107649–01 and 5U01HG009088–02 from the National Institutes of Health. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Competing Financial Interests Statement

B.N. is a member of Deep Genomics Scientific Advisory Board, has received travel expenses from Illumina, and also serves as a consultant for Avanir and Trigeminal solutions. R.D. has received research support from AstraZeneca and Goldfinch Bio.