Simultaneous estimation of bi-directional causal effects and heritable confounding from GWAS summary statistics

2.1 The underlying structural equation model

Let X and Y denote continuous random variables representing two complex traits. Let us assume (for simplicity) that there is one heritable confounder U of these traits. To simplify notation we assume that E(X) = E(Y) = E(U) = 0 and V ar(X) = V ar(Y) = V ar(U) = 1. The genome-wide sequence data for M sequence variants is denoted by G = (G₁, G₂,…, G_M). The aim of our work is to dissect the effects of the heritable confounding factor U from the bi-directional causal effects of these two traits (X and Y). For this we consider a model (see Figure 1) defined by the following equations:

• Download figure
• Open in new tab

Figure 1: Schematic representation of the extended structural equation model (SEM).

X and Y are two complex traits under scrutiny with a latent (heritable) confounder U with causal effects q_x and q_y on them. G represents a genetic instrument, with effects γ_x, γ_y and γ_u, respectively. Traits X and Y have causal effects on each other, which are denoted by α_x→y and α_y→x.

where γ_x, γ_y, γ_u ∈ ℛ^M denote the (true multivariable) direct effect of all M genetic variants on X, Y and U, respectively. All error terms (e_x, e_y and e_u) are assumed to be independent of each other and normally distributed with variances and , respectively.

Note that we do not include in the model reverse causal effects on the confounder (X → U and Y → U). The reason for this is the following: Let s_x and s_y denote those causal effect of X and Y on U. We can see that by reparameterising the original model to and , the genetic effects produced by the extended model with reverse causal effects on U and the simpler model (Figure 1) with the updated parameters are indistinguishable. Thus these extra parameters are not identifiable and the reparameterisation means that α_x→y and α_y→x in our model represent the total causal effects, some of which may be mediated by U.

Note that the model cannot be represented by classical directed acyclic graphs, as the bidirectional causal effects could form a cycle. However, the equations can be reorganised to avoid recursive formulation as follows: Regrouping the terms gives Substituting U into the first two equations yields with where and . Note that i_x is equivalent to the LD score regression intercept^[12].

We model the genetic architecture of these direct effects with a spike-and-slab distribution, assuming that only 0 ≤ π_x, π_y, π_u ≤ 1 proportion of the genome have a direct effect on X, Y, U, respectively and these direct effects come from a Gaussian distribution. Namely, Here, ⊙ denotes element-wise multiplication and ℬ_m(1, q) the m dimensional independent Bernoulli distribution. Further, we assume that all κ_x, κ_y, κ_us are independent of each other and so are all ζ_x, ζ_y, ζ_us. We can refer to as the direct heritability of X, i.e. independent of the genetic basis of U and Y. Similar notation is adapted for and . Note that when q_x = 0 and q_y ≠ 0 (or vice versa), this means that there is no confounder U present, but the genetic architecture of Y (or X) can be better described by a three component Gaussian mixture distribution.

We assume that the correlation (across markers) between the direct effects of a genetic variant on X, Y and U is zero, i.e. cov(γ_x, γ_y) = cov(γ_x, γ_u) = cov(γ_u, γ_y) = 0. Note that this assumption still allows for a potential correlation between the total effect of G on X and its horizontal pleiotropic effect on Y, but only due to the confounder U and through the reverse causal effect Y → X. As we argued above, this is a reasonable assumption, since the most plausible reason (apart from outcome-dependent sampling, which is out of the scope of this paper) for the violation of the InSIDE assumption may be one or more heritable confounder(s).

For simplicity, we also assume that the set of genetic variants with direct effects on each trait overlap only randomly, i.e. the fraction of the genome directly associated with both X and Y is π_x·π_y, etc. This assumption is in line with recent observation that the bulk of observed pleiotropy can be explained by extreme polygenicity with random overlap between trait loci^[13]. Note that uncorrelated effects (e.g. cov(γ_x, γ_y) = 0) do not ensure that the active variant sets overlap randomly, this is a slightly stronger assumption.

2.2 The observed association summary statistics

Let us now assume that we observe univariable association summary statistics for these two traits from two (potentially overlapping) finite samples N_x and N_y of size n_x, n_y, respectively. In the following, we will derive observed summary statistics in sample N_x and then we will repeat the analogous exercise for sample N_y. Let the realisations of X, Y and U be denoted by x, y and . The genome-wide genetic data is represented by and the genetic data for a single nucleotide polymorphism (SNP) k tested for association is . Note the distinction between the k-th column of G_x, which is the k-th sequence variant, in contrast to g_k, which is the k-th SNP tested for association in the GWAS. We assume that all SNP genotypes have been standardised to have zero mean and unit variance. The marginal effect size estimate for SNP k of trait X can then be written as , which is a special case of univariable standard normal linear regression when both the outcome and the predictor is standardised to have zero mean and unit variance^[12]. Note that x′ denotes the transposed of the column vector x. This can be further transformed as By denoting the linkage disequilibrium (LD) between variant k and all markers in the genome with we get with . Given the above-defined genetic effect size distribution the equation becomes Similarly, assuming that the LD structures (ρ_k) in the two samples are comparable, for estimated in the other sample (N_y) we obtain with .

