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

Studies leveraging gene-environment (GxE) interactions within Mendelian randomization (MR) analyses have prompted the emergence of two similar methodologies: MR-GxE and MR-GENIUS. Such methods are attractive in allowing for pleiotropic bias to be corrected when using individual instruments. Specifically, MR-GxE requires an interaction to be explicitly identified, while MR-GENIUS does not. We critically examine the assumptions of MR-GxE and MR-GENIUS in the absence of a pre-defined covariate, and propose sensitivity analyses to evaluate their performance. Finally, we explore the effect of body mass index (BMI) upon systolic blood pressure (SBP) using data from the UK Biobank, finding evidence of a positive effect of BMI on SBP. We find both approaches share similar assumptions, though differences between the approaches lend themselves to differing research settings. Where a suitable gene-by-covariate interaction is observed MR-GxE can produce unbiased causal effect estimates. MR-GENIUS can circumvent the need to identify interactions, but as a consequence relies on either the MR-GxE assumptions holding globally, or additional information with respect to the distribution of pleiotropic effects in the absence of an explicitly defined interaction covariate.

Introduction

Mendelian randomization (MR) is an epidemiological approach applied to observational data, wherein genetic variants are used as instrumental variables (IVs) to estimate the effect of a modifiable exposure on a downstream outcome [1]. MR encompasses a wide range of statistical methods, and typically relies upon three assumptions to test for causality. A suitable genetic IV is strongly associated with the exposure of interest (IV1), independent of confounders of the exposure and outcome as well as confounders of the genetic IV and outcome (IV2), and independent of the outcome when conditioning on the exposure (IV3) [1, 2].

Violation of assumptions IV2–3 can introduce bias into MR effect estimates, and as a consequence methods for identifying and correcting for such bias have formed a central theme within the MR methods literature [2-4]. In many cases, such methods focus upon correcting for bias resulting from associations between a genetic IV and an outcome which are unrelated to the exposure of interest, defined as horizontal pleiotropic pathways [5]. Pleiotropy robust methods frequently use heterogeneity in causal effect estimates across multiple genetic IVs as an indicator of horizontal pleiotropy, though such approaches are less feasible as the number of available genetic IVs decreases [6, 7].

One solution to the problem of limited available genetic IVs is to leverage variation in instrument strength across one or more covariates within a target population, representing a gene-by-covariate interaction [8-10]. Intuitively, were it possible to identify a population subgroup for which a genetic IV and exposure are independent (i.e., a ‘no-relevance group’), it follows that, in the absence of horizontal pleiotropy, the genetic IV and outcome should also be independent. A non-zero instrument-outcome association for such a group would therefore be indicative of pleiotropic bias [8, 11, 12]. It is, however, rare that no-relevance groups of sufficient size are observed in practice.

MR approaches utilising gene-by-covariate interactions, here referred to as interaction-MR, overcome this limitation by using statistical assumptions to extrapolate back to a hypothetical no-relevance group. Two such methods are MR using Gene-by-Environment interactions (MR-GxE) and MR G-Estimation under No Interaction with Unmeasured Selection (MR-GENIUS) [8, 13]. MR-GxE uses an explicitly defined gene-by-covariate interaction to estimate causal effects, and has previously been framed within a summary-level data context [8]. In contrast, MR-GENIUS accommodates both observed and unobserved interactions, provided they induce a dependence between the genetic IV and exposure variance [13]. MR-GENIUS has the advantage of circumventing the need to explicitly identify gene-by-covariate interactions, though the relative strengths and limitations of the approach compared to MR-GxE have previously been unclear.

In this paper we outline the implementation of MR-GxE in the individual level data setting, and critically evaluate the performance of MR-GxE and MR-GENIUS. Specifically, we focus upon the application of MR-GENIUS in the absence of a pre-defined interaction covariate, which is not possible using MR-GxE. Through simulation we demonstrate how both approaches share similar underlying assumptions, and highlight how implicitly leveraging all potential gene-by-covariate interactions using MR-GENIUS can imply more stringent assumptions with respect to the distribution of pleiotropic effects. Throughout we also propose sensitivity analyses to test the assumptions of the MR-GxE. Finally, we conduct applied analyses using MR-GxE and MR-GENIUS to estimate the effect of body mass index (BMI) on systolic blood pressure (SBP). For both approaches we find evidence of a positive causal effect using data from the UK Biobank, comparing results to conventional MR and observational methods.

Materials and methods

The data generating model

Interaction-MR approaches use differences in instrument strength across one or more covariates to estimate and correct horizontal pleiotropic bias [8]. For i ∈ {1, 2, …, N} observations, let G denote a single genetic IV for an exposure X, and let Y represent the outcome of interest. Further, assume there exists an unmeasured confounder U of X and Y, and a set of interaction covariates Z ∈ {Z1, …, Z} 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 G, Z and U.

