Interaction-based Mendelian randomization with measured and unmeasured gene-by-covariate interactions

Overview of MR-GxE and MR-GENIUS

The MR-GxE and MR-GENIUS approaches use differences in instrument strength across one or more covariates to estimate and correct for pleiotropic bias⁷. For i ∈ {1,2, …, N} observations, let G_i denote a single genetic instrument for an exposure X_i, and let Y_i represent the outcome of interest. Further, assume there exists an unmeasured confounder U_i of X_i and Y_i, and a set of interaction covariates Z_i = {Z₁ … Z_K} across which the instrument-exposure association varies. In order to make our ideas concrete, we now define an underlying data generating model for a continuous exposure and outcome, which are themselves a function of genetic variants G_i, Z_i and U_i.

In equations (1-4), the ϵ_(.i) terms represent independent error terms, and relationships with reference to a single interaction covariate Z_ki are illustrated in Figure 1 wherein G_i, Z_ki, and U_i are assumed independent for clarity.

• Download figure
• Open in new tab

Figure 1:

A directed acyclic graph showing the relationship between a genetic instrument G, interaction covariate Z, exposure X, outcome Y, and one or more confounders U. GZ denotes the interaction G × Z, and G, Z, and U are assumed independent for clarity.

The MR-GxE approach requires an interaction covariate (Z_i) to be explicitly observed, in contrast to MR-GENIUS which leverages variance differences for a given exposure (X_i) across subgroups of a genetic instrument (G_i). Such variances differences implicitly rely on the presence of one or more unmeasured interaction covariates. For both approaches, unbiased estimation of causal effects requires three assumptions to be satisfied, summarised as assumptions GxE1-3 below. A suitable interaction (G_iZ_i) is:

GxE1: Strongly associated with the exposure of interest.

GxE2: Independent of confounders of the exposure and outcome.

GxE3: Not directly associated with the outcome of interest.

For MR-GxE each assumption is framed with reference to a single observed interaction covariate. In contrast, MR-GENIUS requires assumptions GxE1-3 to hold globally across the set of observed and unobserved gene-by-covariate interactions in the sample (G_iZ_i). For example, GxE1 is satisfied with respect to MR-GENIUS when the average interaction strength for all possible interactions (γ_3K) is non-zero across the sample.

MR-GxE is implemented by using a gene-by-covariate interaction as an instrument within a two-stage least squares (TSLS) framework. In the first-stage model (equation 5), the exposure is regressed upon a genetic instrument and observed interaction covariate including an interaction term. The second-stage model (equation 6) then regresses the outcome upon the genetic instrument, interaction covariate, and fitted values for the exposure obtained using the first-stage model.

This returns a causal effect estimate , as well as an estimate of pleiotropic effect as the coefficient of the genetic instrument in the second-stage model (see Methods). MR-GENIUS is implemented by first performing a simple regression of the exposure upon the genetic instrument and obtaining a set of residuals . These residuals are then used in conjunction with the genetic instrument and incorporated within a TSLS model as a single instrument for the exposure (see Methods)¹².

MR-GxE sensitivity analyses and simulations

By explicitly defining an interaction-covariate, it is possible to perform a range of sensitivity analyses which evaluate assumptions GxE1-3. MR-GxE is reliant upon a strong first-stage interaction (GxE1), and where this is not the case estimates will exhibit weak instrument bias (see Methods). The first-stage F-statistic for the interaction term in the first-stage model can be used to identify candidate gene-by-covariate interactions, in a similar fashion to utilising GWASs to identify genetic variants associated with a phenotype of interest. Identifying candidate interactions is performed by fitting the first-stage MR-GxE model for each candidate interaction covariate Z_ki and calculating the F-statistic with respect to G_iZ_ki (see equation 5 and Methods)⁹. Applying a Bonferroni multiple testing correction, it is then possible to plot and identify interaction covariates for which the MR-GxE model is identified¹³.

