Robust inference of bi-directional causal relationships in presence of correlated pleiotropy with GWAS summary data.

To infer a causal relationship between two traits, several correlation-based causal direction (CD) methods have been proposed with the use of SNPs as instrumental variables (IVs) based on GWAS summary data for the two traits; however, none of the existing CD methods can deal with SNPs with correlated pleiotropy. Alternatively, reciprocal Mendelian randomization (MR) can be applied, which however may perform poorly in the presence of (unknown) invalid IVs, especially for bi-directional causal relationships. In this paper, first, we propose a CD method that performs better than existing CD methods regardless of the presence of correlated pleiotropy. Second, along with a simple but yet effective IV screening rule, we propose applying a closely related and state-of-the-art MR method in reciprocal MR, showing its almost identical performance to that of the new CD method when their model assumptions hold; however, if the modeling assumptions are violated, the new CD method is expected to better control type I errors. Notably bi-directional causal relationships impose some unique challenges beyond those for uni-directional ones, and thus requiring special treatments. For example, we point out for the first time several scenarios where a bi-directional relationship, but not a uni-directional one, can unexpectedly cause the violation of some weak modeling assumptions commonly required by many robust MR methods. We also offer some numerical support and a modeling justification for the application of our new methods (and more generally MR) to binary traits. Finally we applied the proposed methods to 12 risk factors and 4 common diseases, confirming mostly well-known uni-directional causal relationships, while identifying some novel and plausible bi-directional ones such as between body mass index and type 2 diabetes (T2D), and between diastolic blood pressure and stroke.

1 Introduction

It is of great interest to infer causal relationships between pairs of complex traits or diseases such as for treatment/intervention development and drug repurposing [1, 2], which however is quite challenging and had barely been touched until recently. The availability of large-scale GWAS summary data and the use of SNPs as instrumental variables (IVs) in Mendelian randomization (MR) have made it possible for such inference [3-5]. However, most MR methods and analyses are based on a critical and strong assumption that there is only a uni-directional relationship between two traits and the direction is known; that is, by treating one trait as the exposure and the other as the outcome, one assumes that the causal relationship, if exists, can be only from the exposure to the outcome. To infer the causal direction between two traits (under the uni-directional assumption), recently several methods based on comparing correlations between SNPs/IVs and each trait have been proposed, including Steiger’s method based on a single SNP (that is assumed to be a valid IV) [6], CD-Ratio and CD-Egger based on multiple SNPs, which can be more powerful than Steiger’s method [7]. CD-Egger, similar to Egger regression in MR [8], is also more robust than the other two methods by allowing invalid IVs under the InSIDE assumption; that is, CD-Egger allows invalid IVs with uncorrelated pleiotropy, but not correlated pleiotropy [9]. The first goal here is to develop a correlation-based causal direction (CD) inference method based on constrained maximum likelihood, called CD-cML, and show its higher power and robustness than the above methods, especially in the presence of correlated pleiotropy. Given the wide-spread pleiotropy [10, 11], it is of utmost importance for any method to be robust to pleiotropy, especially correlated pleiotropy that is more challenging to deal with. However, the above CD methods are applicable to infer only uni-directional, but not bi-directional, causal relationships. In a bi-directional relationship, each of the two traits may be causal to the other at the same time. In practice, there may be bi-directional causal relations between some traits (e.g. between insomnia and some major psychiatric disorders [12]), or at least we may not be able to exclude a priori the possibility of such bi-directional relationships.

Alternatively, reciprocal (also called bidirectional) MR can be applied by treating each of the two traits as the exposure while the other as the outcome [13, 14]. However, as shown in [7, 15], bidrectional MR (with the use of many MR methods) does not perform well due to some reasons, including the below one: assuming the true causal direction is from X to Y for two traits X and Y, if an SNP is causal to X (and the sample sizes are large enough), the SNP is associated with both X and Y, and thus may be considered as a candidate IV for both traits; when the SNP is used as an IV for direction X to Y, it will confirm the causal association; however, if it is used as an IV for direction Y to X, it will also yield a non-zero estimate of the causal effect of Y on X, leading to an incorrect conclusion. A naive remedy is to remove any SNP associated with both traits, but it leads to not only loss of power (with fewer SNPs as IVs), but also biased inference (e.g., towards the causal direction X to Y if the truth is Y to X and if the GWAS sample size or power for X is much larger than for Y) [16]. Here we adopt a simple but effective screening/filtering rule based on a simple heuristic: no SNP will be used as an IV for both traits, because no SNP can be valid for both traits. If an SNP is associated with both traits, by Steiger’s method, we use it only for the trait with which its absolute correlation is larger than that with the other trait (because it is more likely to be a valid IV for the chosen trait) [16]. Furthermore, there are some new MR methods, such as constrained maximum likelihood (MR-cML) [17], that are more robust to both uncorrelated and correlated pleiotropy. Our second goal here is to show that, by incorporating MR-cML and the IV screening rule in reciprocal MR, the resulting method, still called MR-cML for simplicity, performs well, in fact almost identically to CD-cML if their modeling assumptions hold; otherwise, CD-cML controls type I error better and is more conservative.

With the two robust and powerful methods, we show their application to infer bi-directional relationships, which has been largely neglected in the literature. It is notable that inferring bi-directional causal relationships is far more challenging than uni-directional ones: for example, for the first time we point out that a bi-directional causal relationship generates a few new scenarios, in which either the InSIDE assumption or the plurality condition required by many existing robust MR methods will be violated (e.g. IVW (random effect), Egger regression [8] and RAPS [18] for the former; our cML methods, MR-ContMix [19], MR-Mix [20], MR-Lasso [21] and MR-Weighted Mode [22] for the latter). We applied the methods to 48 risk factor-complex disease pairs with 12 cardiometabolic risk factors, 3 cardiometabolic diseases (T2D, Stroke and CAD), and asthma (more as a negative control), identifying some interesting bi-directional causal relationships, such as between diastolic blood pressure and Stroke, and between body mass index and T2D [23].