Therefore, the joint effect size estimates can be written as Following the same rational as the cross-trait LD score regression^[14], the noise term distribution is readily obtained where r_x,y is the observational correlation between variables X and Y and n_x∩y is the size of the overlapping samples for X and Y. Since both n_x∩y and r_x,y cannot be estimated, we simply denote as the only estimated parameter and parameterise the covariance term as . Note that i_x,y is the cross-trait LD score regression intercept.

The bivariate probability density function (PDF) of these summary statistics cannot be obtained analytically, but in the following we demonstrate that the characteristic function can be derived. Let us first compute the characteristic function of this two-dimensional random variable, knowing that and are independent, hence the characteristic function can be factorised: In the following we will work out each of the characteristic functions on the right hand side.

2.3 Characteristic functions of and

It is reasonable to assume that linkage disequilibrium (LD) fades off beyond 1Mb distance. Thus, without loss of generality we can assume that non-zero LD does not extend beyond m₀ markers around the focal variant. Hence we can assume that the length of ρ_k is m₀ and only consider γ_x, γ_y and γ_u to be of length m₀ instead of m. Let us first approximate the distribution of ρ_k values following a spike and slab Gaussian mixture, i.e. proportion π_k of the m₀ SNPs have non-zero LD, coming from a Gaussian distribution and the remaining (1 − π_k) fraction of the LD values is zero. In mathematical notation Therefore can be written of the form The PDF of the product of two zero-mean Gaussians (r_k and ζ_u) is a modified Bessel function of the second kind of order zero (K₀(ω)) ^[15], more precisely and its characteristic function ^{[16, 17]} is Next, the characteristic function of the product of (r_k ⊙ ζ_u)_j and a Bernoulli distributed (κ_k,u)_j is Hence the characteristic function of the sum of m₀ independent random variables is the product of them, we have Finally, we apply a first order Taylor series approximation (around 1) of the log of the characteristic function in order to speed up computation and improve numerical accuracy Analogously, the approximation of the logarithm of the characteristic functions of and is Since the characteristic function of a centred multivariate Gaussian with variance-covariance matrix Σ is exp(−(1/2) · t′· Σ · t) we have

2.4 From characteristic function to probability density function

The final form of the logarithm of the joint characteristic function of the transformed summary statistics is Using the inversion theorem for characteristic functions we can express the joint distribution of as This integral can be efficiently computed by Fast Fourier Transformation (FFT, see ^[18] and references within). To speed up computation, we bin SNPs according to their π_k and σ_k values (10 × 10 bins with equidistant centres) and for SNPs in the same bin the PDF function is evaluated over a fine grid (2⁷ × 2⁷ combinations) using the FFT.

To reduce the number of parameters we define t_x := σ_u · q_x and t_y := σ_u · q_y since σ_u and q_x are separately not identifiable, but only their product is. Similarly π_u is unidentifiable, and is set to an arbitrary value of 0.1. For improved interpretability, we slightly reparameterise the likelihood function by using . Since different SNPs are correlated we have to estimate the over-counting of each SNP. We choose the same strategy as LD score regression^[12] and weigh each SNP by the inverse of its restricted LD score, i.e. , where r_jk is the correlation between GWAS SNPs k and j. The log-likelihood function is, thus, of the form where is the log-likelihood function value for SNP k. Parameters {n_x, n_y, m, σ_k=1,…,K, π_k=1,…,K} are known and the other 11 parameters are to be estimated from the observed association summary statistics. In order to speed up computation, we can estimate the 11 parameters in two separate steps: the first estimates for each trait the parameters π_x, i_x and π_y, i_y (SNP polygenicity and LD-score intercept) and the total heritability (unlike the direct heritability obtained by the full-model of LHC-MR) by using a simplified model with only the trait of interest, without a second trait or confounder, e.g. we fit only π_x, and i_x using and assume that π_x and i_x do not change when two traits are taken into account. Note that π_x may change slightly (decreasing from the total-to direct polygenicity), but its value has little impact on the likelihood function. The estimates from the first step can then be fixed for the parameter estimation of trait pairs. Since only π_x, i_x and π_y, i_y are fixed, the remaining parameters to estimate are now: It is key to note that our approach does not aim to estimate individual (direct or indirect) SNP effects, as these are handled as random effects. By replacing U with −U we swap the signs of both t_x and t_y, therefore these parameters are unique only if the sign of one of them is fixed. Thus, we will have the following restrictions on the parameter ranges: , t_x are in [0, 1], t_y, α_x→y, α_y→x, i_x,y are in [−1, 1].

2.5 Likelihood maximisation and standard error calculation

Our method, termed Latent Heritable Confounder Mendelian Randomisation (LHC-MR), maximises this likelihood function to obtain the maximum likelihood estimate (MLE). Due to the complexity of the likelihood surface, we initialise the maximisation using 50 different starting points, where they come from a uniform distribution within the parameter-specific ranges mentioned above. We then choose parameter estimates corresponding to the highest likelihood of the 50 runs. Run time depends on the number of iterations during the maximisation procedure, and is linear with respect to the number of SNPs. It takes ∼ 0.25 CPU-minute to fit the complete model to 50,000 SNPs with a single starting point.