To illustrate the role of interaction strength in MR-GxE and MR-GENIUS analyses, we present two simulated examples using the data generating models given in equations 1-5. Scenario 1 uses a single simulated dataset to demonstrate how sufficiently strong interactions can in principle be identified and visualised by evaluating the F-statistic for the interaction term in the first-stage model (equation 5). We use a sample size of 10,000 including K = 100 potential interaction covariates. Of the 100 interaction covariates, 10 were designated to have a first-stage interaction, assigning a value for γ_3k sampled from a normal distribution with mean 2 and standard deviation 0.5. Remaining interaction covariates were assigned values for γ_3k sampled from a normal distribution with mean 0 and standard deviation 0.01. The complete set of interaction covariates Z_Ki were generated so as to be independent of G_i, such that π_k1 = 0 in equation (1). Scenario 2 illustrates how weak instrument bias results in biased estimates of causal effect. This is achieved using the same data generating model as in simulation 1, varying the strength of the first-stage interaction to demonstrate the effect of weak instrument bias. Causal effect estimates represent the mean estimate for a given mean F-statistic across a total of 5000 iterations, and were obtained using a separate sample of 10,000 observations.

The results of simulations 1 and 2 are presented in Table 1 and Figure 2 below. For simulation 1, Figure 2a shows how a scatter plot can be constructed in a similar fashion to a Manhattan plot in GWAS analyses. In this case, a Bonferroni multiple testing correction can be visualised using a solid horizontal line. From Table 1 it can be seen that using individual interaction covariates within the MR-GxE framework provides results comparable to MR-GENIUS when the assumptions of both approaches are satisfied. Simulation 2 demonstrates how causal effect estimates exhibit both bias and a loss of precision as interaction strength decreases, as shown using a forest plot in Figure 2b and presented in Table 1.

• Download figure
• Open in new tab

Figure 2:

Panel A shows a scatter plot of − log₁₀(p − value) for the first-stage F-statistic across the set of 100 potential interaction covariates for a single simulated data set. Panel B shows a forest plot of mean causal effect estimates and confidence intervals under varying mean interaction strength. A solid horizontal line is included representing the Bonferroni correction threshold for statistical significance in panel A, while the dotted vertical line in panel B represents the true causal effect β₁ = 1.

MR-GENIUS requires the set of unobserved interactions to be sufficiently strong globally to provide sufficiently precise estimates of causal effect. This can be shown in a further simulated example under scenario 1 above, with the addition that the proportion of non-zero interactions is varied within the sample. Table 2 contrasts estimates obtained using a single valid and explicitly defined interaction using MR-GxE with MR-GENIUS estimates implicitly using the entire interaction set. It can be seen that as the mean interaction strength increases across the set of K covariates the precision of the MR-GENIUS approach improves. This highlights the utility of the MR-GxE approach where interaction covariates are readily identifiable, and the flexibility of the MR-GENIUS approach in the absence of observed gene-by-covariate interactions.

Assumption GxE2 requires the gene-by-covariate interaction to be independent of common causes and confounders of the exposure and outcome. This is equivalent to instrument exogeneity in the conventional MR setting, though it should be noted that associations between the instrument G_i and the interaction covariate Z_i do not necessarily introduce bias into MR-GxE estimates. This is most likely the result of problematic confounding structures between G_i, Z_i, and U_i, and consequently it is possible to estimate the correlation between G_i and Z_i, with independence serving as evidence of GxE2 being satisfied. However, it should be highlighted such independence does not necessary imply that GxE2 is satisfied, but an observed correlation can highlight potential issues that warrant further consideration (see Methods).

The third assumption GxE3 requires pleiotropic effects of G_i upon Y_i to remain constant across values of Z_i, with the gene-by-covariate interaction being independent of Y_i when conditioning on X_i. Where this is not the case estimates of causal effect will exhibit bias in the direction of β₄ in a similar fashion to horizontal pleiotropic bias in univariate MR analyses. By reframing MR-GxE within a TSLS framework, it is possible to apply tests of over-identification to evaluate the constant pleiotropy assumption, though this is not possible where only one instrument is available, for example, a single genetic variant. In cases where the single instrument is comprised of many instruments, such as a polygenic risk score, it is possible to examine different configurations of instruments iteratively using MR-GxE and assess heterogeneity in the set of MR-GxE estimates obtained from each iteration. These subsets of instruments are hereafter referred to as sub-instruments.

In this scenario, a Sargan test can be used to compare different MR-GxE estimates of the same causal parameter (the coefficient of X_i in equation 6 – i.e., β₁), assuming we have more instruments than we need to consistently estimate the parameter¹⁴. However, it is important to note that in applying this test it is crucial for each of the sub-instruments to be sufficiently strong to overcome weak instrument bias, though practically the test can be applied where weak interactions are present if assessing the strength of individual instruments of interest.