2 Results

2.1 Overview of methods

Given two traits, X and Y, and two independent GWAS datasets for the two traits, our goal is to infer their possibly bi-directional causal relationship. One of the most challenging issues is that we have a hidden (i.e. unobserved) confounder (or equivalently, an aggregate of many hidden confounders) denoted by U, which is associated with both X and Y with effect sizes θ and θ respectively. There are three possible sets of candidate SNPs to be used as IVs: (1) {g} is the set of valid IVs for X, having direct effects α’s only on X; (2) {g} is the set of valid IVs for Y, having direct effects β’s only on Y; (3) {g} is the set of invalid IVs directly influencing both X and Y, and possibly U, with direct effects γ’, η’s, and ξ’s respectively. Fig 1 illustrates a true causal model, in which θ and θ are the causal effects between the two traits, the unknown parameters of interest.

Fig 1

The true causal model with two traits X and Y of interest.

U is a hidden confounder (or an aggregate of hidden confounders). The IVs in {g} and {g} are valid for X and Y respectively, and the IVs in {g} are invalid ones. The true classification of {g}, {g} and {g} is unknown and needs to be estimated. The arrows give the directions of the direct causal effects (with the effect sizes shown). In particular, θ and θ are the causal effects from X to Y and from Y to X respectively, and are parameters of interest.

Define the (population) Pearson correlations between each candidate SNP/IV g and two traits as ρ = corr(X, g) and ρ = corr(Y, g). It is shown in Methods that with and .

If there exists causal direction from X to Y, we have |K| < 1 (under a suitable condition shown in Methods) and K ≠ 0, which can be used to infer the causal direction of X to Y. Similarly, we use |K| < 1 and K ≠ 0 to infer the causal direction of Y to X. The idea is similar to that used by other CD methods, i.e. Steiger’s method based on a single valid IV, CD-Ratio on multiple valid IVs and CD-Egger on multiple possibly invalid IVs without correlated pleiotropy (i.e. when the InSIDE assumption holds) [6, 7].

Since K and K are unknown, we propose a constrained maximum likelihood, called CD-cML, to infer the two parameters and thus the causal direction. Briefly, based on the given GWAS (summary) data, we calculate the sample (Pearson) correlations between each candidate SNP/IV and each trait, say r and r, which are asymptotically normal and consistent for the (population) correlations ρ and ρ; with Eq (1), we can write down the normal-based log-likelihood under the constraint that the number of invalid IVs is equal to a given integer, say m ≥ 0. We try each possible value of m, then consistently select the best one based on the Bayesian Information Criterion (BIC). The resulting constrained maximum likelihood estimates (cMLEs), say and , are consistent for K and K respectively, and are asymptotically normal. Hence we can construct a normal-based confidence interval for K and K respectively, thus drawing inference on the two possible causal directions from X to Y and from Y to X.

A similar method, called MR-cML, has been used to estimate θ or θ in the framework of MR [17]. It is noted that MR-cML performs well under correlated pleiotropy. Here we also propose applying MR-cML to reciprocal MR to infer both θ and θ, and thus infer a possibly bi-directional causal relationship between X and Y; for simplicity, we still call such a reciprocal MR as MR-cML.

It is noted that MR and CD methods are related but different: For example, for direction of X to Y, MR methods are based on inferring whether θ = 0; in contrast, CD methods are based on whether both K = 0 and |K| < 1. Accordingly, because of the second constraint, we expect that sometimes CD-cML will be more conservative than MR-cML in terms of yielding smaller type I error and lower power.

We also propose a simple but yet effective method for SNP/IV screening based on a simple heuristic: none of SNPs can be a valid IV for both X and Y. Thus, if an SNP is (marginally) associated with both traits, we will use it only for the trait with which it is more correlated than with the other trait. This screening rule is just a simple application of Steiger’s method [6], and was mentioned in [16]. This screening rule is especially useful in the presence of bi-directional causal relationships. For example, when inferring whether there is a causal direction from X to Y in Fig 1, if θ ≠ 0, all IVs in set {g} are associated with trait X and thus are candidate IVs, though they are all invalid IVs; the screening rule will eliminate them as IVs (because they will be more highly correlated with trait Y than with trait X; see Methods). Now we consider what happens if they are indeed used as IVs. Assuming the set size |{g}| is larger than that of the valid IV set, |{g}|, the invalid IV set {g} forms the largest (i.e. plurality) group in (incorrectly) estimating θ as 1/θ (asymptotically); in other words, they lead to the violation of the plurality condition required by cML (and several other MR methods, such as MR-ContMix [19], MR-Mix [20], MR-Lasso [21] and MR-Weighted Mode [22]). In addition, they will also lead to the violation of the InSIDE assumption: it is easy to verify that for any g ∈ {g}, its effect size on trait X is β = θ β, which is clearly correlated with β, its direct effect on Y.

There is another source leading to the violation of the InSIDE assumption, in addition to the more widely recognized one with ξ ≠ 0 (i.e. some IVs are correlated with the hidden confounder) and the one pointed above. It is due to bi-directional causal relationships, again demonstrating that it is more challenging to deal with bi-directional relationships than with uni-directional ones. Consider the causal direction of X to Y: even if ξ = 0 but if θ ≠ 0, any SNP g ∈ {g} would lead to the violation of InSIDE, because its association strength with X and its direct effect size on Y respectively would be γ + ηθ and η, which are clearly correlated. The violation of the InSIDE assumption will lead to biased inference by several popular random-effects model-based methods (that treat direct effects as random effects), such as IVW (random effect), Egger regression [8] and RAPS [18].