Given the particular nature of the underlying directed graph, two different sets of parameters lead to an identical fit of the data, resulting in two global optima. The reason for this is the difficulty in distinguishing the ratio of the confounder effects (t_y/t_x) from the causal effect (α_x→y), as illustrated in Figure S2 by the slopes belonging to different SNP-clusters. More rigorously, it can be show that if {h_x, h_y, α_x→y, α_y→x, t_x, t_y} is an optimum, then so will be , where with w = t_y/t_x (for further derivations, see Supplementary Section 1.1). This allows us to directly obtain both optima, even if the optimisation only revealed one of them. It happens very often that one of these parameter sets are outside of the allowed ranges and hence can be automatically excluded. If not, we keep track of both parameter estimates maximising the likelihood function. Note that, we call the one for which the direct heritability is larger than the indirect one, i.e. , the primary solution. We show that for real data application this solution is far more plausible than the alternative optimum. Finally, note that such bimodality can be observed at different levels: (i) For one given data generation, using multiple starting points leads to different optima; (ii) LHC-MR applied to multiple different data generations for a fixed parameter setting can yield different optima. Both of these situations are signs of the same underlying phenomenon and most often co-occur.

We implemented the block jackknife procedure that is also used by LD score regression to calculate the standard errors. For this we split the genome into 200 jackknife blocks and compute MLE in a leave-one-block-out fashion yielding estimates. The variance of the full SNP MLE is then defined as .

2.6 Decomposition of genetic correlation

Given the starting equations for X and Y we can calculate their genetic correlation. Denoting the total (multivariate) genetic effect for X and Y as δ_x and δ_y, we can express them as follows Substituting the second equation to the first yields Similarly, Thus the genetic covariance is and the heritabilities are Therefore the genetic correlation takes the form These values can be compared to those obtained by LD score regression.

2.7 Computation of the LD scores

We first took 4,773,627 SNPs with info (imputation certainty measure) ≥0.99 present in the association summary files from the second round of GWAS by the Neale lab^[19]. This set was restricted to 4,650,107 common, high-quality SNPs, defined as being present in both UK10K and UK Biobank, having MAF > 1% in both data sets, non-significant (P_diff > 0.05) allele frequency difference between UK Biobank and UK10K and residing outside the HLA region (chr6:28.5-33.5Mb). For these SNPs, LD scores and regression weights were computed based on 3,781 individuals from the UK10K study^[20]. To estimate the local LD distribution for each SNP (k), characterised by π_k, , we fitted a two-component Gaussian mixture distribution to the observed local correlations (focal SNP +/− 2’500 markers with MAF≥ 0.5% in the UK10K): (1) one Gaussian component corresponding to zero correlations, reflecting only measurement noise (whose variance is proportional to the inverse of the reference panel size) and (2) a second component with zero mean and a larger variance than the first component (encompassing measurement noise plus non-zero LD).

2.8 Simulation settings

First, we tested LHC-MR using realistic parameter settings with a mild violation of the classical MR assumptions. These standard parameter settings consisted of simulating m = 234,000 SNPs for two non-overlapping cohorts of equal size (for simplicity) of n_x = n_y = 50,000 for each trait. X, Y and U were simulated with moderate polygenicity (π_x = 5×10⁻³, π_y = 1×10⁻², π_u = 5×10⁻²), and considerable direct heritability . U had a confounding effect on the two traits as such, q_x = 0.3, q_y = 0.2 (resulting in t_x = 0.16, t_y = 0.11), and X had a direct causal effect on Y (α_x→y = 0.3), while the reverse causal effect from Y to X was set to null. Note that in this setting the total heritability of each of these traits is principally driven by direct effects and less than 10% of the total heritability is through a confounder and in case of Y less than an additional 8% of its total heritability is through X.

It is important to note that for each tested parameter setting, we generated 50 different data sets, and each data generation underwent a likelihood maximisation of Eq. 2 using 50 starting points, and produced estimated parameters corresponding to the highest likelihood (simplified schema in Figure S3).

In the following simulations, we changed various parameters of these standard settings to test the robustness of the method. We explored how increased sample size (n_x = n_y = 500, 000) or differences in sample sizes ((n_x, n_y) = (50, 000, 500, 000) and (n_x, n_y) = (500, 000, 50, 000)) influence causal effect estimates of LHC-MR and other MR methods. We also simulated data with no causal effect (or with no confounder) and then examined how LHC-MR estimates those parameters. Next, we varied our causal effects between the two traits by lowering α_x→y to 0.1, and in another setting by introducing a reverse causal effect (α_y→x = −0.1). In addition, we tried to create extremely unfavourable conditions for all MR analyses by varying the confounding effects. We did this in several ways: (i) increasing q_x and q_y (q_x = 0.75, q_y = 0.50), (ii) having a confounder with causal effects of opposite signs on X and Y (q_x = 0.3, q_y = −0.2). We also drastically increased the proportion of SNPs with non-zero effect on traits X, Y and U (π_x, π_y and π_u = 0.1, 0.15, 0.2 respectively). We also simulated data whereby the confounder has lower (π_u = 0.01) polygenicity than the two focal traits.

Finally, we explored various violations of the assumptions of our model (see Section 2). First, we introduced two confounders in the simulated data, once with causal effects on X and Y that were concordant in sign, and another with discordant effects , while still fitting the model with only one U. Second, we breached the assumption that the non-zero effects come from a Gaussian distribution. By design, the first three moments of the direct effects are fixed: they have zero mean, their variance is defined by the direct heritabilities and they must have zero skewness because the effect size distribution has to be symmetrical. Therefore, to violate the normality assumption, we varied the kurtosis (2, 3, 5, and 10) of the distribution drawn from the Pearson’s distribution family. Third, we tested the assumption of the direct effects on our traits coming from a two-component Gaussian mixture by introducing a third component and observing how the estimates were effected. In this simulation scenario we introduced a large effect third component for X while decreasing the polygenicity of U (π_x1 = 1 × 10⁻⁴, π_x2 = 1 × 10⁻², , π_u = 1 × 10⁻²).