To demonstrate the impact of GxE3 violation, as well as the utility of employing an adapted Sargan test as a sensitivity analysis, we present a further simulation shown in Table 3. In this case, a score analogous to a PRS was used as a single IV, comprised of 10 individual sub-instruments of approximately equal strength. Mirroring the previous simulated example, the true causal effect was defined as β₁ = 1 with a horizontal pleiotropic effect β₂ = 0.05. Sub-instruments violating assumption GxE3 were estimated to have a value β₄ = 0.2, varying the proportion of invalid sub-instruments across the set of simulated models. The mean F-statistic for the simulations is 698.5 (Breusch-Pagan 253.97 p < 0.001), and MR-GENIUS estimates are presented for comparison.

As shown in Table 3, both MR-GxE and MR-GENIUS produce biased causal effect estimates when the constant pleiotropy assumption is violated. Violation of the constant pleiotropy assumption is also detected by applying a Sargan test, provided all sub-instruments do not identically violate GxE3 such that β_4k is constant. As the Sargan test relies upon at least one instrument being valid, identical violation of assumption GxE3 would also violate the assumptions of the conventional Sargan approach.

Estimating the effect of adiposity on systolic blood pressure within the UK Biobank

To demonstrate each of the sensitivity analyses previously described, we performed MR analyses estimating the causal effect of adiposity (measured using BMI) on SBP using data from the UK Biobank. This serves as a re-examination of the original applied example in Spiller et al (2019) who first proposed the MR GxE model⁷. Here we go further by evaluating each underlying assumption using the diagnostic tools described above, and contrasting the results with MR-GENIUS^{7, 12}. After performing quality control, removing participants with missing data, and restricting the sample to unrelated individuals of European ancestry, a total of 358,928 participants were included in the analyses.

MR-GxE was implemented by constructing a weighted PRS informed using genetic variants previously identified from the GIANT consortium¹⁵. As the GIANT consortium represents a subset of the most recent UK Biobank release, subsequent analyses have been conducted in a one-sample framework. A total of 95 independent genetic variants were used after performing linkage disequilibrium (LD) pruning, and removing tri-allelic or palindromic variants. Finally, we standardized BMI, SBP, and the weighted PRS using a z-score transformation prior to performing analyses. In previous work we found evidence of a positive association between BMI and SBP using OLS and TSLS regression approaches^{7, 16-18}.

Initially, a discovery subset (N=100,000) was randomly sampled from the UK Biobank data for use in identifying interactions for MRGxE analyses. Causal effect estimates and sensitivity analyses were performed using the remaining data. Candidate gene-by-covariate interactions were detected by estimating the first-stage F-statistic for 576 candidate interaction covariates within the UK Biobank. After applying a multiple testing correction, the 20 interaction covariates with the strongest association were selected and utilised in subsequent analyses. Table 4 shows MR-GxE estimates of causal effect and corresponding sensitivity analyses with respect to each interaction covariate. The strength of each interaction across the set of candidate interaction covariates is illustrated in Figure 3, where annotations give the UK Biobank field ID for each interaction covariate.

• Download figure
• Open in new tab

Figure 3:

A scatter plot showing the first-stage F-statistics for instrument-by-covariate interactions using data from UK Biobank. A horizontal line is included representing the Bonferroni correction for statistical significance. For clarity, blue points represent interactions identified after multiple testing. The 20 strongest interactions have been annotated using their UK Biobank field identification number.

To assess assumption GxE2, we created 9 sub-instruments sampling from the 95 SNPs used to create the initial PRS instrument. Fitting the MR-GxE model using multiple sub-instruments allows for overidentification tests to be performed, testing the extent to which causal effect estimates differ when using individual sub-instruments. In each case, a failure to reject the null can be considered to be evidence of interaction exogeneity as previously outlined. To implement this approach, the set of SNPs were randomly assorted into 9 sub-instruments of approximately equal strength, quantified using the F-statistic with respect to BMI. Repeating this procedure using sub-instruments containing differing SNPs yielded similar results. We also present the mean F-statistic across the set of sub-instruments to emphasise their strength.