Finally we apply the data perturbation (DP) scheme of [17] for better finite-sample inference: it accounts for uncertainty in selecting invalid IVs in CD- and MR-cML, leading to better control of type I errors. We suffix a method with “-DP” and “-S” to refer its use of data perturbation and IV screening respectively. See Section 4.6 for a summary of different methods.

2.2 Simulations

2.2.1 Main simulations

We generated simulated data following the true causal model in Fig 1 for two continuous traits X and Y. There were 15, 10 and 10 SNPs/IVs in sets {g}, {g} and {g}, respectively, with effect sizes α1 to α15, β1 to β10, γ1 to γ10, and η1 to η10 ranging in (−0.3, −0.2) and (0.2, 0.3) (from the corresponding uniform distributions) respectively. For more correlated pleiotropy, we generated ξ’s from a uniform distribution in the range of (−0.2, 0.2); otherwise, we set ξ’s at 0. The MAFs of the SNPs were 0.3. The random errors ϵ and ϵ were independently drawn from N(0, 1), and ϵ was from N(0, 2). We considered various combinations of the true causal effect sizes of θ and θ in the range from 0 to 0.3. For each set-up, we generated 500 pairs of two independent GWAS samples, one for each of the two traits and each of sample size n = N1 = N2 = 50000. We also studied the scenarios with at least one of X and Y being binary. To generate binary traits, we generated the continuous X and Y first, then dichotomized one or both of them by setting the largest 30% of their values to be 1 and the other 70% as 0. For each dataset, we first generated individual level data, then calculated the summary statistics with marginal linear regression for a continuous trait and marginal logistic regression for a binary trait. We set the significance cutoff at 0.05/35 to select relevant IVs for both traits before applying any CD and MR methods.

We summarize the main results in terms of (empirical) type I error and power in Figs 2–5. Figs 2 and 4 show the results for both X and Y being continuous, while Figs 3 and 5 for both X and Y being binary; the results were similar regardless of the traits being continuous or binary, though it was slightly more powerful to use the continuous traits than the binary traits (due to the loss of information by dichotomizing a continuous trait). The top-left panels (for “θ = 0, X to Y”) show (empirical) type I error for the direction of X → Y; when θ = 0 in the right panels, it shows type I error for the direction Y → X; otherwise, it is for (empirical) power. In general, MR-cML-DP-S and CD-cML-DP-S performed almost identically: they could control type I error satisfactorily while having high power. In contrast, all other methods, namely CD-Ratio, CD-Egger and combining (single-SNP-based) Steiger’s method over multiple IVs (by majority voting, “-MV”), could not control type I error, and might have low power. Interestingly, CD-Egger had largely inflated type I error rates even when none of the IVs were correlated with the hidden confounder (i.e. ξ = 0) unless both θ = θ = 0 as shown in Figs 2 and 3. This might sound surprising, but convincingly showed that detecting bi-directional causal relationships is much more challenging than detecting uni-directional ones. As explained in Methods, even if ξ = 0, in the presence of correlated pleiotropy the InSIDE assumption required by Egger regression could be violated. On the other hand, when the SNPs from {g} were correlated with the hidden confounder (with ξ ≠ 0), the InSIDE assumption would always be violated, leading to inflated type I error of CD-Egger as shown in Figs 4 and 5. Notably, both MR-cML-DP-S and CD-cML-DP-S did not suffer from any of these problems.

Fig 2

Empirical type-I error and power (y-axis) for both X and Y continuous, ξ = 0 (i.e. no correlated pleiotropy), θ = 0 (top panels) and θ = 0.3 (bottom) for various values of θ (x-axis).

The left panels show results for the direction X → Y, the right ones show results for the direction Y → X.

Fig 5

Empirical type-I error and power (y-axis) for both X and Y binary, ξ from a uniform distribution (i.e. with correlated pleiotropy, implying InSIDE violated), θ = 0 (top panels) and θ = 0.3 (bottom) for various values of θ (x-axis).

The left panels show results for the direction X → Y, the right ones show results for the direction Y → X.

Fig 4

Empirical type-I error and power (y-axis) for both X and Y continuous, ξ from a uniform distribution (i.e. with correlated pleiotropy, implying InSIDE violated), θ = 0 (top panels) and θ = 0.3 (bottom) for various values of θ (x-axis).

The left panels show results for the direction X → Y, the right ones show results for the direction Y → X.

Fig 3

Empirical type-I error and power (y-axis) for both X and Y binary, ξ = 0 (i.e. no correlated pleiotropy), θ = 0 (top panels) and θ = 0.3 (bottom) for various values of θ (x-axis).

The left panels show results for the direction X → Y, the right ones show results for the direction Y → X.

Empirical distributions of the estimates of the correlation ratio K for CD methods and the causal effect θ for MR method with both X and Y continuous, true θ = 0 and θ = 0.3.

Figs 6–9 show the empirical distributions of the parameter estimates from the methods under various simulation set-ups for both X and Y being continuous. It is confirmed that both MR-cML-DP-S and CD-cML-DP-S always gave (almost) unbiased estimates of the causal parameter (θ) and the correlation ratio parameter K) respectively. In contrast, CD-Egger was biased except when the InSIDE assumption held (i.e. ξ = 0 and no causal effect from the other direction), while CD-Ratio was in general biased in the presence of invalid IVs. Furthermore, CD-cML-DP-S and MR-cML-DP-S were much more efficient in yielding estimates with much smaller variances than those of CD-Ratio and CD-Egger.

Fig 6

Empirical distributions of the estimates of the correlation ratio K for CD methods and the causal effect θ for MR method with both X and Y continuous, true θ = 0 and θ = 0.