2.9 Application to real summary statistics

Once we demonstrated favourable performance of our method on simulated data, we went on to apply LHC-MR to summary statics obtained from the UK Biobank and other meta-analytic studies (Table S1) in order to estimate pairwise bi-directional causal effect between 13 complex traits. The traits varied between conventional risk factors (such as low education, high body mass index (BMI), dislipidemia) and diseases (including diabetes and coronary artery disease among others). SNPs with imputation quality greater than 0.99, and minor allele frequency (MAF) greater than 0.5% were selected. Moreover, SNPs found within the human leukocyte antigen (HLA) region on chromosome 6 were removed due to the abundance of SNPs associated with autoimmune and infectious diseases as well as the complicated LD structure present in that region. For traits with total heritability below 2.5%, the outgoing causal effect estimates were ignored since instrumenting such barely heritable traits is questionable.

In order to perform LHC-MR between trait pairs, a set of overlapping SNPs was used as input for each pair. The effects of these overlapping SNPs were then aligned to the same effect allele in both traits. To decrease computation time further (while only minimally reducing power), we selected every 10th QC-filtered SNP as input for the analysis. We calculated regression weights using the UK10K panel, which may be sub-optimal for summary statistics not coming from the UK Biobank, but we have previously shown^[21] that estimating LD in a ten-times larger data set (UK10K) outweighs the benefit of using smaller, but possibly better-matched European panel (1000 Genomes^[22]).

We also ran IVW for each trait pair in both directions to estimate bi-directional causal effects as well as LD score regression to get the cross trait intercept term. We then added uniformly distributed (∼ U (−0.1, 0.1)) noise to these pre-estimated parameters to generate starting points for the second step of the likelihood optimisation. These closer-to-target starting points did not change the optimisation results, simply sped up the likelihood maximisation and increased the chances to converge to the same (primary) optimum. The LHC-MR procedure was run for each pair of traits 100 times, each using a different set of randomly generated starting points within the ranges of their respective parameters. For the optimisation of the likelihood function (Eq. 2), we used the R function ‘optim’ from the ‘stats’ R package^[23]. Once we fitted this complete model estimating 11 parameters in two steps , we then ran block jackknife to obtain the SE of the parameters estimated in the second step: }.

To support the existence of the confounders identified by LHC-MR, we used EpiGraphDB^{[24, 25]} to systematically identify those potential confounders. The database provided for each potential confounder of a causal relationship, a causal effect on trait X and Y (r1, and r3 in their notation), the sign of the ratio of which (sign(r₃/r₁)) was compared to the sign of the LHC-MR estimated t_y/t_x values representing the strength of the confounder acting on the two traits. We restricted our comparison to the sign only, since the r1, r3 values reported in EpiGraphDB are not necessarily on the same scale.

2.10 Comparison against conventional MR methods and CAUSE

We compared the causal parameter estimates of the LHC-MR method to those of five conventional MR approaches (MR-Egger, weighted median, IVW, mode MR, and weighted mode MR) using a Z-test^[26]. The ‘TwoSampleMR’ R package^[27] was used to get the causal estimates for all the pairwise traits as well as their standard errors from the above-mentioned MR methods. The same set of genome-wide SNPs that were used by LHC-MR, were used as input for the package. SNPs associated with the exposure were selected to various degrees (for simulation we selected SNPs over a range of thresholds: absolute P-value < 5 × 10⁻⁴ to < 5 × 10⁻⁸), and SNPs more strongly associated with the outcome than with the exposure (P-value < 0.05 in one-sided t-test) were removed. The default package settings for the clumping of SNPs (r² = 0.001) were used and the analysis was run with no further changes. We tested the agreement between the significance and direction of our estimates and that of standard MR methods, with the focus being on finding differences in statistical conclusions regarding causal effect sizes.

We compared our causal estimates from all our simulation settings to the causal estimates obtained by running MR-RAPS ^[11] also using the ‘TwoSampleMR’ R package, once by using the entire set of SNPs, and another by filtering for SNPs with a significance threshold of < 5×10⁻⁴. We also compared both our simulation as well as real data results against those of CAUSE^[10]. We first generated simulated data under the LHC model and used them as input to estimate the causal effect using CAUSE. We then generated simulated data using the CAUSE framework and inputted them to LHC-MR (as well as standard MR methods) to estimate the causal parameters. Lastly, we compared causal estimates obtained for the 78 trait pairs (156 bi-directional causal effects) from LHC-MR to those obtained when running CAUSE.

3.1 Overview of the method