As shown in Table 4, there exists substantial disagreement across the range of selected interaction covariates, suggesting that one or more violate underlying assumptions of the MR-GxE approach.

Considering assumption GxE2, several of the identified gene-by-covariate interactions are proxy measures of adiposity, specifically waist circumference, weight in kilograms, and comparative body size at age 10. Such interaction covariates are often problematic, as associations between the genetic variants and the interaction can result in collider bias where the interaction covariate is downstream of the exposure (see Methods). In this case, higher estimates of ρ (G, Z) for these variables supports this interpretation, and their subsequent exclusion from further analyses. A similar argument can also be made with respect to interaction covariates downstream of BMI, including diabetes diagnosis, vascular/heart problem diagnosis, and diastolic blood pressure (DBP).

By applying Sargan tests, a number of interaction covariates related to cognition, physical activity, and age appear to violate assumption GxE3. This could be explained by the gene-by-covariate interactions relating to one or more underlying risk factors, which are not adjusted for in the corresponding MR-GxE models.

After applying sensitivity analyses, three interaction covariates can be identified as appropriate choices for estimation using MR-GxE. This selection was made using Sargan test and correlation p-value thresholds of < 0.0025, applying a multiple testing correction. Selected covariates include alcohol intake frequency and physical activity, both days walked and moderate levels of exercise. Considering alcohol intake and physical activity, the lack of a substantial correlation between each interaction covariate and the PRS suggests that violation of GxE2 is unlikely.

In previous work Townsend deprivation index (TDI) was selected as an interaction covariate in a summary MR-GxE analysis and returned estimates in agreement with both alcohol consumption and physical activity measures identified above. However, it is important to note that TDI shows evidence of a non-zero instrument-interaction covariate correlation, potentially highlighting a violation of assumption GxE2. This can be explained by TDI being plausibly downstream of both BMI and the instrument, representing situation in which the correlation does not invalidate estimates of causal effect. For reference, this would resemble confounding scenario (c) (see Methods Figure 5c) with respect confounding structures invalidating assumption GxE2.

Crucially, adopting alcohol and physical activity as interaction covariates yields causal effect estimates which appear biologically plausible, and support evidence from both observational and MR studies suggesting a positive association between BMI and SBP. Estimates using each interaction covariate are presented in Figure 4.

• Download figure
• Open in new tab

Figure 4:

A forest plot showing MR-GxE causal effect estimates using the interaction covariates presented in Table 4. Observation f.4253.0.5 has been omitted for clarity. Red points indicate analyses for which assumptions may likely be violated, while blue points show potentially valid interaction covariates using accompanying sensitivity analyses. Observational, two-stage least squares (TSLS), and MR-GENIUS estimates are also shown as black points.

As a final analysis, we implemented MR-GENIUS using the PRS, BMI, and SBP measures from UK Biobank. This resulted in a more precise estimate in comparison to MR-GxE, however, the effect estimate appears to strongly disagree with evidence from MR-GxE and alternate approaches. Given MR-GENIUS implicitly relies upon analogous assumptions to MR-GxE, it seems reasonable to assume that such a discrepancy could arise from bias due to violations stemming from one or more unmeasured interactions. This is further supported by MR-GxE estimates of similar direction and magnitude which appear to show evidence of bias, such as vigorous physical activity which shows evidence of GxE3 violation.

The MR-GxE Approach

The extent to which instrument strength varies across strata of a given interaction covariate Z_k is quantified by the magnitude of the first-stage interaction between G_i and Z_ki (γ_3k) shown in equation (3). The MR-GxE approach also relies upon β_4k and θ_3k being equal to zero, analogous to instrument validity in conventional MR. In previous work, MR-GxE was implemented using an approach analogous to MR-Egger regression in two-sample summary MR^{4, 7}. This is achieved by initially obtaining strata specific instrument-exposure and instrument-outcome associations, after which the instrument-outcome associations are regressed upon the instrument-exposure associations including an intercept.