The top and bottom panels show the results for ξ = 0 (i.e. no correlated pleiotropy) and ξ from a uniform distribution (i.e. with correlated pleiotropy, implying InSIDE violated), respectively. The left panels show the estimates for the causal direction of X → Y, while the right ones show that for Y → X. The horizontal dashed lines are for the true values of K or θ. The two rows of the numbers under each panel give the sample mean and standard deviation of the estimates from each method.

Fig 9

Empirical distributions of the estimates of the correlation ratio K for CD methods and the causal effect θ for MR method with both X and Y continuous, true θ = 0.3 and θ = 0.3.

The top and bottom panels show the results for ξ = 0 (i.e. no correlated pleiotropy) and ξ from a uniform distribution (i.e. with correlated pleiotropy, implying InSIDE violated), respectively. The left panels show the estimates for the causal direction of X → Y, while the right ones show that for Y → X. The horizontal dashed lines are for the true values of K (in blue) and θ (in black). The two rows of the numbers under each panel give the sample mean and standard deviation of the estimates from each method.

We also illustrate the significant role of the IV screening rule: without it both MR-cML-DP and CD-cML-DP could be biased. As shown in Figs 8 and 9, when considering the direction of Y to X with θ ≠ 0, without screening rule, some SNPs would be used as invalid IVs, thus estimating θ as 1/θ incorrectly. It is also noted that, in these situations, the estimated K’s were much larger than 1, leading to much smaller type I error (and lower power) of CD-cML-DP than MR-cML-DP.

Fig 8

Empirical distributions of the estimates of the correlation ratio K for CD methods and the causal effect θ for MR method with both X and Y continuous, true θ = 0.3 and θ = 0.

The top and bottom panels show the results for ξ = 0 (i.e. no correlated pleiotropy) and ξ from a uniform distribution (i.e. with correlated pleiotropy, implying InSIDE violated), respectively. The left panels show the estimates for the causal direction of X → Y, while the right ones show that for Y → X. The horizontal dashed lines are for the true values of K (in blue) and θ (in black). The two rows of the numbers under each panel give the sample mean and standard deviation of the estimates from each method.

More results for the methods with or without data perturbation and/or IV screening are provided in the S1 Text.

2.2.2 Secondary simulations

The recently proposed Latent Heritable Confounder MR (LHC-MR) method [24] aims to infer bi-directional causal relationships as well, so we compare our proposed methods with LHC-MR through simulations. Since LHC-MR requires the use of genome-wide GWAS summary data, not just those for significant SNPs as required by our and most other methods, we cannot apply LHC-MR in our previous simulations. Instead, we generated new simulated data for LHC-MR as described in the original LHC-MR paper [24]. In each simulation, we generated 50000 LD-independent SNPs, and two independent samples each of size n = 50000. We set polygenicity parameters as π = 0.0015, π = 0.002, π = 0.0025, and set variances of non-zero direct effects as , , . Then for i = 1, ⋯, 50000, we generated the effects from each SNP to X, Y and U independently as in Equations (14), (15), (16) in [24],

We set the effects from U to X and Y as θ = 0.3 and θ = 0.25 respectively, then calculated the total effects of the SNPs on X and Y as in Equations (10) and (11) in [24], denoted by β and β as

Instead of generating GWAS individual level data, following [20], we directly generated GWAS summary statistics and as and . Then we applied LHC-MR and our proposed methods to simulated data. We tried 4 different combinations (θ, θ) ∈ {(0,0), (0,0.3), (0.3,0), (0.3,0.3)}; for each combination we did 200 simulations (due to the slow speed of LHC-MR: as a comparison, in [24] the simulations were replicated only 50 times for each setup).

Table 1 gives the simulation results: in each cell we show the proportion of the significant results with p-value less than 0.05; when θ (or θ) is 0, it shows the empirical type-I error for X to Y (or Y to X); when θ (or θ) is 0.3, it shows the empirical power for X to Y (or Y to X). We could see that MR-cML-DP-S and CD-cML-DP-S again performed almost identically: they were a little conservative with the type-I error rates smaller than the nominal level 0.05; however, they were less conservative and much more powerful than LHC-MR.

Table 1

Empirical type I error and power in the secondary simulations comparing MR-cML-DP-S, CD-cML-DP-S and LHC-MR.

(θXY, θYX)Methods	MR-cML-DP-S		CD-cML-DP-S		LHC-MR
	X → Y	Y → X	X → Y	Y → X	X → Y	Y → X
(0,0)	0.01	0.015	0.015	0.015	0	0
(0,0.3)	0.015	1	0.015	1	0	0.15
(0.3,0)	1	0.02	1	0.02	0.115	0
(0.3,0.3)	1	1	1	1	0.3	0.165

2.3 Real data examples

We studied possibly bi-directional causal relationships between 12 cardiometabolic risk factors and 4 common diseases, and between pairs of the 4 common diseases, including 3 cardiometabolic diseases, namely coronary artery disease (CAD) [25], Stroke [26], type 2 diabetes (T2D) [27], and asthma (largely serving as a negative control) [28]. The 12 risk factors are triglycerides (TG), low-density lipoprotein cholesterol (LDL), hight-density lipoprotein cholesterol (HDL) [29], Height [30], body-mass index (BMI) [31], body fat (BF) [32], birth weight (BW) [33], diastolic blood pressure (DBP), systolic blood pressure (SBP) [34], fasting glucose (FG) [35], Smoke and Alcohol [36]. Based on the existing literature, [9] partitioned the 48 risk factor-disease pairs into four groups, representing likely “causal”, “correlated”, “unrelated” and “non-causal” relationships.

We used R package TwoSampleMR to extract and pre-process the GWAS summary statistics. For each pair of traits X and Y, at the p-value cutoff 5 × 10−8, we first extracted all SNPs significant with X, denoted by I, and all SNPs significant with Y, denoted by I. Let I = I ∪ I be the combined set of the significant SNPs for at least one of the two traits, then we extracted the summary statistics for all SNPs in I. For each SNP in I, we defined its combined p-value across the two traits as p = p ⋅ p, where p and p were the p-values for the two traits. We applied function clump_data for clumping, using its default setting with distance 10000kb, r2 = 0.001, the European population as the reference panel, and p’s as p-values for the SNPs. After clumping, we obtained a set of approximately independent SNPs, for which we had their GWAS summary statistics , SE() and , SE().