We fitted an 11-parameter structural equation model (SEM) (Figure 1) to genome-wide summary statistics of two studied complex traits in order to estimate bi-directional causal effects between them (for details see Methods). Additional model parameters represent direct heritabilities for X and Y, confounder effects, cross-trait and individual trait LD score intercepts and the polygenicity for X and Y. All SNPs associated with the heritable confounder (U) are indirectly associated with X and Y with effects that are proportional (ratio t_y/t_x). SNPs that are directly associated with X (and not with U) are also associated with Y with proportional effects (ratio 1/α_x→y). Finally, SNPs that are directly Y-associated are also X-associated with a proportionality ratio of 1/α_y→x. These three groups of SNPs are illustrated on the β_x-vs-β_y scatter plot (Figure S2). In simple terms, the aim of our method is to identify the different clusters, estimate the slopes and distinguish which corresponds to the causal- and confounder effects. In this paper, we focus on the properties of the maximum likelihood estimates (and their variances) for the bi-directional causal effects arising from our SEM.

3.2 Simulation results

We started off with a realistic simulation setting of 234,000 SNPs on chromosome 10 (LD patterns used from the UK10K panel) and 50,000 samples for both traits. Traits X, Y and confounder U had average polygenicity (π_x = 5 × 10⁻³, π_y = 1 × 10⁻², π_u = 5 × 10⁻²), with substantial direct heritability for X and , mild confounding (t_x = 0.16, t_y = 0.11) and a causal effect between X and Y (α_x→y = 0.3, α_y→x = 0). Note that with these settings, SNPs associated with U would violate the InSIDE assumption but might still be used by conventional MR methods. Under this standard setting, there were no genome-wide significant SNPs for standard MR methods, and estimates derived using SNPs with a p-value < 5 × 10⁻⁶ showed a downward bias for all MR methods (Figure 2 panel a). MR-RAPS using filtered SNPs (p-value < 5 × 10⁻⁴) was similarly downward biased whereas MR-RAPS using the entire set of SNPs was upward biased with the least amount of variance compared to all methods including LHC-MR. LHC-MR in this scenario slightly over estimated the causal effect in comparison but had the smallest RMSE after MR-RAPS (0.13 vs 0.06, Supplementary Table S2).

• Download figure
• Open in new tab

Figure 2: Simulation results under various scenarios.

These Raincloud boxplots^[28] represent the distribution of parameter estimates from 50 different data generations under various conditions. For each generation, standard MR methods as well as our LHC-MR were used to estimate a causal effect. The true values of the parameters used in the data generations are represented by the blue dots/lines. a Estimation under standard settings (π_x = 5 × 10⁻³, π_y = 1 × 10⁻², π_u = 5 × 10⁻², , t_x = 0.16, t_y = 0.11). b Addition of a reverse causal effect α_y→x = −0.2. c Confounder with opposite causal effects on X and Y (t_x = 0.16, t_y = −0.11).

We ran all our simulation scenarios with a smaller and a larger sample size (50,000 and 500,000) and observed that the relative performance of the methods were in some cases sample size specific. Smaller sample sizes often meant that standard MR methods had little to no IVs reaching genome-wide (GW) significance and hence we were forced to use IVs from less stringent thresholds (< 5 × 10⁻⁴ and < 5 × 10⁻⁶). Therefore, the causal effects were estimated with a substantial downward bias due to weak instrument bias (and winner’s curse). LHC-MR in these cases was able to estimate the causal effect with less bias but with a larger variance compared to most standard MR methods – still outperforming them in terms of RMSE in most settings. In the larger sample size setting, standard MR methods had IVs for every threshold cutoff. However, a pattern also observed with smaller sample sizes, but to a lesser extent, the causal estimates of some methods changed (either in mean or in variance, most noticeably observed in weighted median and IVW) as the threshold became more stringent. This is of particular concern and highlights that while in this simulation setting the 5 × 10⁻⁸ threshold may have optimally cancelled out the different biases for IVW (downward bias due to winner’s curse and weak instrument bias, upward bias due to genetic confounding), its estimate remains strongly setting-dependent. LHC-MR has performed reasonably well, exhibiting lower RMSE than most other methods, except for IVW and MR-RAPS for the 5 × 10⁻⁴ threshold (Figure S4 panel a). However, we observed that the performance of MR-RAPs is particularly setting- and threshold dependent.

Furthermore, unequal sample sizes for the two traits showed an underestimation of the causal effects for almost all MR methods, while LHC-MR remained the most accurate in the case where n_x (50,000) was smaller than n_y (500,000). However, the performances in the reverse scenario, where n_x was larger in size, were akin to the large sample size standard setting, where only IVW and filtered MR-RAPS (< 5 × 10⁻⁴) showed superior performance to LHC-MR both in terms of bias and variance (see Figure S5).

When testing scenarios in the absence of a causal- or a confounder effect (imitating the classical MR assumptions), with a smaller causal effect (α_x→y = 0.1), or with both forward- and reverse causal effects, we note that LHC-MR outperforms the standard MR methods as well as MR-RAPS in all these scenarios.

When there was no causal effect (α_x→y = 0), LHC-MR had the smallest bias out of all the methods in both sample sizes (0.004 in both, Figure S6 panel a and Figure S7 panel a). The variance of the LHC-MR estimates in the larger sample size was much lower (0.0001 vs 0.01), similarly the other methods had a smaller variance in larger sample sizes and had more clearly seen upward biased estimates. The increased upward bias of standard MR methods is due to the fact that confounder-associated SNPs could only be detected in larger sample size and those lead to positive bias (due to the concordant effect of the confounder on the two traits). Note that the variance of standard MR methods are low simply because, in these settings, we were forced to lower the instrument selection threshold, hence artificially included many (potentially invalid) instruments, which lowers the estimator variance while increasing bias. MR-RAPS greatly overestimates the causal effects when the sample size is larger.