The interpretation of estimates and corresponding plots using this approach essentially mirrors MR-Egger regression, though each observation represents a different strata-specific subgroup as opposed to a unique genetic instrument⁴. However, whilst this approach does allow for estimation using publicly available GWAS summary data it has two primary limitations. Initially, it relies upon specification of an interaction-covariate without a readily applicable measure of interaction strength. Second, ambiguity surrounding the definition of interaction-covariate strata can potentially have a substantial impact on causal effect estimates⁷. Applying MR-GxE in an individual level data setting avoids these issues and allows for additional information to be incorporated into effect estimation. Further, such an approach is increasingly applicable with the emergence of large-scale studies with genetic data such as the UK Biobank.

A reduced form model for Y_i given G_i and Z_i incorporating equations (5-6) can be written as

Note that a G_i × Z_1i term is omitted from the second-stage model given in equation (6) due to its role as an instrument, whilst the inclusion of G_i allows for estimation of a horizontal pleiotropic effect on the outcome, denoted by β₂. Using equation (5) and equation (7) we can define the MR-GxE estimand as

The MR-GENIUS Approach

The MR-GENIUS approach is an adapted form of Robins’ G-estimation which is robust to additive confounding and pleiotropic bias^{12, 21, 22}. This essentially involves leveraging differences in the variance of a given exposure X_i across subgroups of a genetic instrument G_i, which are likely the consequence of one or more gene-by-covariate interactions. In the case of a binary instrument and exposure, and using notation from equations 1-5, the MR-GENIUS estimator can be written as: where and ¹².

MR-GENIUS is implemented by first regressing X_i upon G_i and obtaining a set of residuals . These residuals are then used to create an instrument which is incorporated within a TSLS model, as a single instrument for X_i¹². Estimates of remain unbiased, provided the instrument G_i is associated with the exposure of interest, the effect does not change across values of the unmeasured confounders, and the MR-GENIUS model is identified such that the change in variance across levels of the instrument is non-zero¹². In the binary exposure case, this means that the MR-GENIUS model is identified when cov(G_i, var(X_i|G_i)) ≠ 0, and for a continuous exposure the MR-GENIUS model when the residual error ϵ_X is heteroskedastic, that is, not constant across levels of G_i¹². This can be evaluated using a Breusch-Pagan test for heteroskedasticity^{12, 19}.

These conditions also restrict the degree of joint effect modifiers of both X_i and Y_i. Importantly, it should be noted that the interaction covariate need not be explicitly identified using MR-GENIUS, illustrated by the absence of any Z_ki in equation 9. However, identification of the MR-GENIUS model implicitly relies upon the presence of one or more gene-by-covariate interactions to induce the desired dependence between G_i and var(X|G).

GxE1: Interaction strength

The MR-GxE estimator can be viewed as an extension of the Wald ratio, including an adjustment for the direct effects of G_i and Z_ki. Thus, in the special case where G_i and Z_ki are marginally independent of the exposure and outcome (but their interaction via a single covariate Z_ki is not), the MR-GxE estimator collapses to:

MR-GxE is clearly reliant upon a strong first-stage interaction, such that γ_k3 ≠ 0 in order to make the denominator of (10) non-zero (GxE1). As in the conventional MR setting it is possible to assess interaction strength using the F-statistic for the interaction term in the first-stage model, however, the use of a single instrument precludes relating the F-statistic to the magnitude of relative bias (at least three instruments would be required in this case for the asymptotic formula to be valid). Therefore, whilst an F-statistic of 10 may satisfy the standard threshold for sufficient instrument strength, it is not possible to relate this to a 10% relative bias towards the OLS estimate. Further, as in the conventional MR setting, interaction strength does not mitigate bias from violations of assumptions GxE 2-3²³. Where possible, candidate interactions should be identified in separate samples to avoid issues related to Winner’s curse, as is the case with instrument selection in conventional MR.

At this point it is important to highlight several features of gene-by-covariate interactions which require careful consideration prior to performing MR-GxE. Firstly, interactions are scale dependent, and as a result applying transformations can create spurious associations²⁴. As such spurious associations can exist as an artefact of the data, estimates leveraging such information can potentially be unreliable. A related concern is that the interaction may not be linear, as is assumed in equation 5. A potential solution to this issue is to fit flexible models (e.g., fractional polynomial models, which include varying exponents with respect to G_iZ_ki) to allow for non-linear interactions to be identified²⁵.

As MR-GENIUS does not require gene-by-covariate interactions to be identified, testing for identification is performed globally by evaluating heteroskedasticity with respect to the residuals ϵ_Xi. Specifically, MR-GENIUS relies upon the residual error in a regression of the exposure upon the instrument to be heteroskedastic, such that .