We input the summary statistics into each CD or MR method, drawing conclusions according to the decision rules (see Methods). For cML, we used 10 random starting points to find cMLEs. For data perturbation, we set the number of perturbations T = 100. For all methods except Steiger’s method, for both directions, we used the Bonferroni-corrected significance level 0.05/96 ≈ 0.0005 to construct confidence intervals. For Steiger’s method we show the results as (Proportion, Majority Vote). For example, as shown in S1 Text, for TG to CAD 50.5% SNPs were significant, for CAD to TG 18.4% SNPs were significant, and the rest 31.1% SNPs were not significant for either direction; Steiger-MV would conclude that TG to CAD as TRUE, and CAD to TG as FALSE. See S1 Text for detailed results.

For LHC-MR, we downloaded the whole genome GWAS summary data for the 16 traits from IEU GWAS database [37], which are the same as the data included in R package TwoSampleMR, ensuring a fair comparison. Then we applied LHC-MR with R package lhcMR as described in [24] to make inference about bi-directional causal relationships between any pairs of traits.

As shown in Fig 10, many of the well-accepted causal relationships from a risk factor to a disease were confirmed by most of the methods, such as from TG to CAD, and LDL to CAD. There were also some interesting findings about causal effects from diseases to risk factors, for example, T2D to BMI, and T2D to FG identified by MR-cML-DP-S, CD-cML-DP-S, CD-Ratio-S and LHC-MR; and Stroke to DBP by MR-cML-DP-S, CD-cML-DP-S and LHC-MR. Between MR-cML-DP-S and CD-cML-DP-S, they gave mostly consistent results except for two pairs: MR-cML-DP-S, but not CD-cML-DP-S, identified two well-accepted causal relationships from BMI and FG to T2D. In both cases, the two methods gave the relatively wider CIs with that of CD-cML-DP-S covering K = 1, indicating perhaps its lack of power and/or its being more conservative as discussed earlier. Finally, both CD-Ratio-S and LHC-MR suggested a few more causal relationships, including some from a disease to a risk factor, which however need to be further investigated, especially for those given by the former due to its strong assumption of all valid IVs.

Fig 10

Causal relationship of 48 risk factor-disease pairs.

As shown in the legends, the size, color and direction of a triangle represent the statistical significance, method being used, and the direction/sign of the estimated effect. In each cell, the first row is for the causal direction from the risk factor to the disease, and the second row is for the reverse; the background color represents the classification of the relationship from the risk factor to the disease.

Following [9], we can classify each of the 48 risk factor-disease pairs (from risk factor to disease) into four groups, including one as established causal relations. Fig 11 shows the numbers of significant causal effects from risk factors to diseases detected by different methods at the p-value cutoff 0.0005. It is clear that LHC-MR detected most (13) causal relations as well as incorrectly including more non-causal relations (4) as causal; CD-Ratio-S included the same number of non-causal relations as LHC-MR, but fewer true causal relations (10); in contrast, our proposed MR-cML-DP-S detected second most (11) true causal relations with few non-causal ones (2); CD-cML-DP-S detected fewer (9) true causal relations but the same number (2) of non-causal ones. On the other hand, perhaps due to its low power, CD-Egger-S detected the smallest numbers of both causal (8) and non-causal (1) ones.

Fig 11

Numbers of significant causal effects from risk factors to diseases detected by different methods at p-value cutoff 0.0005.

In the S1 Text we show the detailed results between any two of the four diseases. Based on the Bonferroni-adjusted significance cutoff of 0.05/12 ≈ 0.004, all five methods identified two causal relationships: from CAD to Stroke, and from T2D to CAD. All methods except LHC-MR identified a causal relationship from T2D to Stroke [23]. On the other hand, CD-Ratio-S, CD-Egger-S and LHC-MR suggested a reverse causation from Stroke to CAD, which might be questionable. Finally, we note that MR-cML-DP-S and CD-cML-DP-S yielded consistent results.

3 Discussion

Inference for bi-directional causal relationships is far more challenging than for uni-directional ones, which has been the focus in the literature. In particular, in most MR analyses, a uni-directional causal relationship is assumed to be known and thus pursued. There are a few exceptions, though many considered bi-directional ones under the unrealistic assumption of having only valid IVs [7, 38]. The new LHC-MR method of [24] aimed for the same problem as considered here. However, there are some major differences. Foremost, although the two true causal models are similar, the technical approaches are substantially different: they estimate all the parameters related to various direct and indirect effect sizes among IVs, the hidden confounder and the two traits in the model, requiring much stronger modeling assumptions, such as the normality assumption on the distribution of the effect sizes, leading to not only a more complex estimation procedure (with a Bayesian approach), but also some non-identifiability for parameter estimation. In contrast, except for the few key parameters of interest, such as the causal effect sizes, we treat the majority of the other parameters as nuisance ones that are combined into two intercept parameters in our model (1), leading to a much simpler estimation procedure without any distributional assumption on these nuisance parameters. In addition, LHC-MR requires the availability of genome-wide GWAS summary data, while cML requiring only that for the subset of significant SNPs. Although we have shown some numerical comparisons, more studies are needed to investigate their relative performance in practice. We also note that the major assumption in our cML method is the “plurailty condition” (that the valid IVs form the largest/plurality group in estimating the same causal parameter as their estimand), which is relatively weak because it does not even require the majority of the IVs to be valid. Nonetheless, for bi-directional causal relationships, it will always be violated for one trait if no screening is applied: since the valid IVs for one trait are invalid for the other trait, we have either an equal or even a larger number of invalid IVs to or than the number of valid IVs; the former set of invalid IVs correspond to the same (incorrect) causal estimand, implying that the set of valid IVs will not form the largest group in estimating the same (correct) causal estimand; that is, the plurality condition will be violated. In addition, this problem will be catastrophic to other methods based on modeling direct effects of invalid IVs as random effects, such as Egger regression, IVW-RE and RAPS, because it leads to the violation of the InSIDE assumption required by these methods. The simple screening rule based on the Steiger’s method is surprisingly effective in eliminating (or at least alleviating) such a problem.