In the absence of a confounder effect, there is not much of a difference between the two sample sizes; standard MR methods have a large variance and are downward biased, LHC-MR is less biased compared to them but MR-RAPS performs best with the least bias and variance when all the SNPs are used as instruments (Figure S6 panel b and Figure S7 panel b). Trying a smaller causal effect led to an upward bias for all MR methods including both filterings of MR-RAPS in the larger sample size. Alternately, when n_x = n_y = 50, 000, the MR methods are downward biased (Figures S6 panel c and Figures S7 panel c). Lastly, when a (negative) reverse causal effect is introduced, all MR methods and MR-RAPS are negatively biased in their estimation of the causal effect (see Figure 2 panel b). LHC-MR has a much smaller bias for the forward causal effect estimate in this case, and a generally small bias for the reverse causal effect in both sample sizes (0.05 for n = 50, 000 and 0.03 for n = 500, 000, Figure S4 panel b).

Increasing the indirect genetic effects, by intensifying the contribution of the confounder to X and Y (t_x = 0.41, t_y = 0.27), led to a general over estimation of the causal effects by all methods including LHC-MR, though more drastically seen in standard MR methods and MR-RAPS in larger sample sizes, when there is sufficient power to pick up these confounder-associated SNPs. The causal effect estimates of standard MR methods in the smaller sample size were much less affected by the presence of a strong confounder compared to LHC-MR and MR-RAPS (Figure S8). The reason for this is that the confounder-associated SNPs remain undetectable at lower sample size and hence instruments will not violate the classical MR assumptions.

Further testing the effects of the confounder trait on the causal estimation, we tested the impact of confounders with opposite effects on X and Y. We observe a major underestimation of the causal effects for standard MR methods as well as MR-RAPS, whereas LHC-MR performs better for both sample sizes (RMSE = 0.01 and 0.1 for larger and smaller n respectively), see Figures 2 panel c and S4 panel c.

Our LHC-MR method is influenced by the unlikely scenario of extreme polygenicity for traits X, Y and U, and it suffers from increased bias and variance regardless of sample size (see Figure S9). Standard MR methods as well as filtered MR-RAPS underestimated the causal effect when n = 50, 000. Some also underestimated α_x→y when n = 500, 000, with the exception of IVW, Mode and filtered MR-RAPS, that outperformed the rest. Decreasing the proportion of confounder-associated SNPs to 1% only, does not seem to affect our method and shows similar results to the standard setting (Figure S10).

Furthermore, we simulated summary statistics, where (contrary to our modelling assumptions) the X−Y relationship has two confounders, U₁ and U₂. When the ratio of the causal effects of these two confounders on X and Y ( and respectively) agreed in sign, the corresponding causal effects of standard MR methods were over-estimated in larger sample sizes and, conversely, underestimated in smaller sample sizes (Figures S11 and S12, panels a). LHC-MR and weighted median performed better however in larger sample sizes and had a bias of 0.03 and 0.07 respectively. However, when the signs were opposite (, for U₁ and , for U₂), conventional MR methods and MR-RAPS in this case almost all underestimated the causal effect regardless of sample size. LHC-MR outperformed them both in the larger sample size (bias of 0.007) and in the smaller sample size (bias of −0.003), see Figures S11 and S12, panels b.

Finally, we explored how sensitive our method is to different violations of our modelling assumptions. First, we simulated summary statistics when the underlying non-zero effects come from a non-Gaussian distribution. Interestingly, we observed that, for smaller sample sizes, the variance of the causal effect estimate was dependent on the kurtosis for most MR methods. LHC-MR estimations yielded slightly more pronounced upward bias than IVW, while still exhibiting the lowest RMSE among all methods (Figure 3 panel a). Similar results are seen in larger sample size with smaller variance for all methods under all degrees of kurtosis except for IVW, which showed a better performance than LHC-MR (Figure S13 panel a). Second, we simulated effect sizes coming from a three-component Gaussian mixture distribution (null/small/large effects), instead of the classical spike-and-slab assumption of our model. The smaller sample size estimates mirror those of the standard setting with n also equal to 50, 000 (see Figure 3 panel b). However, in the larger sample size, LHC-MR overestimates the causal effect. This bias could be due to the merging of true effect estimates with confounder effect leading to an overestimation of α_x→y (Figure S13 panel b). MR-Egger, IVW and filtered MR-RAPS have the smallest RMSE in this case.

• Download figure
• Open in new tab

Figure 3: Simulation results under various scenarios.

These Raincloud boxplots^[28] represent the distribution of parameter estimates from 50 different data generations under various conditions. For each generation, standard MR methods as well as our LHC-MR were used to estimate a causal effect. The true values of the parameters used in the data generations are represented by the blue dots/lines. a The different coloured boxplots represent the underlying non-normal distribution used in the simulation of the three γ_x, γ_x, γ_u vectors associated to their respective traits. The Pearson distributions had the same zero mean and skewness, however their kurtosis ranged between 2 and 10, including the kurtosis of 3, which corresponds to a normal distribution assumed by our model. The standard MR results reported had IVs selected with a p-value threshold of 5×10⁻⁶. b Addition of a third component for exposure X, while decreasing the strength of U. True parameter values are in colour, blue and red for each component (π_x1 = 1 × 10⁻⁴, π_x2 = 1 × 10⁻², ).

