Abstract
Mendelian randomization (MR) is an approach to causal inference that utilises genetic variants to obtain estimates of the causal effect of an exposure on an outcome in the presence of unobserved confounding. MR relies on a set of assumptions to obtain unbiased effect estimates, one of these assumptions is that there is no pathway from the genetic variants to the outcome that does not act through the exposure. Increasing genome-wide association study (GWAS) sample sizes for the exposure enables discovery of instrumental variables with smaller effect sizes. We illustrate through simulations how smaller effect sizes could arise from genetic variants that act through traits that have greater liability to confound an exposure-outcome relationship. When such genetic variants are selected as instruments this can bias the MR effect estimate obtained from that instrument in the same direction as the confounded observational association but with larger magnitude. Through simulation we illustrate how the total bias of the MR estimates increases across a range of standard MR estimation methods increases as the proportion of the genetic instruments that are associated with the confounder increases. However, if such heritable confounders are known and can be instrumented, the confounder free effect estimate can be obtained through applying a pre-estimation filtering to standard MR methods, removing instruments that explain more variation in that confounder than the exposure, or by estimating effects through multivariable MR. We highlight the potential for SNPs identified in GWAS to be associated with potential confounders through examination of a recent GWAS of C-Reactive Protein. Finally, we illustrate our approach through estimation of the causal effect of age at menarche on type 2 diabetes, hypothesising that the MR effect estimate may be biased by confounding due to the inclusion of genetic variants associated with early life adiposity as instruments.
Introduction
Confounding in epidemiological studies can bias observed estimated associations between traits so that they cannot be interpreted as the causal effect of the exposure trait on the outcome. Adjusting for confounders will remove this bias provided that all confounders are known and measured without error, which is rarely plausible.[1] Mendelian Randomisation (MR) is an alternative approach to causal effect estimation which utilises genetic variants to infer causality of the exposure on the outcome that is not biased by that confounding.[2, 3] MR is most often implemented by using genetic variants associated with the exposure within an instrumental variable (IV) estimation framework.[4]
IV estimation relies on three key assumptions to be able to test for a causal effect of the exposure on the outcome. These assumptions are (IV1) the instrument is associated with the exposure, (IV2) there is no confounding of the instrument and the outcome and (IV3) there is no effect of the instrument on the outcome that doesn’t act through the exposure. These assumptions are illustrated in Figure 1A. In MR genetic variants associated with the exposure of interest are used as IVs. These are usually single nucleotide polymorphisms (SNPs) identified through genome-wide association studies (GWAS).[5] GWAS estimate the association between millions of genetic variants and a trait of interest. Variants associated in a GWAS with the exposure trait at a predetermined level are then selected to be used as IVs in the MR estimation. The level chosen is typically set to select all independent variants where the p-value for the estimated association between the variant and the exposure is smaller than 5×10−8. This level is chosen as it is the Bonferroni corrected 5% significance level, given the effective number of independent tests being conducted in a GWAS.
The polygenic nature of complex traits, for which there is substantial empirical support[6], posits that the majority of heritable variation is driven by a very large number of genetic factors, each with very small effect sizes. To discover such variants GWAS has been growing in sample size rapidly in recent years, especially through the emergence of large scale biobanks. The most recent GWAS of body mass index (BMI) included approximately 700,000 individuals and identified 941 independent SNPs associated with BMI.[7] The most recent GWAS of educational attainment included approximately 3 million individuals and identified nearly 4000 independent SNPs.[8] Similarly, this growing availability of large sample sizes means that for some traits where previously no SNPs were identified as genome-wide significantly associated with the trait (i.e. no SNPs had a p-value < 5×10−8) large studies can now identify enough SNPs to conduct MR estimation.
For MR analyses involving exposures that are themselves complex traits, it is well understood that genetic variants associated with the exposure are mediated through other phenotypic layers.[9] Such mediation opens the possibility of horizontal pleiotropy if those mediating traits also have effects on the outcome through pathways independent of the exposure. If some, or all, of the variants associated with the exposure exhibit horizontal pleiotropy the MR effect estimates can be biased through violation of assumption IV3. There are a number of different forms that horizontal pleiotropy can take and so bias MR.[4] Correlated pleiotropy is used to describe pleiotropy that causes bias in MR which is correlated with the magnitude of association between the SNPs and the exposure. Correlated pleiotropy can occur through heritable confounding, misspecification of the primary phenotype or reverse causation, illustrated in Figure 1B. The distinction between heritable confounding and misspecification of the primary phenotype (i.e. the forms of pleiotropy illustrated by Figure 1B(ii) and 1B(iii)) is not always clear, especially as one may lead to the other. If the genetic variant for a heritable confounder is taken to be an instrument for the exposure that the heritable confounder influences then this misspecification of the primary phenotype will clearly lead to erroneous inference, with the extreme case being when the upstream phenotype is the outcome under investigation, as in Figure 1B(iv). As each of these mechanisms lead to the same bias in the resulting MR estimates and misspecification of the primary phenotype could be considered to be an extreme form of heritable confounding. Note however that for there to be bias in the MR estimation the path from the confounder to the outcome is required. If there is no effect of the confounder on the outcome that does not act through the exposure then MR estimates will not be biased. Throughout the rest of this paper we consider correlated pleiotropy driven by genetic variants associated with a confounder of the exposure and outcome being included within the set of instruments used in an MR analysis and refer to this as heritable confounding.
There remains a question of whether increasing GWAS sample sizes will lead to increased horizontal pleiotropy through heritable confounding, i.e. are smaller genetic effects for the exposure more likely to also have an independent path to the outcome for any given exposure-outcome relationship? One mechanism by which this could arise is if genetic variants that have effects on the exposure that are mediated through other traits tend to be those that are more phenotypically connected in the overall genotype-phenotype map. Such variants would be likely to have smaller effects on a trait and so only become discoverable as GWAS sample sizes increase. One implication of this is that the genetic variants with the largest effect sizes in the GWAS may be the least likely to be pleiotropic instruments.[10] Although the extent to which this will hold will depend on the genetic architecture of the trait. If all of the pathways from those upstream mediating traits to the outcome act through the exposure then the MR effect estimates will not be biased and the increased GWAS sample size will increase the power of the MR estimate to detect causal effects. However, when one (or more) of the upstream traits is a confounder of the exposure and outcome this will lead to violation of the IV3 assumption and biased MR effect estimates.
There are a number of MR methods available that are robust to violations of assumption IV3, each of which assume particular forms for the structure of the pleiotropy[4]. The most commonly used methods are MR Egger[11], weighted Median[12] and weighted Mode[13]. MR Egger assumes that the size of any pleiotropic effect of each SNP on the outcome is uncorrelated with the SNP-exposure association. Under correlated pleiotropy this assumption will not hold as the size of the pleiotropic effect will be associated with the size of the SNP-exposure association. Weighted Median assumes that more than half of the SNPs (weighted by their inverse variance weights) are not pleiotropic. Weighted Mode assumes that the modal effect estimate obtained by the weighted SNPs is not biased by pleiotropy and so assumes that a plurality of the SNPs are not pleiotropic. The assumptions made by these methods mean that the weighted mode and weighted median approaches will be robust to a small proportion of the SNPs included in the estimation being associated with a confounder of the exposure and outcome but will both be biased as the proportion of SNPs associated with confounders of the exposure and outcome increases. IVW and MR Egger will be biased if any of the SNPs are associated with the confounder, although the level of the bias will depend on the proportion of the SNPs used as instruments that are associated with the confounder, and the presence/absence of another confounder balancing the bias in the opposite direction. This result has been illustrated elsewhere.[14]
Alternative methods of MR estimation are available that adjust for correlated pleiotropy under particular assumptions of the structure the pleiotropy takes.[14-21] These approaches are often computationally intensive, using Bayesian estimation approaches applied to a large number of genetic variants. They also impose an assumption of a structure of the relationship between the pleiotropic variants and other variants, often assuming that the plurality of the instruments are valid.
Here we propose two alternative approaches to estimation in the presence of correlated pleiotropy when the confounder or confounders are hypothesised and measured; including the confounder in a multivariable MR (MVMR) and removing the SNPs associated with the confounder through confounder-based Steiger filtering. Each approach relies on knowledge of the likely confounders, however when this is unknown the other traits associated with the genetic instruments, and so most likely to be confounders, can be identified through PheWAS analysis.[22] The advantage of these approaches are that they are simple to implement and already widely used in MR analyses. They therefore provide a simple approach to explore the potential for bias from heritable confounding in MR studies.
The remainder of this paper is structured as follows; we first compare the bias from confounding in IV estimation to linear regression. We explore the question of whether increasing GWAS sample sizes could increase the rate of correlated pleiotropy through a simulation of a network of traits explained by set of genetic variants. We then present a simulation study to illustrate the bias induced by correlated pleiotropy and how MVMR and confounder-based Steiger filtering can correct the bias. We consider the example of C-reactive protein (CRP) and explore how many potential confounders are likely to be associated with SNPs identified in the most recent large GWAS of CRP. Finally, we consider estimation of the causal effect of age at menarche on risk of type 2 diabetes and show that in this analysis the univariable IVW MR estimate is likely to be biased by correlated pleiotropy from SNPs associated with childhood adiposity.
Results
Bias in IVW MR estimates compared to linear regression
In this section we consider a simple model of a relationship between an exposure (X) and an outcome (Y) which is confounded by a single unobserved confounder (U). The exposure is associated with a single instrument (Z) that also has an effect on the outcome via the unobserved confounder. This single instrument could be considered to be a genetic risk score constructed from all SNPs identified as strongly associated with the exposure in a GWAS. We additionally assume that all relationships are linear and there are no interactions. We focus on the setting with a single heritable confounder for simplicity and clarity however in any application there would also be a number of non-heritable confounders as well as potentially more than one heritable confounder.
The results given here therefore only compare the results for bias from a single heritable confounder and not the overall bias that would be expected. The model we consider can be written as; Where vy, vx and vu are all uncorrelated random error terms. This model is illustrated in Figure 2A. In this model the estimate of the effect of the exposure on the outcome from linear regression is given by;
As we have assumed no other confounders in this simple model vy is independent of the exposure X and so E(∑ Xvy)= 0. We can additionally define δux as the estimate that would be obtained from a regression of X on U. Note: this will not be the same as βux due to the confounding of U and X by Z but will estimate the correlation between X and U. Therefore; This can alternatively be expressed as; Where is the variation in X explained by U. i.e. the true causal effect of the exposure on the outcome plus a bias term that depends on the effect of the confounder on the outcome and the correlation between the confounder and the exposure.
The estimate of the effect of the exposure on the outcome from IV estimation is given by; The instrument Z is independent and exogenous so taking expectations of this expression, E(∑ Zvy )= 0, E ∑ Zvx) = 0 and E ∑ Zvu 1 = 0. Therefore; This is the true causal effect of the exposure on the outcome plus a bias term that depends on the ratio of the pleiotropic effect of the instrument on the outcome (βuy πzu) to the effect of the instrument on the exposure (πzx + πux πzu). i.e. the larger the pleiotropic effect as a proportion of the effect of the instrument on the exposure the larger the bias of the IV estimator. The bias in the MR estimate from a single heritable confounder will be less than the bias of the linear regression estimator from the same confounder when; The smaller the pleiotropic effect as a proportion of the effect of the instrument on the exposure the more likely this is to hold. Under a simpler model of heritable confounding, where U is a confounder of X and Y and has its own independent genetic profile in which Z influences only one path, i.e. πzx = 0 the bias term for the IV estimator reduces to . The implication of this is that if the genetic variants influence the exposure x only through the confounder U the bias in the MR estimate from that confounder will always be greater than the bias in the observational association from the same confounder.
To illustrate the relationship between the bias of linear regression and IV estimation across different levels of pleiotropy we simulated the model given in Figure 2A varying the effect of the instrument on the confounder but holding all other parameters constant. We estimated βxy using linear regression and the two-stage least squares IV estimator. The mean bias of each estimator across 1000 repetitions for varying levels of confounding is given in Figure 2B. These results illustrate how the bias of the IV estimator increases as the pleiotropic effect via the confounder increases and can be substantially more than the bias from linear regression. Note that the exact position and slope of these lines will depend on the other parameters in the model. Figures 2C and 2D show the same bias plots as Figure 2B with an additional non-heritable confounder added to the model. The direction of the confounding from this additional confounder is varied from positive (Figure 2C) to negative (Figure 2D). These plots highlight how the overall bias of each estimator will depend on other factors as well as the level of heritable confounding and the bias of the estimate from MR may not always act in the same direction as the bias of the estimate from linear regression even in the presence of heritable confounding due to other confounding.
Increasing GWAS sample size may increase correlated pleiotropy
GWAS sample sizes of complex traits are continually growing, pursuing a well-worn strategy to uncover hundreds or thousands of smaller effect sizes that in aggregate will explain a large polygenic component of phenotypic variation. A genetic variant inducing correlated pleiotropy must influence a trait that is a parent to both exposure and outcome. We hypothesised that such a genetic variant would generally be associated with more traits upstream of the exposure than valid instruments, and therefore would have relatively smaller effects on the exposure than valid instruments. The consequence of such a relationship would be that increasing GWAS sample sizes to identify more instruments of smaller effect sizes could lead to increasing bias in MR estimates.
To evaluate this hypothesis we performed simulations in which we consider that all hypothesised exposure-outcome relationships exist in a broader genotype-phenotype network. The network represents a directed acyclic graph where all nodes are traits except those that have no parents, which are genetic variants. The simulation proceeds by randomly choosing two traits, one as exposure and one as outcome, and then evaluating a) which variants are instruments for the exposure, and whether or not the instruments act through confounders in the graph, and b) the path length from instrument to exposure. The length of the edge between any pair of nodes has been fixed to be the same length for all edges.
Figure 3 illustrates that in this model instruments that act via confounders tend to have smaller effects due to more traits lying on the path from the genetic variant to the exposure. For traits where our simulation is a realistic model of the underlying biology an implication of this result is that as sample sizes grow for GWAS, the rate of discovery of genetic instruments that are liable to act via confounders will grow. However, it is possible that this will not hold for all traits, for some traits it may be that the strongest effects are via confounders (e.g. because those confounders have very strong effects on the exposure) and the genetic variants with smaller effects will be the unconfounded ones. It is therefore important to consider if this model is likely to hold for each exposure trait considered.
Simulating the impact of genetic confounding on MR estimates
We simulated a model with an exposure and outcome which are confounded by a third (potentially unobserved) heritable trait. Here exposure has no causal effect on the outcome however, linear regression of the exposure on the outcome would give a positive estimate of the association between the exposure and outcome due to confounding. The exposure is influenced by two traits, one simply mediates genetic effects onto the trait (i.e. influences on X), whereas the other is a heritable confounder of the exposure and the outcome. Both mediator and confounder are each associated with a set of genetic variants. This model is illustrated in Figure 4, which corresponds to the scenario illustrated in Figure 1B(ii). We have included only a single heritable confounder in these models, no non-genetically influenced confounders and have assumed that all other IV assumptions hold. The results therefore illustrate the bias resulting from that confounder and not the overall bias that would be expected in either the linear regression or the MR estimates.
The genetic variants associated with the mediator are valid instruments for estimating the effect of the exposure on the outcome whereas genetic variants associated with the confounder are not as there is a path from them to the outcome that does not act via the exposure. Therefore, MR using just the set of genetic variants associated with the mediator (G1) would estimate the true null causal effect of the exposure on the outcome. However, if the researcher mistakenly includes genetic variants in G2 as instruments for X, bias in the MR will be introduced through the pleiotropic pathway. We use this simulation to illustrate how MR is biased under correlated pleiotropy and explore the ability of commonly used pleiotropy robust methods to identify this bias when different proportions of the total set of genetic variants are in sets G1 and G2.
Results from this simulation are given in Figure 4 and Supplementary Tables 1 and 2. When the total proportion of the variation in the exposure is explained by genetic variants that act on that exposure through the confounder is moderate IVW and MR Egger are biased across all p-value selection thresholds, however both weighted median and weighted mode correctly estimate the causal effect of the exposure on the outcome (Figure 4, Supplementary Table 1). When almost all of the SNPs are associated with the confounder all univariable methods of estimation are biased and give results that are consistent with each other (Figure 4, Supplementary Table 1). With a stricter p-value cut-off, equivalent to a smaller sample size, MR Egger is very imprecise. This is expected as a stricter p-value cut off will mean fewer SNPs are included in the estimation and MR Egger has low power when the number of SNPs included in the estimation is small. In all cases both MVMR and the univariable MR approaches with Steiger filtering between the exposure and confounder applied gave unbiased estimates of the causal effect of the exposure on the outcome. However, each of these methods requires the confounder to be known (/suspected) and have GWAS summary statistics available.
In these simulations there is a single confounder of the exposure and outcome of interest biasing the estimates obtained. These results therefore show how the MR effect estimates are biased by that confounder in the same direction as the confounded observational association and may be substantially more biased than the observational association. In practice it is likely that there may be many potential confounders, some of which will be associated with a subset of the genetic variants selected as instruments. The overall bias in the linear regression will depend on the direction and magnitude of the bias from all of the confounders and so is likely to be substantially more than the bias in the MR estimation from confounding. The results obtained here would apply to each of the heritable confounders, and both Steiger filtering and MVMR could be applied to each of them in turn to identify and reduce that bias.
Table 1 shows the mean number of genetic variants in each simulation that were selected for the univariable MR estimation methods (IVW, MR Egger, weighted median and weighted mode). When the mediator is associated with as many SNPs as the confounder in the data generating process the final proportion of SNPs that are associated with the confounder in the MR estimation is small due to the effects of those SNPs on the exposure being smaller. Weighted median and weighted mode MR estimation methods are therefore unbiased as their assumptions regarding the proportion of SNPs unaffected by pleiotropy are satisfied. When most of the SNPs are associated with the confounder fewer SNPs are selected for estimation. However, most of them are associated with the confounder and so violate the IV assumptions leading to bias in all of the univariable MR methods of estimation.
Example: Heritable confounding in genetic variants associated with C-Reactive Protein
To illustrate how some GWAS may be detecting genetic variants associated with pleiotropic effects, we considered a single potential exposure that has been widely used as an exposure in MR studies; C-Reactive Protein (CRP), a marker of inflammation.[23-25] These studies have considered CRP using either all GWAS significant loci or using a cis MR approach, making use only of the significant loci that lie close to the CRP gene.[26] Cis-MR has the advantage of reducing the probability of high levels of pleiotropy however it has lower power than standard genome-wide MR approaches and limited ability to investigate how sensitive the results are to that pleiotropy. A recent GWAS has identified a large number of genetic variants associated with CRP across the genome, many of which may be associated with CRP due to their influence on traits which are upstream of CRP. We therefore explore the extent of heritable confounding that these additional SNPs may add to MR studies. We made use of the FUMA platform[27] to identify other traits that have been associated with the genomic loci in a recent large GWAS of CRP.[28]
The genetic signature of CRP was comprised of 728 lead SNPS (i.e., SNPs associated with CRP levels at p< 5×10−8 and LD r2 < 0.1) in 266 genomic loci and we identified 33,573 SNP -> trait associations among these lead SNPs or related SNPs in strong LD (LD r2 > 0.6) in the GWAS Catalog. The full list of these associated traits is given in the Supplementary material Section 2. Many of the SNPs identified as associated with CRP in the GWAS are likely to be associated with CRP through the traits we have identified. Therefore, these traits are potential heritable confounders which may bias any MR estimation which includes CRP as an exposure or outcome. We therefore selected 5198 traits for follow up analysis. These traits can be broadly categorised into: anthropometry measures related to height and weight, bone mineral density, glycaemic measures, proteins, metabolites and smoking behaviour traits (full list given in Supplementary Material Section 2). For each of these traits we estimated whether MR showed evidence of a causal effect of the trait on CRP. We applied Wald ratio for traits with a single SNP available as an instrument, inverse variance weighting for all traits with two or more SNPs available as instruments. Where possible (>3 SNPs) we also applied MR Egger, Weighted Median and Weighted Mode. Table 2 summarises the proportion of the analyses conducted that indicated evidence of a causal effect using a heuristic threshold to identify associations that indicate a potential causal effect (p<0.05). These results indicate that approximately 20% of the traits identified as associated with CRP-SNPs had an effect on CRP and so are potential confounders in any analysis including CRP. These effects were consistent across different pleiotropy robust MR estimators. These traits have the potential to bias any MR estimate in which CRP is an exposure through heritable confounding and indicate the importance of considering this form of bias.
Example: Estimation of the effect of age at menarche on risk of type 2 diabetes
We illustrate the potential for genetic variants associated with a confounder to bias MR estimation through a example estimating the causal effect of age at menarche on risk of type 2 diabetes. Here, we expect that adiposity confounds age at menarche and type 2 diabetes, and with a well-powered GWAS for age at menarche, some of the instruments will arise via body mass index and therefore MR estimation of this effect is likely to be biased. The hypothesised DAG is given in Figure 5 and results are given in Table 3. Univariable MR analysis shows an estimated causal effect of age at menarche on type 2 diabetes using IVW [odds ratio for a category increase in age at menarche on the risk of type 2 diabetes: 0.76, 95% confidence interval: 0.65 to 0.89]. However, this result is not consistent across the robust MR methods.
We also applied Steiger filtering in which 9 of 181 instruments for age at menarche appeared to primarily influence childhood adiposity were removed, leading to an attenuated effect estimate [OR: 0.90, 95% CI 0.79 to 1.02] (Figure 6). Including childhood adiposity in a MVMR more substantially attenuates the effect estimate obtained [OR: 0.99, 95% CI: 0.78 to 1.09]. As childhood adiposity is an established risk factor for both age at menarche and – through adulthood BMI - type 2 diabetes, these results are likely to reflect the model given in Figure 2b and suggest that the result seen in the naïve IVW estimation is likely to be biased by confounding from SNPs associated with childhood adiposity. As has been shown elsewhere, there is an overlap in the SNPs associated with childhood and adulthood adiposity.[29] As adiposity is a well-established risk factor for Type 2 Diabetes the SNPs associated with childhood adiposity will also be associated with Type 2 Diabetes through adulthood adiposity.
Discussion
In this paper we have shown how correlated pleiotropy due to heritable confounding can bias MR effect estimates when the genetic variants selected as instruments are associated with a confounder of the exposure and outcome of interest. We show that the bias in MR from including instruments associated with a heritable confounder could be greater than the bias in the linear regression from that confounder. In our network simulations we show, under a plausible model of biological effects, how genetic variants associated with confounders of an exposure and outcome may be associated with more traits upstream of the exposure than genetic variants which are associated with the exposure but not confounders of the exposure and outcome. The implication of this is that, as GWAS sample sizes increase these confounding genetic variants are increasingly likely to be detected as strongly associated with the exposure. This will increase the presence of correlated pleiotropy in MR effect estimates which use genetic variants strongly associated with the exposure in a GWAS as instruments for that exposure. Through simulations where there is a mixture of valid instruments or those acting via confounders, we show that bias caused by a specific confounder in MR estimation could be larger than the bias in the observation association caused by the same confounder and that this problem is likely to increase as sample sizes increase. Although it is likely that there will also be non-heritable confounding biasing the results from linear regression these results highlight how MR results may be increasingly biased as sample sizes for GWAS grow.
If a large proportion of the genetic instruments are associated with a confounder of the exposure and outcome, correlated pleiotropy will bias all of the standard pleiotropy robust methods of estimation. In this case each pleiotropy robust method is biased in the same way as the IVW estimator, potentially leading to false confidence in the results obtained. These results are likely to particularly apply to exposures that are downstream of a phenotype (such as BMI) which is associated with many genetic variants and has an influence on a wide range of traits, or where there are relatively few genetic variants directly associated with the exposure of interest. In these cases a high proportion of the genetic variants used may be associated with one or more heritable confounders, increasing the overall bias in the MR estimates.
We show that this bias from the inclusion of genetic variants associated with a confounder could be corrected through the application of either MVMR including the confounder, or Steiger filtering to remove any genetic variants that explain more variation in the confounder than the exposure. We therefore recommend that researchers applying MR with summary statistics generated from large GWAS studies check the association of their SNPs with the main potential confounders of the exposure and the outcome through a PheWAS analysis. Any genetic variant-trait associations found should then be checked for whether the trait involved is likely to be a confounder.[30] Steiger filtering, removing SNPs that explain more variation in the confounder than the exposure, and/or MVMR also including the confounder can then be applied to the MR analysis as an additional robustness test. A substantial difference between the results obtained is an indication of bias in the original MR analysis. However, in practice it is likely that confounders will not be identifiable/measured and there will be a number of confounders acting in this way and so such an approach does not guarantee the remove of all bias from heritable confounding.
Recently novel MR methods have been proposed that attempt to adjust for correlated pleiotropy.[14, 16, 18, 19, 31, 32] Each of these methods depends on different assumptions about the nature of the correlated pleiotropy and applies a different approach to correct for that pleiotropy. When standard MR methods are suspected to be biased by correlated pleiotropy then applying one or more of those approaches will help to identify whether correlated pleiotropy is present.
We have previously shown that adjustment for a potential confounder in the GWAS for the exposure will not mitigate bias in the MR estimates.[33, 34] This is because the adjustment for the confounder generates an inverse association between those genetic variants and the exposure. This biases the results obtained from the GWAS the opposite direction to the direction of effect of the genetic variant on the confounder. Previous work showed that this could be corrected by inclusion of the confounder the GWAS was adjusted for in a MVMR estimation[33] however all of the results described here would equally apply in that setting.
Where the exposure has a clear biological mechanism by which some genetic variants affect the exposure, comparing the results obtained for the MR estimation with all genetic variants to the results obtained for the genetic variants which have a known biological pathway can help to support the results or highlight the potential for bias in the MR estimate. Here we used the example trait of CRP to explore how many of the genetic variants detected in a recent GWAS were likely to be associated with confounders. These results showed a high level of pleiotropy among the variants identified, with many of the traits identified showing evidence of having an effect on CRP and so being potential confounders. These results highlight how prevalent heritable confounding bias is likely to be for many traits. CRP is an exposure for which alternative ‘cis-MR’ approaches could be used, making use only of the genetic variants identified within the CRP gene region. Where possible such an approach is likely to be more robust to heritable confounding than using the full set of genome-wide significantly associated variants and so both approaches should be used if possible.
In our illustrative application we estimated the effect of age at menarche on risk of Type 2 diabetes. The standard IVW effect estimate showed a negative causal effect with later age at menarche estimated to be protective for risk of type 2 diabetes. This effect however attenuated in all of the robust methods as well as when Steiger filtering or MVMR was applied to adjust for childhood adiposity. Childhood adiposity is an established risk factor for both age at menarche[22, 35, 36] and type 2 diabetes, through adulthood adiposity, [29, 37-40] and these analyses support the effect observed in the IVW estimation being driven by a set of SNPs that are associated with the confounder childhood adiposity being included as instruments for age at menarche.
The results given here highlight why it is important to consider whether the IV assumptions are likely to hold even when a range of MR estimates give consistent results for the effect of the exposure on the outcome. This is particularly true when there is a strongly genetically predicted potential confounder for the exposure and outcome, such as BMI, that has not been accounted for in the MR estimates or where the total number of genetic variants associated with the exposure of interest is likely to be small. This checking could be done in systematically through PheWAS analysis or through expert subject knowledge of the likely confounders.[30]
While these results show that the MR results can be more biased than linear regression from the exclusion of a single heritable confounder there are many other sources of bias for both types of estimation that should be considered when determining how much weight to give to the results from each type of analysis. It is likely that there will be a number of confounders of the exposure and outcome, if only some of these are heritable in the way described here then the MR estimates may still be less biased from confounding than the observational linear regression estimates. Additionally, MR can overcome bias and effect attenuation from measurement error that occurs in linear regression. However, MR is subject to other forms of bias such as population stratification and pleiotropy of other forms which should also be considered in any analysis.
A limitation of the results here is how we have set up the structure of relationships between genetic variants and phenotypes. In our network simulations we assumed that every edge has the same strength of effect. In our simulations of the bias of the different MR estimators we have assumed the same distribution in size of the association between the relevant genetic variants and the exposure and confounder. These assumptions are relatively strong and are not likely to be realistic for every trait. For some traits the confounders may have much stronger genetic effects than the exposure, which will accentuate the biases described here. For other traits the direct effect of the genetic variants on the exposure will be much larger than for the confounders, meaning that bias of the type described here will be less prevalent. Which of these is likely to hold, along with how strongly the confounders are likely to affect the outcome, will determine how much of a problem the results described here are likely to be in any particular case. Each MR study conducted should therefore consider individually how likely there is to be bias from heritable confounding in that setting and what the main confounders causing that bias may be.
Methods
Network simulations
We use the randomDAG function in the R/dagitty package to simulate a DAG with an arbitrary number of nodes, and a fixed probability that any set of nodes is connected. Nodes with no parents in the simulated DAG are identified as genetic variants, and 500 trait pairs are selected at random, in each pair assigning one to be the hypothesised exposure and one to be the hypothesised outcome. We then identify the path length for every genetic variant for the exposure, and also determine whether they act via a confounder, or otherwise (i.e. directly or through other traits that mediate the association between the genetic variant and the exposure and have no paths to the outcome). The effect size of the genetic variant on the exposure is βL where L is the path length from genotype to exposure, and for simplicity every edge in the DAG has the same effect β= 0.2. These simulations were conducted for a set of parameters of graph size N=75, 100, 125, 150 and graph density based on probability of an edge occurring ranging from p = 0.01, 0.02, …, 0.1. Each parameter combination was repeated 10 times.
Simulations evaluating MR methods
We generated the data in this model for two separate samples. In each sample 400 SNPs were generated using a binomial distribution. The proportion of these SNPs associated with the mediator and the confounder varied across the simulations, varying between 50% associated with each trait to 12.5% associated with the mediator and 87.5% associated with the confounder. The SNP-trait (i.e. mediator/confounder) association was drawn from a normal distribution with mean 0 and standard deviation of 0.04 for each SNP, orientated so all SNPs had a positive association with the relevant trait.
The confounder was generated as the sum of all the confounder - SNP effects plus a random error term. The mediator was generated as the sum of all mediator – SNP effects plus a random error term. The exposure was generated as the sum of the mediator, the confounder and a random error term. As shown earlier, SNPs associated with a mediator may be more strongly associated with the exposure than SNPs associated with a confounder. Therefore, we down weighted the confounder – exposure association by 25% relative to the mediator – exposure association (by multiplying them by 0.75) to mimic this more distant association. Results with this down weighting varied (to 50% and 0%) are given in the supplementary material. Finally, the outcome was generated as the effect of the confounder and a random error term.
In one of the samples generated we estimated the association between all of the SNPs and each of the exposure and confounder to generate GWAS summary statistics. This was repeated in the second sample for the SNPs and the outcome. We then select all SNPs (from either set) that were significantly associated with the exposure in the first sample at a defined p-value threshold in these summary statistics to estimate the effect of the exposure on the outcome using summary-data MR methods. Following what is currently standard practice in summary-data MR the same sample was used for the SNP-exposure association in estimation as had been used to identify the SNPs associated with the exposure. We applied four common univariable MR methods of estimation IVW[41], MR Egger[11], weighted median[12] and weighted mode[13]. We additionally estimated the model using multivariable MR (MVMR)[42, 43] including the confounder as an additional exposure, including all SNPs significantly associated with the confounder at the same p-value threshold. All univariable methods of estimation were repeated with Steiger filtering[44] applied between the exposure and confounder to remove SNPs that explain more variation in the confounder than the exposure.
With all other things fixed increasing the sample size will decrease the standard error in the SNP-trait association proportionally across all of the SNPs. This will lead to more SNPs being selected based on a fixed p-value threshold cut-off. Therefore, to model varying the sample size in the GWAS while keeping all other aspects of the data constant we varied the p-value threshold used to select SNPs for use as instruments in the MR across our simulations with larger p-values being equivalent to a larger sample size.
Simulation details
A set of 400 SNPs were generated as; For individual i and SNP j. The SNP exposure effects were generated as; The confounder was generated as; Where vc is a randomly normally distributed error term with mean 0 and standard deviation 1. njc is an SNP effects vector for the confounder where the elements of nj associated with the mediator (rather than the confounder) have been set to 0.
The mediator was generated as; Where vc is a randomly normally distributed error term with mean 0 and standard deviation 1. njrn is an SNP effects vector for the confounder where the elements of nj associated with the confounder (rather than the mediator) have been set to 0.
The exposure was generated as; Where vx ∼ N (0, 11).
The outcome was generated as; Where vy ∼ N (0, 11).
This data was generated for two different samples and summary statistics were generated by estimating the linear regressions; and For all j in the first sample and For all j in the second sample, to obtain estimates of βcj, βxj and βyj.
All models were simulated with a sample size of 100,000 for each of the exposure and outcome samples and 1000 repetitions.
Heritable confounding in genetic variants associated with C-Reactive Protein
C-reactive protein (CRP) is a widely studied acute phase hepatic protein and inflammation marker. Early observational studies consistently found associations between elevated CRP levels with increased cardiovascular disease (CVD) risk, motivating interest in targeting CRP as a therapeutic for CVDs.[45] In contrast with the observational findings, human genetics studies leveraging SNPs in the CRP locus revealed no direct causative links between CRP and cardiovascular diseases (CVDs),[28, 46, 47] indicating that CRP may not be the inflammatory driver of heart disease, which diminished the interest in developing CRP-targeted therapies.[46, 47] It appears that the associations between CRP and CVD observed in observational studies were likely influenced by confounding factors from other biomarkers and risk factors causally related to CVDs.[45, 48] More recent GWASs of circulating CRP levels have been conducted,[28, 49-51] uncovering hundreds of variants across the genome linked to CRP levels. Many of the loci are involved in metabolic and immune response pathways,[28] and the CRP GWAS demonstrates strong genetic correlation with cardiometabolic risk factors and biomarkers such as LDL cholesterol (LDL-C), high-density lipoprotein cholesterol (HDL-C), and triglycerides.[51] Correspondingly, MR studies instrumenting CRP levels with SNPs located throughout the genome have identified consistent, adverse relationships CRP and various CVD outcomes,[52-55] replicating the biased observational estimates and suggesting that heritable confounding with other CVD biomarkers and risk factors may be present.
We therefore sought to explore the pleiotropic relationships of the genomic loci associated with CRP to assess whether there is evidence for heritable confounding. SNPs associated with CRP levels were identified from in a recent GWAS of >500,000 individuals with European genetic ancestry.[28] This GWAS identified 266 loci as associated with CRP, 211 of which were novel. We made use of the FUMA platform[27] to identify other traits which have also been associated the 266 CRP genomic loci. We first identified the lead SNPs the CRP loci (independent SNPs LD r2 < 0.1, and GWAS P-values < 5×10−8) and performed a look up in the GWAS Catalog[56] (look up was performed on April 10 th,2024) on 728 lead SNPs and 3,077 independent SNPs in strong LD with the lead SNPs (GWAS P-values < 5×10−8 and LD r2 >0.1 with the respective lead SNP).
The resulting list of associated traits was curated to prioritise biomarkers and risk factors over diseases. From this curated list we conducted a summary-data MR analysis of each trait on to CRP levels. Where only one SNP was available as associated with the exposure trait of interest the Wald ratio was calculated. For traits with 2 or more SNPs the IVW MR effect estimate was obtained and for traits with 3 or more SNPs the MR Egger[11], Weighted Mode[13] and Weighted Median[12] estimates were also obtained.
Estimation of the effect of age at menarche on type 2 Diabetes
All data used were downloaded from MRC IEU OpenGWAS[57]. SNP – exposure associations were extracted for age at menarche for data on recalled age at which periods started for women in UK Biobank (n = 243,944) obtained through the IEU GWAS pipeline[58, 59]. SNP – exposure associations were obtained for childhood adiposity for a previously conducted GWAS using data on self-reported childhood adiposity for men and women in UKBiobank.[29] SNP – outcome associations for Type 2 diabetes was taken from a GWAS from FinnGen[60] to avoid sample overlap with UKBiobank.
We applied the univariable summary data estimation approaches; IVW, MR Egger, weighted median and weighted mode and additionally MVMR adjusting for BMI. For each approach we selected all SNPs associated with the relevant exposure at genome-wide significance (p-value for the SNP-exposure association < 5×10−8). For the MVMR we selected all SNPs associated with either age at menarche or childhood adiposity at the same threshold. Each univariable method was applied twice, once including all SNPs and once with Steiger filtering applied to remove SNPs that explained more variation in childhood adiposity than age at menarche. Steiger filtering was not applied to the MVMR estimation due to the inclusion of childhood adiposity as an exposure. Selected SNPs were then clumped to retain only independent SNPs (r2 < 0.001 based on a reference LD panel). All analyses were done in R using the ‘TwoSampleMR’ package with default parameters for SNP selection and clumping applied.[61]
Data Availability
All code used is available online at https://github.com/eleanorsanderson/confounding. Summary data used is available at https://gwas.mrcieu.ac.uk.
Data availability
Summary statistics used in all analyses are available at OpenGWAS study IDs ukb-b-3768, finn-b-E4_DM2, ukb-b-4650.
Code availability
The code for the simulations and analyses is available at https://github.com/eleanorsanderson/confounding.
Acknowledgements
ES, TP, GH, GDS, KT work in a unit funded by the MRC [MC_UU_00032/01, MC_UU_00032/02].