GxE2: Interaction exogeneity

In previous work we show how assumption GxE2 is most likely violated when certain confounding structures exist, specifically, where G_i and Z_i are simultaneously downstream of a confounder U_i or where there is an open path between the two variables through U_i⁷. These confounding structures are illustrated in Figure 5, where scenarios (a), (b), and (d) represent pathways which can result in bias. In contrast, scenario (c) represents a blocked path which does not necessarily result in bias.

• Download figure
• Open in new tab

Figure 5:

A panel showing a set of 4 DAGs wherein the instrument G and interaction covariate Z are simultaneously associated with one or more confounders U. In this case, DAGs (a), (b), and (d) would likely result in biased causal effect estimates using MR-GxE, in contrast to confounding structure (c).

To understand why scenarios (a)-(c) induce bias into MR-GxE estimates, we can extend the MR-GxE estimand to allow for violation of GxE 2, by including covariance terms between U_i and (Z_i, G_iZ_i)U_i, such that were it possible to include U_i in the TSLS model, the resulting estimate could be written as: where each β^∗indicates a multivariable regression estimate pertaining to the second subscript variable when regressed upon the first, including the unmeasured confounder U_i. As it is not possible to directly measure and adjust for U_i, we rely upon independence between U_i and G_iZ_i for equation (11) to be equivalent to the MR-GxE estimator in equation (8).

Violation of GxE2 and subsequent bias can result from specific configurations of G_i, Z_i, and U_i associations. For example, a pathway from G_i to Z_i through U_i (see Figure 5a) would result in bias, though the converse is less likely where genetic variants are used as instruments⁷.

A more serious concern is the potential for collider bias when estimating fitted values in the first-stage MR-GxE model. As shown in equation (5), it is necessary to include the interaction covariate in the first-stage model. However, in cases where G_i and U_i are both simultaneously upstream associated with Z_i, conditioning on Z_i will induce collider bias in the first-stage MR-GxE model, such that the estimate of pleiotropic effect and subsequent adjustment will be inaccurate. This case is illustrated in Figure 6.

• Download figure
• Open in new tab

Figure 6:

A DAG showing a situation in which conditioning on Z when G and X are simultaneously upstream associated with Z would induce collider bias in a regression of X upon G and Z outlined in equation 5.

We present two strategies for limiting the impact of such bias. As an initial test, it is possible to estimate the correlation between G_i and Z_i, with observed independence serving as evidence against violation of GxE2. Specifically, with reference to Figure 5 scenarios (a), (b), and (d) can in principle lead to observed correlations between G_i and Z_i which can be identified. However, it is important to emphasise that independence cannot necessarily be interpreted as GxE2 being satisfied. This would primarily be the case where a three-way interaction exists between the instrument G_i, interaction covariate Z_i, and one or more confounders U_i. Rather than removing the possibility, an observed correlation between G_i and Zi can highlight a potential issue in the analysis which warrants further consideration.

A second and potentially more robust approach would be to adopt a genetic proxy variable for the interaction covariate Z_i, as this would share the same benefits with regard to causal direction as G_i with respect to environmental confounders. For example, when estimating the association between alcohol and SBP using education as an interaction covariate, adopting a PRS for education would in principle utilise the explained variation in education excluding environmental confounders such as socio-economic status.

GxE3: Constant pleiotropy

The third MR-GxE assumption requires pleiotropic effects of G_i upon Y_i to remain constant across values of Z_i, with the gene-by-covariate interaction being independent of Y_i when conditioning on X_i. Where this is not the case estimates of causal effect will exhibit bias, such that the degree of bias in the MR-GxE estimate will be equal to:

To illustrate how over-identification tests can be applied in the context of MR-GxE, consider an extension of equations 5 and 6 to include an arbitrary number of sub-instruments, wherein a single instrument G_i is comprised of m ∈ {1,2, …, M} sub-instruments. Where G_mi denotes the m^th sub-instrument in G_i, we can rewrite the data generating models presented in equations (3-4) as:

Sargan tests can then be used to compare different MR-GxE estimates of the same causal parameter provided we have more instruments of sufficient strength than we need to consistently estimate the parameter¹⁴.