Although we have mainly focused on the robustness of the cML methods to correlated and uncorrelated pleiotropy (i.e. violation of IV Assumptions 2 and 3), as shown in Theorem 1 and its following discussion, we also allow the presence of some invalid IVs that are irrelevant (with the IV Assumption 1 violated) as long as the plurality condition holds. Hence, the cML methods are robust to a large extent to the violations of all three IV assumptions. In addition, we have shown that the proposed cML methods performed well for both quantitative traits and binary traits. In particular, we have proposed a general framework justifying the application of both CD and MR methods to binary traits as both exposures and outcomes.

The recent popularity of Steiger’s method may suggest that correlation-based CD methods should be more effective than MR in determining causal directions. Surprisingly, as shown here, equipped with the same screening rule and the same robust cML estimation method, we do not see better performance of CD-cML over MR-cML when their modeling assumptions, mainly of the plurality condition, hold. A possible explanation is the following: CD methods exploit the key relationship characterizing the correlations between each trait and each SNP/IV as defined by the K parameter, which however is derived from the same causal model used by MR. In fact, since K is proportional to the causal parameter θ, according to the invariance property and the asymptotic optimality of the maximum likelihood estimator, there should be no advantage in estimating either one over the other, thus their large-sample inference and the corresponding conclusions should be (asymptotically) equivalent if their modeling assumptions hold. However, since CD methods exploit the key condition |K| < 1, they are less likely to make type I errors (albeit with lower power) than MR when their modeling assumptions are violated, as shown in our simulations (Figs 8 and 9). Therefore, as evidenced in our real data application, given no guarantee of the validity of all modeling assumptions, including the plurality condition, for any given problem, CD-cML can serve as a more conservative alternative to MR-cML.

There are a few limitations with the current study. First, we have only considered the two-sample design with two independent GWAS datasets for the two traits. It will be useful to extend the methods to the case with overlapping individuals in the two GWAS datasets, or to the one-sample design with only one GWAS dataset for both the traits [9, 24]. Second, in our real data examples the SNPs/IVs were selected from the same two GWAS datasets for the two traits, which might introduce selection bias. Alternatively, we can use an independent GWAS dataset for SNP selection for each trait [39], or analytically account for selection effects [40-42]. Third, as usual in MR we only considered linear models, while non-linear modeling may be more flexible and gain power [43, 44]. Finally, and most importantly, any analysis method comes with its modeling assumptions, especially in causal inference with observational data. Although we feel that our modeling assumptions for our proposed CD- and MR-cML, mainly the plurality condition, are relatively weak, they may or may not hold in practice. It would be useful to develop model checking techniques, or apply alternative methods under different modeling assumptions for triangulation [45]. Only more real data applications can shed light on how likely these modeling assumptions, especially the plurality condition, would hold in practice.

4 Methods

4.1 Model

The true causal model depicted in Fig 1 can be expressed as where the random errors ϵ, ϵ and ϵ are independent with each other, and for simplicity of notation we omit the intercepts (which are used in practice). The reduced form of the true model is

First, for an IV g in {g}, its (population) correlations with X and Y are

We have

Under the key assumption , we have |K| < 1. More discussions on why this assumption is often reasonable and how to empirically check this assumption are given in [7].

Second, for an IV g in {g}, its correlations with X and Y are

We have

Again under the key assumption , we have |K| < 1.

Third, for an IV g in {g}, its correlations with X and Y are

We have where K and K are defined as before, and satisfying |K| < 1 and |K| < 1 under each of the two key assumptions respectively.

Note that ξ ≠ 0 induces correlations between ρ and b, and between ρ and b, through the shared term ξ, leading to the violation of the InSIDE assumption, as expected. However, even if ξ = 0, that θ ≠ 0 or θ ≠ 0 would still induce a correlation between ρ and b (through the shared term η), or between ρ and b (through shared γ), respectively, again leading to the violation of InSIDE. Furthermore, for g ∈ {g}, if θ ≠ 0, from Eq (9) we have which means that if g ∈ {g} is used as an IV to infer the causal direction from X to Y, it would incorrectly estimate K as 1/K. In particular, if the set size |{g}| is larger than that of the valid IV set, |{g}|, it implies that the plurality condition will be violated when the SNPs in {g} are used as candidate IVs; this problem, completely due to considering a bi-directional relationship (i.e. when considering direction X to Y while allowing a causal direction of Y to X), can be avoided or alleviated by applying the IV screening rule to be discussed later. In summary, these various scenarios showcase unique challenges with analysis of bi-directional relationships.

Eq (12) can be rewritten as

Hence, combining all possible cases for g ∈ {g} ∪ {g} ∪ {g}, Eq (11) holds and can serve as a general model for statistical estimation.

Finally, it is noted that K ≠ 0 (or K ≠ 0) if and only if θ ≠ 0 (or θ ≠ 0). Accordingly we can infer causal directions based on whether K and K are non-zero and whether their absolute values are less than 1. Similarly, based on whether θ and θ are 0 or not, we can infer causal directions. Since these are unknown parameters, next we propose how to estimate them.

4.2 Constrained maximum likelihood