3.3 Application to association summary statistics of complex traits

We applied our LHC-MR and other MR methods to estimate all pairwise causal effects between 13 complex traits (156 causal relationships in both directions). Our results are presented as a heatmap in Figure 4 (and are detailed in Supplementary Table S3). Further, we calculated the alternate set of estimated parameters that naturally results from our model (for reference see Sections 2.2 and Supplementary materials 1.1). Among trait pairs for which the exposure had sufficient heritability (> 2.5%), the alternate parameters of a 102 trait pairs were within the possible ranges mentioned in methods (i.e. the confounder and the exposure are interchangeable). However, for all of these pairs, the alternative parameter optima lead to lower direct-than indirect heritability, which we deem unrealistic. Therefore, we report only the primary set of estimated optimal parameters in the main results and provide the alternative parameters in the Supplementary Table S4. The comparison of the results obtained by LHC-MR and standard MR methods is detailed below and more extensively in Supplementary Tables S5-S6. In summary, LHC-MR provided reliable causal effect estimates for 132 out of 156 exposure traits (i.e. those exposures had an estimated total heritability greater than 2.5%). These estimates were compared to five different MR methods. Seventy-four causal relationships were deemed significant by LHC-MR. Furthermore, for 117 out of those 132 comparable causal relationships, our LHC-MR causal effect estimates were concordant (not significantly different) with at least two out of five standard MR methods’ estimates.

• Download figure
• Open in new tab

Figure 4: Heatmap representing the bi-directional causal relationship between the 13 UK Biobank traits.

The causal effect estimates in coloured tiles all have a significant p-values surviving Bonferroni multiple testing correction with a threshold of 3.2×10⁻⁴. We did not report an estimated causal effects for exposures with an estimated total heritability less than 2.5%. White tiles show an absence of a significant causal effect estimate. Abbreviations: BMI: Body Mass Index, BWeight: Birth Weight, CAD: coronary artery disease, DM: Diabetes Mellitus, Edu: Years of Education, HDL: High-density Lipoprotein LDL: Low-density Lipoprotein, MI: Myocardial Infarction, PSmoke: # of Cigarettes Previously Smoked, SBP: Systolic Blood Pressure, SHeight: Standing Height, SVstat: Medication Simvastatin

By simply comparing the significance status and the direction of the causal effects between the methods, we see that LHC-MR agrees in sign and significance (or the lack there of) with at least 3 MR methods 77 times. For 31 relationships, LHC-MR results lead to different conclusions than those of standard MR methods. For 28 of those, LHC-MR identified a causal effect missed by all standard MR methods. For the other three, we observed a disagreement in sign: LDL has a negative effect on BMI according to weighted mode and weighted median, whereas we show a positive effect, HDL and LDL show a negative bi-directional causal effect for weighted mode but a positive bi-directional effect with LHC-MR. Despite the conflicting evidence for the causal relationship of LDL on BMI, studies have shown that the relationship between them is non-linear^[29], possibly explaining the discrepancy between the results.

LHC-MR agreed with most MR estimates and confirmed many previous findings, such as increased BMI leading to elevated blood pressure^{[30, 31]}, diabetes mellitus^{[32, 33]} (DM), myocardial infarction^[34] (MI) and coronary artery disease^[35] (CAD). Furthermore, we confirmed previous results^[36] that diabetes increases .

Interestingly, it revealed that higher BMI increases smoking intensity, concordant with other studies^{[37, 38]}. It has also shown the protective effect of education against a range of diseases (e.g. CAD and diabetes^{[39, 40]}) and risk factors such as smoking^{[41, 42]}, in agreement with previous observational and MR studies. Probably reflecting lifestyle change recommendations by medical doctors upon disease diagnosis, statin use is greatly increased when being diagnosed with CAD, (systolic) hypertension, dislipidemia, and diabetes as is shown by both LHC-MR and standard MR methods.

Furthermore, causal effects of height on CAD, DM and SBP have been previously examined in large MR studies^{[43, 44]}. LHC-MR, agreeing with these claims, did not find significant evidence to support the effect of height on DM, but did find a significant protective effect on CAD and SBP. However, unlike the first two, the relationship between height and SBP also revealed the existence of a confounder with causal effects 0.14 (P = 9.2 × 10⁻¹¹) and 0.11 (P = 3.39 × 10⁻⁸) on height and SBP respectively. Another example of a trait pair for which LHC-MR found an opposite sign confounder effect is HDL and its protective effect on SBP. The confounder had a positive effect ratio of t_y/t_x = 0.84, opposing the negative causal effect of supported by observational studies^[45]. This causal effect was not found by any other MR method.

It is important to note that while the effects of parental exposures on offspring outcomes can be seen as genetic confounding, LHC-MR would not be able to distinguish parental and offspring causal effects, because the LHC-MR model assumes that there is no correlation between the genetic effects on the exposure and the genetic effects on the confounder (which is not the case of parental vs offspring traits). Thus, LHC-MR causal effect estimates are just as likely to reflect parental effects as any other MR method ^[46]. This may be the case, for example, for the detrimental effect of increased (parental) BMI on education (supported by longitudinal studies^[47]), the positive effect of (parental) height on birth weight^[48], or on education^[49]. There are also some associations identified only by LHC-MR that might reflect parental effects: the negative causal effect of CAD on education or on birth weight, the positive impact of HDL on birth weight, or DM reducing height. All these pair associations uniquely found by LHC-MR are examples of LHC-MR’s use of whole-genome SNPs instead of GW-significant SNPs only, as our estimates are of larger magnitude than those found by standard MR. Interestingly, for the CAD→birth weight relationship, LHC-MR revealed a confounder of opposite causal effects, which could have masked/mitigated the causal effect of standard MR methods.