In Eqs 1–4, the ϵ(. terms represent independent error terms, and relationships with reference to a single interaction covariate Z are illustrated in Fig 1 wherein G, Z, and U are assumed independent for clarity.

Fig 1

Illustration of data generating model.

A directed acyclic graph showing the relationship between a genetic instrument G, an 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.

An overview of MR-GxE and MR-GENIUS

The MR-GxE and MR-GENIUS approaches rely upon one or more first-stage interactions which induce variation in the association between the genetic IV and exposure. Specifically, the MR-GxE approach requires an interaction covariate (Zi) to be explicitly observed, in contrast to MR-GENIUS which leverages variance differences for a given exposure (X) across subgroups of a genetic IV (G). In this paper we illustrate how both approaches are reliant upon three assumptions, summarised as assumptions GxE1–3 below. A suitable interaction (GZ) is:

Strongly associated with the exposure of interest (GxE1).

Independent of confounders of the exposure and outcome (GxE2).

Not directly associated with the outcome of interest (GxE3).

MR-GxE was originally implemented using an approach analogous to MR-Egger regression in two-sample summary MR [7, 8]. Initially sets of instrument-exposure and instrument-outcome associations are obtained across strata of a pre-specified interaction covariate, after which the instrument-outcome associations are regressed upon the instrument-exposure associations including an intercept [8]. While in principle the approach can be performed using publicly available data from genome-wide association studies (GWAS), the summary level MR-GxE approach has two notable limitations. First, summary MR-GxE does not readily provide a means of evaluating interaction strength, relying on observed heterogeneity across gene-exposure associations across interaction covariate strata [8]. Second, ambiguities surrounding the optimal number of interaction covariate strata can have a substantial impact of effect estimates [8]. To address these issues, we propose an individual-level form of MR-GxE within a two-stage least squares (TSLS) framework.

Individual level MR-GxE is implemented by using a gene-by-covariate interaction as an instrument within a TSLS regression model. In the first-stage model (Eq 5), the exposure is regressed upon the genetic IV and observed interaction covariate including an interaction term (γ3). The second-stage model (Eq 5) then regresses the outcome upon the genetic IV, interaction covariate, and fitted values for the exposure () obtained using the first-stage model.

This returns a causal effect estimate (), as well as a horizontal pleiotropic effect estimate as the coefficient of the genetic IV () in the second-stage model. To define the MR-GxE estimand, a reduced form model for Y given G and Z incorporating Eqs 5 and 6 can initially be written as

Using Eqs 5 and 6 the MR-GxE estimand is then defined as

Note that a G × Z term is omitted from the second-stage model given in Eq 6 due to its role as an instrument, whilst the inclusion of G allows for estimation of a horizontal pleiotropic effect on the outcome, denoted by β2.

The MR-GENIUS approach is an adapted form of Robins’ G-estimation which is robust to additive confounding and pleiotropic bias [13-15]. This essentially involves leveraging differences in the variance of a given exposure X across subgroups of a genetic instrument G, 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 Eqs 1–4, the MR-GENIUS estimator can be written as: where and [13].

MR-GENIUS is implemented by first regressing X upon G 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 [13]. Estimates of remain unbiased, provided the instrument G 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 [13].

In the binary exposure case, the MR-GENIUS model is identified when cov(G, var(X|G)) ≠ 0, and for a continuous exposure when the residual error ϵ is heteroskedastic, that is, not constant across levels of G [13]. This can be evaluated using a Breusch-Pagan test for heteroskedasticity, and these conditions also restrict the degree of joint effect modifiers of both X and Y [13, 16].

Importantly, it should be noted that the interaction covariate need not be explicitly identified using MR-GENIUS, illustrated by the absence of Z in Eq 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 and var(X|G). In the absence of a predefined interaction covariate, MR-GENIUS estimates the total effect of X upon Y, without adjusting for the interaction covariate Z. This contrasts with MR-GxE, which estimates the direct effect of X upon Y adjusting for the interaction covariate in the second stage model.

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 and Z. Thus, in the special case where G and Z are marginally independent of the exposure and outcome (but their interaction via a single covariate Z is not), the MR-GxE estimator simplifies to:

From Eq 10 MR-GxE is clearly reliant upon a strong first-stage interaction, such that γ ≠ 0 in order to make the denominator non-zero (GxE1). When individual-level data are available, the first-stage F-statistic for the gene-by-covariate interaction can be used to quantify instrument strength, though several aspects of this approach warrant consideration. First, when using a single interaction the F-statistic cannot be related to the magnitude of relative bias towards the observational estimate in a one-sample setting and null in a two-sample setting. This is because such a relationship between instrument strength and the direction of bias only holds when multiple instruments, in this case interactions, are used. Therefore, while an F-statistic of 10 may satisfy the standard threshold for sufficient instrument strength, it would not be possible to relate this to a 10% relative bias towards the observational estimate obtained by regressing the outcome on the exposure without incorporating additional interaction covariates. Second, interaction strength does not mitigate bias from violations of assumptions GxE2–3, just as is the case in conventional MR analyses. Finally, where possible candidate interactions should be identified in separate samples to avoid issues related to Winner’s curse, where instruments, in this case interactions, are selected using spurious associations which may be sample population specific [17].

The reliance of MR-GxE upon explicitly defined interactions also invites two potential interaction-specific issues: scale dependency and non-linear interactions. First, as interactions are scale dependent it is possible that applying transformations can create spurious associations [18]. Such spurious associations can exist as an artefact of the data, and consequently estimates leveraging such information are unlikely to be reliable. Gene-by-covariate interactions may also be non-linear, which could potentially be considered by fitting more flexible models (e.g., fractional polynomial models, which include varying exponents with respect to GZ) to allow for non-linear interactions to be identified. It is, however, important to take care to avoid issues of over-fitting [19].

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 ϵ. Specifically, MR-GENIUS relies upon the residual error in a regression of the exposure upon the genetic IV to be heteroskedastic, evaluated using a Breusch-Pagan test for heteroskedasticity [13].

As a means of identifying candidate gene-by-covariate interactions for MR-GxE we propose using the first-stage F-statistic for the interaction term in the first stage, in a similar fashion to utilising GWASs to identify genetic variants associated with a phenotype of interest. Interactions of sufficient strength can be identified by fitting the first-stage MR-GxE model for each candidate interaction covariate Z and calculating the F-statistic with respect to GZ (see Eq 5) [10]. Applying a Bonferroni multiple testing correction, and plotting the −log10(p − value) for the F-statistic then allows for instrument strength to be effectively visualised using a scatter plot, following a similar intuition to the use of Manhattan plots in the presentation of results from GWAS [20]. Note that as it is often the case that multiple independent genetic variants are associated, it is often appropriate to use a polygenic risk score as an instrument to maximise instrument strength.

GxE2: Interaction exogeneity

In previous work we show how assumption GxE2 is potentially violated when certain confounding structures exist, specifically, where G and Z are simultaneously downstream of a confounder U or where there is an open path between the two variables through U [8]. To briefly recapitulate how such associations can induce bias, consider the path diagram shown in Fig 2. In this case, the interaction covariate Z is independent of X and Y, and determined by a confounder U. Further, U is downstream associated with the genetic instrument G.

Fig 2

An example of GxE2 through confounding.

A path diagram illustrating a case in which the instrument G is a determinant of the interaction covariate Z through a confounder U. The bidirectional dashed arrow from Z to Y represents an association induced due to confounding as a result of not adjusting for U in the second stage MR-GxE model.

In Fig 2U serves not only as a confounder of X and Y, but also of Z and Y. As the MR-GxE model only instruments G, it is likely estimates for the effect of Z on Y will exhibit bias. When G is not independent of U, however, the resulting induced association between Z and Y from failing to control for U mimics a pleiotropic association, such that a pathway from G to Y is created through U and Z. Importantly, such associations do not necessarily bias estimates of the effect of X on Y, but can inflate type-I error rates when evaluating instrument validity.

Assumption GxE2 can also be violated when a gene-by-covariate interaction is simultaneously associated with the exposure and one or more confounders of the exposure and outcome, as depicted in Fig 3.

Fig 3

Illustration of general GxE2 violation.

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. In this case, the presence of an association between a gene-by-covariate GZ and U violates assumption GxE2.

In Fig 3 a bidirectional arrow is included to highlight that any direction of association between GZ and U can potentially introduce bias into MR-GxE estimates. Where GZ is upstream of U this can be viewed as instrument strength varying across levels of the confounder, with pleiotropic effects being associated with interaction strength. This issue can be viewed as analogous to the INstrument Strength Independent of Direct Effect (InSIDE) assumption in two-sample summary MR [7]. An association from U to GZ would suggest that a three-way interaction may be present, such that interaction strength γ3k varies across levels of U. This can bias effect estimates by inducing an association between GZ and Y, violating the constant pleiotropy assumption GxE3.

To understand how an association between GZ and U can induce bias into MR-GxE estimates, we can extend the MR-GxE estimand (Eq 8) to incorporate violation of GxE2, by including covariance terms between U and (Zi, GZi)U, such that were it possible to include U 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. As it is not possible to directly measure and adjust for U, the independence U and GZ is relied upon for Eq 11 to be equivalent to the MR-GxE estimator in Eq 8.

A further consideration is the introduction of collider bias when estimating fitted values in the first-stage MR-GxE model. As shown in Eq 5, it is necessary to include the interaction covariate in the first-stage model. However, in cases where G and U are both simultaneously upstream associated with Z, conditioning on Z 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 Fig 4.

Fig 4

Illustration of collider bias when estimating .

A diagram showing a situation in which conditioning on Z when G and U are simultaneously upstream associated with Z would induce collider bias, as shown by the dashed bidirectional arrow.

Relating assumption GxE2 to the MR-GENIUS approach, associations violating GxE2 would imply associations vary across values of the unmeasured confounders violating the second MR-GENIUS assumption [13]. However, this problem can be mitigated by incorporating additional interaction covariates within the MR-GENIUS model, as described in Eric Tchetgen et al., 2021 [13]. This would necessitate the inclusion of specific interaction covariates within the MR-GENIUS model, such that differences in the variance of X would be evaluated across subgroups of G, conditional on one or more interaction covariates Z.

For MR-GxE we present two strategies for addressing GxE2 violation. To evaluate the possibility of collider bias in the first-stage model estimating the correlation between G and Z could serve as an initial test for GxE2 violation. Intuitively, if G and Z are independent, then conditioning on Z would not induce an association between G and U. 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, interaction covariate Z, and one or more confounders U. Rather than removing the possibility, an observed correlation between G and Z can highlight a potential issue in the analysis which warrants further consideration.

A potentially more robust approach would be to adopt a genetic proxy variable for the interaction covariate Z, as this would share the same benefits with regard to causal direction as G with respect to environmental confounders. For example, when estimating the effect of alcohol consumption on SBP using education as an interaction covariate, adopting a polygenic risk score (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 upon Y to remain constant across values of Z, with the gene-by-covariate interaction being independent of Y when conditioning on X (i.e. β4 = 0). Where this is not the case estimates of causal effect will exhibit bias in the direction of β4 in a similar fashion to horizontal pleiotropic bias in univariate MR analyses, equal to:

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 PRS, 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 in Eq 6—i.e., β1), assuming we have more instruments than we need to consistently estimate the parameter [21]. 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 illustrate how over-identification tests can be applied in the context of MR-GxE, consider an extension of Eqs 5 and 6 to include an arbitrary number of sub-instruments, wherein a single instrument G is comprised of m ∈ {1, 2, …, M} sub-instruments. Where G denotes the m sub-instrument in G, we can define a corresponding data generating model as:

A Sargan test can be applied by fitting multiple sub-instruments G in the same TSLS model. Alternatively, a heterogeneity test such as Cochran’s Q-statistic could be used to evaluate heterogeneity in MR-GxE effect estimates using all sets of non-overlapping sub-instruments.

Results

An illustration of assumptions GxE1–3 through simulation

To illustrate the importance of assumptions GxE1–3 with respect to MR-GxE and MR-GENIUS we present six simulation studies, categorised by assumption, using the data generating model presented in Eqs 1–4. Throughout, we demonstrate the utility of the sensitivity analyses proposed, and highlight the relative performance of both MR-GxE and MR-GENIUS. In each simulation gene-by-covariate first-stage effects (γ3) are generated to be positive, to avoid the possibility that the combined effects of all candidate interactions have a mean of zero. This would potentially invalidate the MR-GENIUS approach in the unlikely event that leveraged candidate interactions have effects such that cov(G, var(X|G)) ≈ 0. Code for performing each simulation study and further information is available at https://github.com/WSpiller/GxE_Simulation.

Simulation set 1: Interaction selection and strength

As an illustration of how gene-by-covariate interactions can be identified through the evaluation of their first-stage F-statistics, we generated 1, 000 independent data sets, containing 100, 000 observations for a single instrument G, exposure X, outcome Y, and 100 candidate interaction covariates Z (Simulation 1). All variables were treated as continuous, with observations of exogenous variables randomly sampled from a normal distribution with mean 0 and standard deviation 1. Endogenous variables, determined by one or more additional covariates, were generated following the population models defined in Eqs 1–4, with error terms randomly sampled from a normal distribution with mean 0 and standard deviation 1. The effect of X upon Y was defined as β1 = 1. Of the 100 interaction covariates, 10 were designated to have a non-zero first-stage interaction, assigning a value for γ3 sampled from a normal distribution with mean 2 and standard deviation 2, ensuring that all coefficients for non-zero first stage interactions were greater than 1. The complete set of interaction covariates Z were also generated so as to be independent of G, such that π = 0. Fig 5A shows how a scatter plot can be constructed in a similar fashion to a Manhattan plot in GWAS analyses. Each value on the scatter plot represents the mean −log10(p − value) value for the first-stage F-statistic corresponding to each candidate interaction across the set of 1, 000 simulated data sets. A Bonferroni multiple testing correction is shown using a solid horizontal line.

Fig 5

Plots corresponding to simulations 1–2, identifying interactions and visualising the impact of weak instrument bias for MR-GxE.

Panel A shows a scatter plot of −log10(p − value) for the mean first-stage F-statistic across the set of 100 potential interaction covariates in simulation 1. A solid horizontal line is included representing the Bonferroni correction threshold for statistical significance in panel A. Panel B shows a forest plot of mean causal effect estimates and confidence intervals under varying mean interaction strengths in simulation 2. The dotted vertical line in panel B represents the true causal effect β1 = 1, and arrows are used to indicate confidence intervals exceeding the limits of the forest plot.

In Fig 5A the 10 defined non-zero gene-by-covariate interactions have been identified, with super-imposed numbers indicating the identity of each interaction covariate Z. The corresponding estimates for the 10 identified interactions show no evidence of apparent bias, with a mean MR-GxE estimate of 1.000 (95% CI = 0.996, 1.004) and a mean F-statistic of 3616.15 (p − value < 0.001). Individual mean estimates for each interaction are provided in the Supplementry material (see S1 Table). Using MR-GENIUS resulted in mean estimate of 1.000 (95% CI = 0.994, 1.006), producing an estimate comparable to MR-GxE without the need to explicitly identify an interaction covariate. Performing a Breusch-Pagan test for identification in the MR-GENIUS model yielded a mean value of 1041.74 (p − value < 0.001), suggesting MR-GENIUS estimates are sufficiently strong so as to overcome weak instrument bias.

To further investigate the impact of weak instrument bias using MR-GxE we perform additional simulations, evaluating the performance of MR-GxE using a single non-zero interaction covariate of varying strength (simulation 2). Specifically, first-stage F-statistics of approximately 1, 5, 10, 25, 50, and 100 are considered, generating 1, 000 data sets for each F-statistic value and presenting mean MR-GxE effect estimates and 95% confidence intervals. In each case the genetic instrument G is generated so as to satisfy assumption IV1 (γ1 = 1), with a causal effect of X on Y again equal to 1 (β1 = 1). Fig 5B shows a forest plot including the mean MR-GxE effect estimate and 95% confidence interval for each interaction covariate with a mean F-statistic as indicated on the y-axis. The precision of MR-GxE increases substantially as the mean F-statistic increases, and there does not appear to be evidence of directional bias using weak interactions.

In simulation 1 MR-GENIUS appears to perform well when many gene-by-covariate interactions are present, with the potential to outperform MR-GxE when individual stronger interactions are not observed (see S1 Table). To explore the extent to which MR-GENIUS is reliant on a global non-zero mean first-stage interaction, an additional simulation is conducted varying the proportion of non-zero interactions present in the data (simulation 3). A total of K = 100 interactions were generated, such that the number of non-zero interactions represent 1%, 5%, 10%, 50%, and 100% of all candidate interaction covariates. For each predefined proportion, 1,000 independent data sets were generated, using previous parameter definitions from simulation 1. MR-GxE effect estimates were obtained using a single randomly sampled non-zero interaction covariate, while MR-GENIUS estimates do not specify an observed interaction covariate. The mean MR-GxE and MR-GENIUS estimates for each proportion are presented in Table 1.

Table 1

Simulated results using differing proportions of non-zero interaction covariates (simulation 3).

Proportionγ3K ≠ 0	Mean MR-GxE β^_1(95% CI)	MR-GENIUS β^_1(95% CI)	MeanF-statistic	BP-Testp-value
1%	1.000 (0.99,1.01)	1.002 (-2.34,4.34)	6.14	0.202
5%	1.000 (0.99,1.01)	0.999 (0.96,1.04)	23.68	0.007
10%	1.000 (0.99,1.01)	1.000 (0.98,1.02)	39.25	0.002
50%	1.000 (0.99,1.01)	1.000 (0.99,1.01)	76.51	<0.001
100%	1.000 (0.98,1.02)	1.000 (0.99,1.01)	87.90	<0.001

From Table 1 it appears the precision of MR-GENIUS estimates improves as the mean interaction strength across all leveraged instruments () increases in magnitude, indicated by the increase in mean F-statistic across the set of candidate interaction covariates. This suggests that in cases where few gene-by-covariate interactions of moderate strength are available, MR-GxE can furnish more precise estimates than MR-GENIUS, though MR-GENIUS can outperform MR-GxE as the number of non-zero gene-by-covariate interactions increases. It is also important to highlight that, just as was the case for MR-GxE, weak instrument strength for MR-GENIUS does not appear to induce observed directional bias.

This reliance of MR-GENIUS upon mean interaction strength across candidate interactions has two important implications. First, the precision of MR-GENIUS estimates is a function of both interaction strength and the number of candidate non-zero interactions present, and consequently it is possible for MR-GENIUS to outperform MR-GxE as the number of strong gene-by-covariate interactions increases. Second, it is possible that interactions of similar magnitude acting in opposite directions can counteract each other, such that cov(G, var(X|G)) ≈ 0. This scenario is unlikely to occur beyond a simulation setting, and motivates using a data generating model where first-stage interactions are generated so as to be in the same direction, to ensure .

Simulation set 2: Interaction exogeneity

The MR-GxE approach relies upon the selected gene-by-covariate interaction being independent of all confounders U of the exposure X and outcome Y. As previously discussed, this is most likely to be the case where either the association between G and U varies across levels of Z, or the association between Z and U varies across levels of G. To demonstrate how such associations can introduce bias into both MR-GxE and MR-GENIUS effect estimates, we present a simulation with a similar structure to simulation 3, in this case varying the proportion of interactions for which θ3k ≠ 0 (simulation 4). Specifically, proportions of 0%, 1%, 5%, 10%, 50%, and 99% are considered for which θ3k = 1, generating 1, 000 independent data sets for each proportion. For each data set, the mean MR-GENIUS and MR-GxE estimates were obtained across all interaction covariates. Additionally, a mean estimate using a randomly sampled interaction for which θ3k ≠ 0 is also presented, to illustrate how assumption GxE2 is interaction covariate specific for MR-GxE. The simulation results are given in Table 2.

Table 2

Simulated results using differing proportions of non-zero gene-by-covariate interaction with respect to confounders (simulation 4).

Proportionθ3K ≠ 0	Valid MR-GxEaβ^_1(95% CI)	Mean MR-GxE β^_1(95% CI)	MR-GENIUS β^_1(95% CI)	MeanF-statistic
0%	1.000 (0.998,1.002)	1.000 (0.990,1.010)	1.000 (0.995,1.005)	67.62
1%	1.000 (0.998,1.002)	1.003 (0.990,1.016)	1.001 (0.995,1.008)	67.71
5%	1.000 (0.997,1.004)	1.014 (0.993,1.035)	1.008 (0.997,1.019)	68.17
10%	1.000 (0.995,1.005)	1.027 (0.999,1.055)	1.015 (1.001,1.030)	68.94
50%	1.000 (0.990,1.010)	1.136 (1.088,1.184)	1.070 (1.042,1.097)	72.80
99%	1.000 (0.987,1.014)	1.268 (1.268,1.304)	1.122 (1.090,1.153)	76.49

a Valid MR-GxE is used to indicate MR-GxE estimates obtained using a single interaction covariate for which GxE2 is satisfied.

From Table 2 both MR-GxE and MR-GENIUS exhibit bias in the direction of the interaction coefficient θ3. MR-GENIUS appears to be more robust to GxE2 violation compared to MR-GxE, though estimates decrease in precision as the magnitude of θ3 increases. When selecting a single interaction covariate for which θ3 = 0, MR-GxE provides an unbiased causal effect estimate, in contrast to MR-GENIUS. This can be explained by MR-GENIUS implicitly relying upon when an interaction covariate is not specified. As a consequence, MR-GxE appears to be capable of producing results with markedly less bias using a single GxE2 satisfying interaction, compared to MR-GENIUS where .

Simulation set 3: Constant pleiotropy

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 two simulated examples. Initially, we generated data as in simulations 3–4, instead varying the proportion of candidate interactions with non-zero second stage interactions (β4 ≠ 0) (simulation 5). The first-stage interaction coefficient across all interaction was set to γ3 = 1, with designated invalid interactions having a second stage interaction(β4 = 1). For each proportion the mean MR-GENIUS estimate was obtained, as well as the mean MR-GxE estimate across all candidate interactions. In addition, a mean MR-GxE estimate using a single randomly sampled interaction for which β4 = 0 is also provided, with results presented in Table 3.

Table 3

Simulated results using differing proportions of non-zero gene-by-covariate interaction with respect to the outcome (simulation 5).

Proportionβ4K ≠ 0	Valid MR-GxEaβ^_1(95% CI)	Mean MR-GxE β^_1(95% CI)	MR-GENIUS β^_1(95% CI)	MeanF-statistic
0%	1.000 (0.998,1.002)	1.000 (0.990,1.010)	1.000 (0.995,1.005)	67.62
1%	1.000 (0.998,1.002)	1.004 (0.990,1.017)	1.001 (0.995,1.008)	67.64
5%	1.001 (0.997,1.004)	1.019 (0.995,1.042)	1.008 (0.997,1.019)	67.66
10%	1.000 (0.995,1.005)	1.036 (1.004,1.069)	1.016 (1.001,1.031)	67.69
50%	1.000 (0.990,1.010)	1.185 (1.117,1.254)	1.082 (1.054,1.110)	67.52
99%	0.997 (0.987,1.014)	1.366 (1.300,1.432)	1.163 (1.131,1.194)	67.63

aValid MR-GxE is used to indicate MR-GxE estimates obtained using a single interaction covariate for which GxE3 is satisfied.

In this scenario it can be seen that MR-GxE and MR-GENIUS exhibit bias in the direction of GxE2 violation, while utilising a single interaction for which GxE3 is satisfied provides unbiased estimates. This would suggest that MR-GENIUS is reliant upon across all implicitly leveraged interactions.

To further 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 4 (simulation 6). In this case, a score analogous to a PRS was used as a single IV, comprised of 1, 000 individual sub-instruments of approximately equal strength. Mirroring the previous simulated example, the true causal effect was defined as β1 = 1 with a horizontal pleiotropic effect β2 = 0.05. Sub-instruments violating assumption GxE3 were estimated to have a value β4 = 0.2, varying the proportion of invalid sub-instruments. The mean F-statistic across all iterations simulations was 98.80 (Breusch-Pagan 31.80, p − value = 0.013), and MR-GENIUS estimates are presented for comparison.

Table 4

Simulated results illustrating use of Sargan test to identify GxE3 violation (simulation 6).

Proportiona β^_4k≠0	MR-GxE β^_1(95% CI)	MR-GENIUS β^_1(95% CI)	MeanF-statistic	Sarganp-value
0%	1.000 (0.997, 1.003)	1.000 (0.963,1.036)	98.95	0.480
1%	1.010 (1.006, 1.014)	1.005 (0.956,1.054)	98.95	<0.001
5%	1.050 (1.045, 1.056)	1.022 (0.940,1.104)	98.88	<0.001
10%	1.100 (1.093, 1.110)	1.049 (0.937,1.161)	98.92	<0.001
50%	1.499 (1.485. 1.514)	1.256 (1.012,1.551)	98.95	<0.001
100%	2.000 (1.981, 2.020)	1.505 (1.149,1.861)	98.77	0.456

aProportion β4 ≠ 0 refers to the proportion of sub-instruments which violate assumption GxE3.

As shown in Table 4, 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. 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. The UK Biobank obtained written consent from all participants, and received ethical approval from the Research Ethics Committee (REC reference for UK Biobank is 11/NW/0382). This analysis was approved by the UK Biobank access committee as part of project 8786. Consent was sought by UK Biobank as part of the recruitment process. This serves as a re-examination of the original applied example in Spiller et al. (2019) who first proposed the MR-GxE model [8]. The UK Biobank has approval from the North West Multi-centre Research Ethics Committee (MREC) as a Research Tissue Bank (RTB) approval, and consequently separate ethical clearance was not required for this project which was conducted under the RTB approval. In this study we evaluate each underlying assumption using the diagnostic tools described above, and contrasting the results with MR-GENIUS [8, 13]. 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 [22]. 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 [8, 23–25].

Initially, a discovery subset (N = 100, 000) was randomly sampled from the UK Biobank data for use in identifying interactions for MR-GxE 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 5 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 Fig 6, where annotations give the UK Biobank field ID for each interaction covariate.

Table 5

MR-GxE estimates and sensitivity analyses using each candidate interaction covariate and MR-GENIUS.

Covariate (UK Biobank Field ID)	F-Statistic	β^_1 (p-value)	ρ(G, Z)a(p-value)	Sarganb(p-value)	MeanFc
Waist circumference(f.48.0.0)	182.86	-0.524(<0.001)	0.103(<0.001)	7.531(0.481)	79.40
Weight (kg)(f.21002.0.0)	123.16	-0.687(<0.001)	0.119(<0.001)	9.342(0.314)	48.79
Diabetes diagnosis(f.2443.0.0)	54.22	-0.065(0.470)	0.020(<0.001)	12.19(0.143)	41.22
Alcohol intake frequency(f.1558.0.0)	50.65	0.163(0.006)	0.001(0.526)	5.69(0.682)	41.62
Physical activity (vigorous)(f.904.0.0)	42.10	0.017(0.862)	0.004(0.003)	20.40(0.009)	17.03
Vascular/ heart problem diagnosis(f.6150.0.0)	33.65	-0.446(<0.001)	0.028(<0.001)	7.22(0.513)	16.76
Time number displayed during memory test(f.4253.0.5)	28.42	-2.155(0.333)	0.015(0.002)	14.54(0.069)	13.51
Number of days per week walked 10+ mins(f.864.0.0)	27.87	0.208(0.011)	0.001(0.705)	6.87(0.551)	18.98
DBP (automated, baseline)(f.4079.0.0)	26.45	-0.324(<0.001)	0.020(<0.001)	6.39(0.603)	16.32
Physical activity (moderate)(f.884.0.0)	23.60	0.165(0.107)	0.001(0.324)	3.66(0.886)	14.27
Townsend deprivation index(f.189.0.0)	23.01	0.108(0.489)	-0.016(<0.001)	9.24(0.323)	16.80
Comparative body size at age 10(f.1687.0.0)	20.65	0.283(0.004)	0.048(<0.001)	9.94(0.269)	14.62
Time to complete pair matching activity(f.400.0.2)	20.49	0.052(0.689)	-0.007(<0.001)	29.50(<0.001)	11.84
Pulse rate(f.4194.0.0)	20.45	0.031(0.873)	-0.010(<0.001)	13.78(0.088)	4.77
Time watching television(f.1070.0.0)	20.01	-0.140(0.211)	0.017(<0.001)	14.83(0.063)	15.08
DBP (automated, follow-up)(f.4079.0.1)	19.55	-0.501(<0.001)	0.016(<0.001)	6.37(0.606)	11.72
Own or rent accommodation(f.680.0.0)	18.41	0.078(0.607)	-0.006(<0.001)	19.84(0.011)	10.13
Age at assessment(f.21003.0.0)	18.15	0.697(<0.001)	0.013(<0.001)	17.28(0.027)	14.03
Birthweight known(f.120.0.0)	17.93	0.067(0.851)	-0.014(<0.001)	14.16(0.078)	4.94
Year of birth	15.85	0.710(<0.001)	-0.014(<0.001)	17.12(0.029)	13.94
OLS	-	0.186 (<0.001)	-	-	-
TSLS	7776.52	0.130 (<0.001)	-	-	-
MR-GENIUS	1332.7 (<0.001) d	0.034 (0.009)	-	-	-

aρ(G, Z) represents the correlation between the PRS and interaction covariate,

b Sargan shows the results to over identification tests using sub-instruments,

c The mean F-statistic for sub-instruments,

d BP Heterogeneity Test.

Fig 6

Identified gene-by-covariate interactions with respect to genetically predicted body mass index.

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 GxE3, 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 over identification 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 5, 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 Materials and 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 p-value<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.

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 Fig 7.

Fig 7

A forest plot showing MR-GxE causal effect estimates using the interaction covariates presented in Table 5.

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.

Discussion

In this paper we examine two related interaction-based MR approaches: MR-GxE and MR-GENIUS. Both MR-GxE and MR-GENIUS rely upon similar underlying assumptions, whilst differing based on whether a gene-by-covariate interaction needs to be explicitly incorporated within the estimation model. Specifically, MR-GxE relies upon at least a single measured gene-by-covariate interaction which satisfies assumptions GxE 1–3, whilst MR-GENIUS does not require such an interaction to be observed. However, as a consequence of implicitly leveraging multiple underlying interactions, the MR-GENIUS approach requires assumptions GxE 1–3 to hold globally. Essentially, stronger assumptions are required to mitigate the absence of an observed gene-by-covariate interaction. It should be emphasised however, that evaluation of MR-GENIUS in this paper does not consider the inclusion of observed gene-by-covariate interactions.

Through an examination of the MR-GxE assumptions, several approaches aiming to evaluate assumptions GxE 1–3 have been outlined. Interaction strength (GxE1) can be evaluated using the first-stage F-statistic for the interaction term, analogous to evaluating instrument strength in conventional MR. The corresponding global test for interaction strength using MR-GENIUS and a continuous exposure is the Breusch-Pagan test for heteroskedasticity [13-16].

Assumption GxE2 can initially be evaluated by estimating the correlation between Z and both G and X respectively. Where Z is observed to be correlated with G, it is possible that a confounding relationship exists violating assumption GxE2. Further, the simultaneous association of Z with G and X can result in bias where Z is downstream of X. However, as the existence of such correlations does not necessarily imply that this assumption is violated, a more promising approach may be to adopt an interaction covariate Z which is highly likely to be exogenous (see Materials and methods). For example, one could employ genetic variants which instrument a likely interaction covariate. Future work will explore this possibility.

The constant pleiotropy assumption (GxE3) can be tested in cases where the initial instrument G is a composite instrument, that is, comprised of multiple sub-instruments such as genetic variants within a PRS. Heterogeneity in effect estimates obtained using sub-instruments can be considered as evidence of violation of the constant pleiotropy assumption, analogous to heterogeneity in two-sample summary MR [7, 26]. In principle, a similar approach can be applied using sub instruments with MR-GENIUS, though such an examination is beyond the scope of this paper. A summary of the MR-GxE assumptions and proposed tests is given in Table 6.

Table 6

A summary of the MR-GxE assumptions and proposed sensitivity analyses.

Assumption	Description	Consequence of violation	Tool to assess plausibility
Interaction strength (GxE1)	An observed gene-by-covariate interaction GZ should be selected, such that the association between the instrument G and exposure X varies across levels of Z	Insufficient precision to detect causal effects and directional bias when multiple interactions are used.	Estimating the first stage F-statistic for GZ and adopting an interaction covariate such that F ≥ 10.
Interaction exogeneity (GxE2)	The gene-by-covariate interaction GZ should be independent of confounders of the exposure X and outcome Y.	Inflated type-I error rates when evaluating instrument validity and biased effect estimates for the effect of X on Y.	Estimating the association between the instrument G and interaction covariate Z, selecting an interaction such that G and Z are independent.
Constant pleiotropy (GxE3)	The direct effect of an instrument G on the outcome Y should remain constant across levels of the interaction covariate Z.	Estimates of the effect of the exposure X on Y will be biased in the direction of the effect of GZ on Y.	Using a Sargan test when sub-instruments can be constructed from a composite instrument G.

In the applied analysis the effect of BMI on SBP was estimated using MR-GENIUS and a range of interaction covariates in conjunction with MR-GxE. We identified three suitable interaction covariates, which suggest a positive effect of BMI upon SBP in agreement with previous observational and MR analyses. Importantly, we highlight interaction covariates which violate the MR-GxE assumptions and link these issues to the possibly biased effect estimates obtained using MR-GENIUS.

Several limitations remain with respect to MR-GxE which warrant further explanation. Firstly, reliance upon an observed gene-by-covariate interaction limits the extent to which the method can be applied in contrast to MR-GENIUS. We advocate the use of MR-GENIUS in cases where no interaction covariate is available, though care needs to be taken in justifying the more stringent assumptions MR-GENIUS entails if an interaction covariate is not specified. Second, evaluating GxE2 using the correlation of between Z and G does not provide a clear indication of whether the assumptions hold. It is possible that GxE2 can be violated when Z and G appear to be independent, and assuming the direction of effect between Z and X relies upon a priori knowledge regarding the direction of association. It is therefore critical to identify plausible biological mechanisms underpinning the observed relationships in the MR-GxE model.

Finally, whilst an overidentification test has been presented for evaluating GxE3, there is not at present a method aiming to correct for violation of the constant pleiotropy assumption. It is likely that pleiotropy robust methods, such as median or modal regression, could be utilised to correct for resulting bias, and the application of such methods will be fully explored in future work.

Conclusion

MR-GxE and MR-GENIUS are two interaction-based MR approaches which leverage gene-by-covariate interactions to estimate causal associations, while correcting for instrument invalidity. MR-GxE can be adapted to the individual level data setting and allows for the underlying assumptions of the approach to be tested provided a gene-by-covariate interaction is explicitly identified. In contrast, MR-GENIUS does not require such an interaction to be identified, but instead relies upon a more stringent set of assumptions analogous to MR-GxE. The use of each method should therefore reflect the specific research questions considered, as each approach is especially suited to particular research contexts. However, it is essential that the strengths and limitations of each approach are given sufficient consideration prior to their application.

Simulated results and effect estimates for subset of interaction (denoted Z) identified from Fig 5A (simulation 1).

(PDF)

Click here for additional data file.

Transfer Alert

This paper was transferred from another journal. As a result, its full editorial history (including decision letters, peer reviews and author responses) may not be present.

26 Apr 2022

3 Jun 2022

Please see the formatted letter submitted with the manuscript. The contents of the letter are also included below.

Many thanks,

Response to reviewers

Dear Reviewers,

Many thanks for your very kind and careful consideration of our manuscript. Please find responses to each of your comments below, which we hope will explaining how we have attempted to address the issues highlighted.

Journal requirements

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

Thank you for highlighting these potential issues. The author summary has been removed, and the file names have been changed to correspond to the journal guidelines.

2. Please review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript. If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

I do not believe we have cited any retracted papers, and to my knowledge all papers have been cited correctly.

3. Please provide additional details regarding participant consent. In the ethics statement in the Methods and online submission information, please ensure that you have specified what type you obtained (for instance, written or verbal, and if verbal, how it was documented and witnessed). If your study included minors, state whether you obtained consent from parents or guardians. If the need for consent was waived by the ethics committee, please include this information.

Additional information as been provided regarding written consent obtained from each participant. The study is conducted using data collected by the UK Biobank, and no additional data was obtained for this manuscript.

4. Thank you for stating the following financial disclosure: "Wes Spiller is supported by a Wellcome Trust studentship (108902/B/15/Z)."Please state what role the funders took in the study. If the funders had no role, please state: "The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript."

The suggested statement has been included in the funding section for resubmission.

5. In your Data Availability statement, you have not specified where the minimal data set underlying the results described in your manuscript can be found. PLOS defines a study's minimal data set as the underlying data used to reach the conclusions drawn in the manuscript and any additional data required to replicate the reported study findings in their entirety. All PLOS journals require that the minimal data set be made fully available. For more information about our data policy, please see http://journals.plos.org/plosone/s/data-availability. Upon re-submitting your revised manuscript, please upload your study’s minimal underlying data set as either Supporting Information files or to a stable, public repository and include the relevant URLs, DOIs, or accession numbers within your revised cover letter. For a list of acceptable repositories, please see http://journals.plos.org/plosone/s/data-availability#loc-recommended-repositories. Any potentially identifying patient information must be fully anonymized. Important: If there are ethical or legal restrictions to sharing your data publicly, please explain these restrictions in detail. Please see our guidelines for more information on what we consider unacceptable restrictions to publicly sharing data: http://journals.plos.org/plosone/s/data-availability#loc-unacceptable-data-access-restrictions. Note that it is not acceptable for the authors to be the sole named individuals responsible for ensuring data access. We will update your Data Availability statement to reflect the information you provide in your cover letter.

Many thanks for highlighting potential issues surrounding data availability. As stated in the data availability section, the data used in simulations is available from https://github.com/WSpiller/GxE_Simulation and the code produces the data sets evaluated for each simulation.

The UK Biobank data unfortunately cannot be made publicly available following the legal requirements for access to the data. We have included the application number under which the analyses were performed, in addition to a link to information for researchers who wish to apply for access. Once access is granted through the UK Biobank system, researchers will have the necessary data to replicate our analyses.

Reviewer comments

Reviewer #1

The paper discusses the relative strengths and limitations of interaction-MR approaches for both MR-GxE and MR-GENIUS methods, which is a great contribution to the field.

Many thanks for your kind words regarding our manuscript.

1. For simulation settings, could you please provide more details how the variables in equations are generated (i.e. what are the distributions they come from and the associated parameters)?

Thank you for highlighting this issue. Additional information has now been added to page 3, describing how each variable is generated. Each variable is treated as continuous, with each observation either randomly sampled from a normal distribution with mean 0 and standard deviation of 1, or as a combination of prespecified covariates and an error term randomly sampled from a normal distribution with mean 0 and standard deviation of 1.

2. For simulation results, could you make it clear what type of outcomes they are (seems the results are for continuous exposure and outcome).

This is tied with the previous comment, where there was some ambiguity regarding the generated variables. We have now stated that all variables are continuous for simplicity. Thank you for highlighting this oversight.

3. Have you tried to extend the simulations to other types of outcomes (i.e. binary outcome)?

This is a very interesting question, as often binary outcomes are of interest in clinical settings. In the most simple setting, a researcher could fit the linear model to a binary exposure using a linear probability model, though that can have it’s own problems (interpretability, predicted values less than 0 or greater than 1 etc.)

We believe that the MR-GxE model could be adapted into a two-stage predictor substitution model, which would essentially use a first stage linear model, and then incorporate the fitted values within a second stage logistic regression model. There is no reason to believe this wouldn’t work, but there are two reasons why we haven’t explored this possibility in the paper.

First, MR-GENIUS is typically applied to continuous outcomes, while allowing for binary instruments and exposures. As one of the main focuses of the paper is a comparison between the MR-GxE and MR-GENIUS approaches, it was thought best to conduct all analyses within the continuous outcome setting. Second, we wanted to keep the manuscript to an appropriate length for publication, and evaluating binary outcome MR-GxE models with sufficient detail would increase the length of the paper substantially. As the potential issues and body of work related to binary outcomes in IV analyses (and binary exposures), we thought it best that a detailed discussed be deferred to future work.

Reviewer #2

This is an exceptionally well written manuscript that I thoroughly enjoyed reading. The background material in the Materials and Methods was very useful - providing intuition for understanding the mechanics of the approaches and without sacrificing intellectual rigour. I do not have any major concerns with the paper, but do have a number of suggestions for improving its clarity in certain parts, and some additional points that the author might want to consider.

Thank you for your very kind words and comments on the paper.

Line 110 Page 5: The notation is confusing here. I assume this is meant to be epsilon_hat_subscript_X_i but it could also be read as epsilon_hat multiplied by X_subscript_i.

You are correct, this is a typo which we have now amended.

Line 138 Page 5: “First, when using a single interaction the F-statistic cannot be related to the magnitude of relative bias, as at least three instruments would be required for the asymptotic formula to be valid.”. This sentence is not clear to me.

Thank you for highlighting the need for greater clarification on this point. In conventional IV analyses there is great interest in the consequences of using weak instruments to estimate causal effects. In MR where we can often have multiple instruments, it is possible to determine the direction of weak instrument bias based on whether instrument-exposure and instrument-outcome associations are obtained from non-overlapping samples. This is true, however, only in instances where multiple instruments (or for MR-GxE interactions) are used.

Admittedly this is more of a technical point however, as the use of weak instruments will typically reduce precision to the point that the direction of weak instrument bias is of marginal interest. It is something we mention because it is highlighted as a useful feature of MR with multiple instruments in the two sample setting (producing more conservative effect estimates). The section has now been rewritten to improve readability.

Top of page 6: Scale dependency- obviously scale dependent interactions may also be present in the absence of transformation (i.e. when the variables are on their original scale).

This is absolutely true. There is no guarantee that the interaction observed is not an artifact of the arbitrary means with which it is initially measured. This is often why the most convincing GxE interactions are identified using a combination of sources, and why surprising GxE interactions from large-scale interrogation of studies with genetic data warrant careful consideration before being used.

Line 159 Page 6: “Specifically, MR-GENIUS relies upon the residual error in a regression of the exposure upon the genetic IV to be heteroskedastic, such that (X|G) = E(ϵ2X|G) [13]. This is evaluated using a Breusch-Pagan test for heteroskedasticity [13].” Again, I find this a little unclear especially the equation. My understanding is that in a Breusch-Pagan test the square of the first stage residuals is regressed on the IV (here genotype). I assume you are saying that the expected value of the squared residuals given G doesn’t change across levels of G- but is this the correct way to say this using an equation?

This is an excellent point to make, as we have tried to present instrument strength with regard to MR-GENIUS as close to the representation in their paper as possible. We agree in hindsight that the equation is more confusing than helpful, and so have adjusted this section to highlight that the residual error from regressing X on G should be heteroskedastic, and that a BP-Test can be used to evaluate this criterion.

Line 163-171. I assume that these analyses to identify gene x covariate interactions are based on polygenic risk scores, rather than individual variants. I suggest you make this clear, especially since single variants are unlikely to show large GxE for most traits.

Thank you for the suggestion, and we completely agree. Often it is most appropriate to use a PRS as an instrument G, and we have now stated this explicitly in the paper.

Line 173 (page 6) “In previous work we show how assumption GxE2 is potentially violated when certain confounding structures exist, specifically, where Gi and Zki are simultaneously downstream of a confounder Ui or where there is an open path between the two variables through Ui”. Could the authors please discuss the intuition for this result perhaps with the aid of a path diagram (or description)? Do you need both paths to be present (i.e. a path from U to G AND a path from U to Z, or is only one sufficient to produce bias?). Why? Maybe this could be pointed out to the reader?

We agree that the description of interaction exogeneity needed greater attention, and thank the reviewer for highlighting this issue. We have edited this section to better reflect issues related to confounding in the context of MR-GxE. This was a particularly interesting comment, as it links to the use of MR-GxE as a test for instrument validity in a way not communicated in the paper. Many thanks.

“Relating assumption GxE2 to the MR-GENIUS approach, associations violating GxE2 would imply associations vary across values of the unmeasured confounders violating the second MR-GENIUS assumption [13]. However, this problem can be mitigated by incorporating additional interaction covariates within the MR-GENIUS model, as described in Tchetgen Tchetgen et al, 2021 [13].” Please expand on this.

We have now included a description of how such a model would be constructed, though expanded MR-GENIUS model does not receive much attention in the original cited paper.

“The third MR-GxE assumption requires pleiotropic effects of Gi upon Yi to remain constant across values of Zki, with the gene-by-covariate interaction being independent of Yi when conditioning on Xi.”. I would put in parentheses after this statement (i.e. beta_4 equals zero)

We agree this additional clarification approves readability, and have made the change to the paper.

Some further thoughts:

A table summarizing the core assumptions, the consequences for violating them, and methods of testing them might be useful for readers.

This is an excellent suggestion. We have now included an additional table (Table 6) in the discussion which provides a brief overview of the assumptions, their consequences, and methods for evaluating their plausibility.

One of the difficulties in ordinary MR studies is whether the SNPs chosen to be instruments are primarily associated with the exposure, or the outcome. Various methods (e.g. Steiger filtering) have been proposed to get a handle on this problem. My question is whether the addition of SNP*covariate terms in GZ where the SNP is primarily an outcome associated SNP is problematic? My guess would be yes (it implies beta_4 is not null). Do the authors think that this would be a common problem? If so, could they recommend procedures to guide against this possible source of bias?

This is a common problem, and really highlights the care that needs to be taken when selecting SNPs which form the instrument G. We would recommend researchers use tools like Steiger filtering, as well as interrogate individual SNPs which are used to see if there are likely pleiotropic pathways (using a tool such as PhenoScanner for example). We agree that if the SNPs are primarily related to the outcome there is a greater chance for violation of GxE3, though this would require the effect on the outcome to also vary across levels of the interaction covariate Z. In terms of a procedure, assuming SNPs are selected as in conventional MR (F<10, Steiger filtering, using LD pruning) the only additional steps we would suggest at present would be to perform the tests described in the paper. This is somewhat unsatisfying as a recommendation, but anything more would require greater knowledge of the mechanisms underpinning the adopted GxE interaction, and if that is available then it would be more valuable to examine mechanisms through direct effects on of G on Y could vary across Z from a biological perspective.

Sex seems an obvious candidate for a genotype x “environment” interaction. Was this included in the list of covariates and results not presented because the evidence for interaction was so low?

This is correct, and often sex can serve as an ideal interaction covariate, owing to its plausible independence from genetic instruments which do not lie on sex chromosomes. There were two reasons why these interactions were omitted:

1. The GWAS effect estimates used to create the PRS for BMI were adjusted for participant sex, and so we wanted to avoid issues of inadvertently conditioning on a variable multiple times.

2. It is plausible that the underlying causal effect of BMI on SBP differs by biological sex, and ideally we rely on the underlying effect of BMI on SBP being the same across interaction subgroups (assumption IV4 of a homogenous causal effect applies to MR-GxE just as it does in conventional MR). The extent to which this would be the case is unclear, and a thorough examination of this issue would unfortunately lie outside the scope of this paper.

Even so, biological sex has been used in interaction MR analyses, most notably evaluating the impact of alcohol on cardiovascular disease using historic sex-stratified alcohol consumption in East Asian countries. The results from such studies are certainly plausible, and a deeper analyses of this (and IV4 in general) is something that would be exciting to explore.

You have a “positive empirical control” in this paper (i.e. we know BMI causes SBP) would it be worth also including a negative empirical control? (i.e. two phenotypes where we are pretty sure the exposure doesn’t cause the outcome). I realize that this may potentially be a substantial amount of work so would be happy for the paper to be accepted without it, but i think it is something worth considering.

This is a fantastic idea, and really brings the method back to its roots in negative controls (we ideally want an interaction covariate subgroup for which the instrument and exposure are independent). Unfortunately, given the number of simulations and the applied analysis it would be difficult to include the work while keeping the paper at an accessible length, though it is something we are keen to explore in the near future.

Submitted filename: WS_Reviewer_response.docx

Click here for additional data file.

11 Jul 2022

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

PONE-D-22-06423R1

Dear Dr. Spiller ,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Momiao Xiong

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

You have address all the issues  which the reviewers are concerned.

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: (No Response)

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: (No Response)

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: (No Response)

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: (No Response)

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: (No Response)

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: (No Response)

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

**********

14 Jul 2022

PONE-D-22-06423R1

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

Dear Dr. Spiller:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Prof. Momiao Xiong

Academic Editor

PLOS ONE