From a GWAS summary dataset with sample size N1 for trait X, for each g in {g} ∪ {g} ∪ {g}, we have its estimated effect size with standard error SE(); as discussed in [7], we calculate the sample correlation r between X and g as asymptotically we have , and the standard deviation σ is estimated as

Similarly, for trait Y, we obtain its sample correlation and SE(r).

For causal direction X → Y, based on the consistency and asymptotic normality of the sample correlations r’s and approximating their true standard deviations by SE(r)’s, with general model (14), we write down the log-likelihood as

We apply the constrained maximum likelihood method [17]: where m is a specified number of invalid IVs with pleiotropy (i.e. non-zero b). Since m is unknown, we try , where m is the total number of the IVs being used, then use the Bayesian Information Criterion (BIC) to select the best m. The BIC with m is where n can be an any integer between N1 and N2, the sample sizes for the two GWAS for X and Y respectively, though we recommend using n = min(N1, N2). We select , and estimate the set of invalid IVs as .

Under the plurality condition (i.e. that the valid IVs form the largest group giving the same (asymptotic) estimate of K) and that the two GWAS sample sizes are comparable, as in [17], we can consistently select m (as the true number of invalid IVs), and the resulting constrained maximum likelihood estimate (cMLE) is consistent for the true value of K and (asymptotically) normally distributed, as shown below.

Assumption 1. (Plurality condition.) Suppose that is the index set of the invalid IVs with pleiotropy for direction X → Y, i.e. with b ≠ 0 if and only if ; . For any B ⊆ {1, ⋯, m} and , if , then there does not exist any constant such that for all g ∈ B.

Assumption 2. (Orders of the variances and sample sizes.) There exist positive constants l, l, l and u, u, u such that we have l/N1 ≤ SE(r)2 ≤ u/N1, l/N2 ≤ SE(r)2 ≤ u/N2, and l ⋅ N2 ≤ N1 ≤ u ⋅ N2 for g = 1, ⋯, m.

Theorem 1 Under Assumptions 1 and 2, if , we have and as N1, N2 → ∞. Furthermore, the cMLE is consistent and asymptotically normal: where

In practice, we can consistently estimate V with

The proof for the selection consistency of BIC is similar to that for MR-cML [17], while that for the consistency and asymptotic normality of the cMLE parallels the proof of Theorem 3.3 in [18], for which the conditions are satisfied, including that here we consider only a fixed/finite m. See S1 Text for proof of Theorem 1. As in [17], we can also use the information matrix to consistently estimate the variance of ; this alternative estimator for V was used throughout this paper. At the end, based on the asymptotic normality, we can construct a normal confidence interval (CI) for K (or calculate a p-value). We can similarly estimate K and draw inference.

Note that Theorem 1 states only consistent selection of the invalid IVs with pleiotropy (i.e. with b ≠ 0) violating IV Assumptions 2 and/or 3, but not the invalid IVs that are irrelevant (i.e. violating IV Assumption 1). As discussed in Section 3.4.2 of [18], including irrelevant IVs will not affect the validity of the asymptotic normality of , but it will decrease the estimation efficiency with the variance of increased.

We note that our proposed the cML method, called CD-cML here, is similar to MR-cML in [17]. The main difference is that in MR-cML, the effect size estimates and replace the sample correlations r and r to estimate the causal effect θ or θ, instead of K or K; all other aspects remain the same.

4.3 Data perturbation

Similar to extending MR-cML to MR-cML-DP in [17], we apply CD-cML with data-perturbation to account for invalid IV selection uncertainties with finite sample sizes of GWAS datasets and weak pleiotropic effects. For t = 1, ⋯, T, we generate perturbed samples and for g ∈ {g} ∪ {g} ∪ {g}. Here T is the number of perturbations, we suggest setting it to be at least 100. With perturbed data we solve the constrained problem similar to Eq (16) as and get the corresponding maximum likelihood as

Then we average over the T perturbed estimates to get and estimate standard error of as

We construct and select with the smallest BIC(m), obtaining the corresponding DP estimate together with , then make inference about K. Similarly. we apply data-perturbation and make inference about K. This method is called CD-cML-DP.

Next we show that the proposed data perturbation scheme is consistent for CD-cML; we have a similar result for MR-cML as shown in the S1 Text. The technical details and proof are relegated to the S1 Text.

Theorem 2. Under Assumptions 1 and 2, conditional on the original GWAS summary data, as N1, N2 → ∞.

4.4 IV screening

For direction X → Y we identify the initial set of the significant SNPs associated with X, denoted by I; for direction Y → X we identify the initial set of the significant SNPs with Y, denoted by I. Then for each SNP in the intersection, g ∈ I ∩ I, if |r| ≥ |r| we keep SNP g in I and remove it from I; if |r| < |r| we keep it in I and remove it from I.

4.5 Other methods and decision rules

We will compare the proposed CD-cML with MR-cML and other CD methods: Steiger’s method, CD-Ratio (like MR-IVW) and CD-Egger (like MR-Egger) [7]. Since Steiger’s method is based on a single SNP/IV, we propose a majority voting (MV) method to combine its results across multiple SNPs/IVs. Each IV would conclude with one of three possible results: (1) no causal relationship between X and Y; (2) X has a causal effect on Y; (3) Y has a causal effect on X. We would go with the conclusion reached by the majority of the IVs (and randomly break the ties if any), called Steiger-MV. It is clear that Steiger-MV cannot detect bi-directional relationships. In the S1 Text, we show the proportions of the IVs in the three conclusion groups respectively, called Steiger-Prop.

For any MR method, we first estimate and SE() for X → Y, and get a p-value p; for a given significance cutoff α, if p < α we conclude X has a causal effect on Y, otherwise we do not. Similar we make a conclusion on whether there is a causal relationship of Y → X.