A systematic comparison between IVW and LHC-MR has shown generally good agreement between the two methods, which is illustrated in Figure 5. To identify discrepancies between our causal estimates and those of the standard MR results, we grouped the estimates into several categories, either non-significant p-value for both or either, significant with an agreeing sign for the causal estimate, or significant with a disagreeing sign. The diagonal (seen in Figure 5) representing the agreement in significance status and sign between the two methods, is heavily populated. On the other hand, 34 pairs have causal links that are significantly non-zero according to LHC-MR, but are non-significant for IVW, while the opposite is true for seven pairs. We believe that many of these seven pairs may be false positives, since four of them are picked up by no other MR method, two are confirmed by only one other method and the last one by two methods. Further comparisons of significance between LHC-MR estimates and the remaining standard MR methods can be found in Table S7.

• Download figure
• Open in new tab

Figure 5: A scatter plot of the causal effect estimates between LHC-MR and IVW.

To improve visibility, non-significant estimates by both methods are placed at the origin, while significant estimates by both methods appear on the diagonal with 95% CI error bars.

LHC-MR identified a confounder for 16 trait pairs out of the possible 78. In order to support these findings, we used EpiGraphDB^{[24, 25]} to systematically identify those potential confounders. EpiGraphDB could identify reliable confounders for ten out of the 16 trait pairs. Notably, for the birth weight - diabetes pair, the average epigraph confounder-effect ratio (r₃/r₁) clearly agreed in sign with our t_y/t_x ratio, indicating that the characteristics of the confounder(s) evidenced by LHC-MR agree with those found in an exhaustive confounder search, and are mainly obesity-related traits (Figure S18 panel a). Six other trait pairs showed mixed signs of different confounders, indicating the possibility of having heterogeneous confounders (Figure S18 panels b-e). Finally, three trait pairs showed a disagreement between our estimated confounder effect ratio and the bulk of those found by epighraphDB as seen in Supplementary Figure S18 panels f-j. However, at least one of the top ten potential confounders showed effects that are in agreement with our ratio for each of these pairs. Note that since the reported causal effects of the confounders on X and Y reported in EpiGraphDB are not necessarily on the same scale, we do not expect the magnitudes to agree.

As described in the methods (Eq. 3), genetic correlation can be computed from our estimated model parameters. To verify that the fitted LHC-MR model leads to a genetic correlation similar to the one obtained from LD score regression^[50] (LDSC), we compared whether the two approaches produce similar genetic correlation estimates. We did this by taking the estimated parameters obtained from the 200 block jackknife to estimate the genetic correlations between traits (and their standard errors), and plotted them against LD score regression values as seen in Figure 6. As expected, we observe an overall good agreement between the estimates of the two methods, with only six trait pairs differing in sign. Of these six, only 2 were nominally significantly different between the two methods (LDL→Asthma and LDL→DM). Further de-composition of the genetic covariance into heritable confounder-led or causal effect-led revealed that most of the genetic covariance between traits can be attributed to bi-directional causal effects. A reason for this could be that confounders would need to have very strong effects to substantially contribute to the genetic correlation (≈ t_x · t_y) compared to the bi-directional causal effects .

• Download figure
• Open in new tab

Figure 6: Scatter plot comparing the genetic correlation for each trait obtained from LDSC against the value calculated using parameter estimates from the LHC-MR model.

A 95% CI is shown for each point. Values from both methods are reported in Supplementary Table S9.

As for the comparison of LHC-MR against CAUSE for real trait pairs, we ran CAUSE on all 156 trait pairs (bi-directional), and extracted the parameter estimates that corresponded to the methods winning model. The p-value threshold was corrected for multiple testing and was equivalent to 0.05/156. Based on that threshold, the p-value that compared between the causal and the sharing model of CAUSE was used to choose one of the two. Then the parameters estimated from the winning model, γ (only for causal model), η and q, were compared to their counterparts in LHC-MR. A visual comparison of LHC-MR’s causal estimates and those of CAUSE can be seen in Figure S19.

Whenever the causal effect estimates were significant both for CAUSE and LHC-MR (30 causal relationships), they always agreed in sign (Table S8) with a high Pearson correlation of 0.592. Calculating the correlation for their estimates regardless of significance yielded a smaller value of 0.377. When compared to the causal effect estimate from IVW, LHC-MR was strongly correlated (0.585), whereas CAUSE had a slightly weaker correlation (0.471) using all estimates.

Similarly, the significant confounder effect ratio of LHC-MR (t_y/t_x) can be compared to the significant confounder effect estimate of CAUSE (η) when a sharing model is chosen. These 12 confounding quantities by CAUSE and LHC-MR disagreed in sign for all but one trait pair (Height→MI), with a Pearson correlation compatible with zero (−0.357 (95% CI [-0.77, 0.27])).