For any CD method except Steiger’s, for direction X → Y, we estimate and SE(), for X → Y, then for a given significant cutoff α then we construct a (1 − α) confidence interval of K as , here z is the upper α/2 quantile of the standard normal distribution. If CI is completely within (-1,0) or (0,1), we conclude that X has a causal effect on Y, otherwise X does not. Similarly, we make a conclusion about Y → X.

4.6 Summary of different methods

In summary, we applied 15 different methods for inferring bi-directional causal relationships. They are in two groups denoted as CD and MR respectively. In the first group, Steiger-Prop and Steiger-MV are based on Steiger’s Method [6], CD-Ratio and CD-Egger are from [7], CD-Ratio-S and CD-Egger-S are the two methods with IV screening, CD-cML and CD-cML-DP are our newly proposed methods described in Section 4.2, so are CD-cML-S and CD-cML-DP-S with IV screening. In the second group, MR-cML (which does not use data-perturbation) and MR-cML-DP (which uses data-perturbation) are from [17], MR-cML-S and MR-cML-DP-S are the two methods with IV screening introduced here, and LHC-MR is from [24].

4.7 Main simulation setups

We independently generated 15 IVs in {g} with α1, ⋯, α15 from a uniform distribution Unif((−0.3, −0.2) ∪ (0.2, 0.3)); 10 IVs in {g} with β1, ⋯, β10 from Unif((−0.3, −0.2) ∪ (0.2, 0.3)); and 10 IVs in {g} with γ1, ⋯, γ10 from Unif((−0.3, −0.2) ∪ (0.2, 0.3)), and η1, ⋯, η10 from Unif((−0.3, −0.2) ∪ (0.2, 0.3)). We generated ξ’s in two ways: i) set ξ’s = 0 for no correlations with the confounder; ii) generated ξ’s from Unif(−0.1, 0.1) or Unif(−0.2, 0.2) for correlations with the confounder.

The MAFs of the SNPs/IVs were all set as 0.3. We generated continuous traits X and Y following the true causal model in Fig 1; ϵ, ϵ were independently drawn from N(0, 1), and ϵ from N(0,2). We also studied the scenarios with at least one of X and Y being binary. To generate binary traits, we generated the continuous X and Y first, then dichotomized one or both of them by setting the largest 30% of their values to be 1 and the other 70% as 0. We tried different combinations of (θ, θ) ∈ {0, 0.02, 0.1, 0.2, 0.3 } × {0, 0.02, 0.1, 0.2, 0.3}.

We generated two independent samples each of size n = N1 = N2 = 50000 from the reduced form (5) of the causal models for the two traits. Then we calculated and used subsequently the summary statistics for the two traits X and Y respectively. We set the significance cutoff at 0.05/35 to select relevant SNPs/IVs for both directions, and applied all methods for comparison.

4.8 A justification for binary traits

As shown in our simulations, both MR-cML and CD-cML performed well for binary traits, which was not coincident. This might seem surprising because our derivations for CD-cML are based on linear models. Although linear models are applicable to binary traits, often one would like to apply logistic regression as done in our simulations. Here we offer some explanations on why CD-cML, and more generally MR, work for binary traits in a general framework of bi-directional causal relationships in the presence of invalid IVs, which, to our knowledge, appears to be new. A key assumption is that a binary trait is obtained by dichotomizing a latent Normal (liability) trait. In addition, we assume that each SNP’s genetic effect is small, which is reasonable for complex traits.

Now consider a latent quantitative trait X as specified in Eq (5) with a shortened notation: where the error term ϵ (combining ϵ, ϵ and ϵ) is independent of SNP g; we further assume . A binary trait X* is defined as: X* = 1 if X > C for some constant cut-off C, and X* = 0 otherwise. Hence, conditional on the SNP g, we have where Φ(.) is the cumulative distribution function of the standard Normal N(0, 1); the first approximation is well known as given in [18] with H(.) = expit(.) being the inverse logit function, and the second approximation is based on a Taylor expansion of μ at , the value of μ when all the genetic effects are 0. Accordingly, we can express X* as a linear model sharing the same functional form as X in Eq (5) (except with an intercept term μ0 and a scaling factor C being newly added). Therefore, we can derive ρ and K’s as before.

The above formulation also explains why it is fine to fit either a linear model or a logistic regression model to estimate the marginal association between SNP g and the binary trait X*. To be concrete, let’s consider a valid IV g ∈ {g} for the causal direction of X to Y:

By Eqs (13) and (14), we see that we would obtain the same estimate r and its SE no matter whether we use or , obtained from the marginal logistic regression or linear regression respectively.

Although MR has been widely applied to both quantitative traits and binary traits, the usual justification for MR is based on linear models for quantitative traits; while there are some discussions on the use of a binary exposure or a binary outcome [18, 46], we are not aware of any systematic treatment of both a binary exposure and/or a binary outcome in the general framework of bi-directional relationships with possibly invalid IVs. Our above formulation offers a justification for the use of MR to binary traits in the general context. Again consider the causal direction from binary X* to continuous Y: with a valid IV g ∈ {g}, we have then the causal parameter is defined either as θ = β/[1.7/σ ⋅ β] or θ = β/[C ⋅ β], depending on whether a logistic regression model or a linear model is used for X*. In fact, θ may be also defined based on the association parameters between g and the latent X (and between g and Y). This is in agreement with [46]: if there is a discrepancy between the regression model being used in the definition of the causal parameter and the actual use of the model to analyze the binary exposure X*, the resulting MR estimate will be biased; nevertheless, there will be no discrepancy if there is no causal effect of X* on Y: we will always have θ = 0, regardless of its specific definition, if β = 0, thus it is always consistent to test for the presence of a causal relationship using any of the definitions or regression models. We can have a similar treatment for a binary Y*.

Supplementary file with additional simulation results, detailed real data analysis results, and more theoretical results and proofs.

(PDF)

Click here for additional data file.