Multivariable Mendelian Randomization: The Use of Pleiotropic Genetic Variants to Estimate Causal Effects

A conventional Mendelian randomization analysis assesses the causal effect of a risk factor on an outcome by using genetic variants that are solely associated with the risk factor of interest as instrumental variables. However, in some cases, such as the case of triglyceride level as a risk factor for cardiovascular disease, it may be difficult to find a relevant genetic variant that is not also associated with related risk factors, such as other lipid fractions. Such a variant is known as pleiotropic. In this paper, we propose an extension of Mendelian randomization that uses multiple genetic variants associated with several measured risk factors to simultaneously estimate the causal effect of each of the risk factors on the outcome. This “multivariable Mendelian randomization” approach is similar to the simultaneous assessment of several treatments in a factorial randomized trial. In this paper, methods for estimating the causal effects are presented and compared using real and simulated data, and the assumptions necessary for a valid multivariable Mendelian randomization analysis are discussed. Subject to these assumptions, we demonstrate that triglyceride-related pathways have a causal effect on the risk of coronary heart disease independent of the effects of low-density lipoprotein cholesterol and high-density lipoprotein cholesterol.

Mendelian randomization employs genetic variants as instrumental variables to estimate the causal effect of a risk factor on an outcome using observational data, even in the presence of unmeasured confounding (1,2)

. A genetic variable is a valid instrumental variable if
• the variant is associated with the risk factor of interest, • the variant is not associated with any confounder of the risk factor-outcome association, and • the variant is conditionally independent of the outcome given the risk factor and confounders (3,4).
These assumptions can be illustrated using a causal directed acyclic graph, displaying a causal effect of one variable on another by an arrow and the absence of a direct causal effect by the lack of an arrow ( Figure 1) (5). Although a genetic variant need not be causally associated with the risk factor to be a valid instrumental variable, we assume that there exists a causal variant for which the measured variant is a proxy (6). In order to avoid violations of the second and third instrumental-variable assumptions, Mendelian randomization experiments have generally relied on genetic variants which are associated with a single risk factor. In practice, however, many variants are pleiotropic-that is, associated with multiple risk factors. Indeed, in some cases, there may be no variants which are solely associated with the risk factor of interest, and a Mendelian randomization analysis cannot be performed without considering pleiotropic variants. In any case, it may desirable to include information on pleiotropic variants in order to provide a more powerful analysis, provided that this does not prejudice its validity. It may also be that multiple quantitative traits relating to the same risk factor are of interest; for example, in cardiovascular disease, the concentration of lipoprotein(a) and the size of lipoprotein(a) particles (7). In this case, the relative proportions of risk reduction associated with interventions separately targeting lipoprotein(a) concentrations and the size of lipoprotein(a) particles may be of interest, and the traits may be regarded as independent risk factors, even if the same genetic variants influence both traits. The possibility of including multiple risk factors in an instrumental-variable analysis is discussed in many econometric textbooks (8), and applied instrumental-variable analyses involving multiple risk factors have been performed (9,10), but we are unaware of any application of the approach in genetic epidemiology.
The context of this paper is that there are measurements on multiple genetic variants and several associated risk factors, the causal effect of at least 1 of which on the outcome is of interest. We assume that the genetic variants do not influence the outcome via any pathway except those fully mediated by one of the measured risk factors or by some combination of the measured risk factors. Questions about variants with potentially unmeasured or unknown pleiotropic associations are reserved for the Discussion section. We initially discuss how pleiotropic associations may arise and the methods and assumptions necessary for estimating causal effects with several risk factors. We demonstrate the use of these methods in an applied example and then construct a simulation study with parameters chosen to be similar to those in the example to investigate how the methods perform. Finally, we discuss the application of the methods in epidemiologic practice and the interpretation of the applied example.

Mechanisms for association with multiple risk factors
There are several causal mechanisms by which a genetic variant may be associated with multiple risk factors (11).
We divide the possible mechanisms into 2 cases ( Figure 2): 1) vertical pleiotropy, where a variant is associated with multiple risk factors due to the causal effect of the primary risk factor on a secondary trait, and 2) functional pleiotropy, where the genetic variant is associated with multiple pathways. These 2 cases are not mutually exclusive; it is possible for both of them to exist for the same variant.
In the case of vertical pleiotropy, genetic variants associated with the primary risk factor would be expected to show consistent associations with the secondary trait. For example, genetic variants associated with higher body mass index (weight (kg)/height (m) 2 ) would be expected to show a consistent association with higher blood pressure. In this case, the causal effect of body mass index on the outcome would include an indirect effect, mediated through blood pressure, as well as a direct effect comprising all other pathways from body mass index to the outcome that are not operating via blood pressure. If it is assumed that the only pathway by which the genetic variants are associated with the secondary trait is via the primary risk factor, then a simple Mendelian randomization analysis would consistently estimate the causal effect of the primary risk factor on the outcome in spite of the apparent pleiotropic association.
In the case of functional pleiotropy, we suppose there are multiple genetic variants (at least as many variants as there are risk factors) which have different magnitudes of effect on the risk factors. These genetic variants can be used to estimate the causal effects of each risk factor even if none of the variants are specifically associated with any 1 particular risk factor. Since Mendelian randomization is analogous to a randomized trial (12), the use of genetic variants to assess the causal effects of multiple risk factors in a single study is analogous to a factorial randomized trial (see Web Figure 1, available at http://aje.oxfordjournals.org/), where multiple randomized interventions are simultaneously assessed (13). We refer to such an analysis as "multivariable Mendelian randomization."

Assumptions
For a multivariable Mendelian randomization analysis to be valid, the genetic variants must satisfy a similar set of assumptions as a conventional instrumental variable, but in this case the variants must be exclusively associated not with a single risk factor but with a set of measured risk factors. It is not necessary for each variant to be associated with every risk factor in the set, but a variant cannot have associations with the outcome except via the risk factors of interest. Specifically, for each variant, we assume that • the variant is associated with 1 or more of the risk factors, • the variant is not associated with a confounder of any of the risk factor-outcome associations, and • the variant is conditionally independent of the outcome given the risk factors and confounders.
In order to define and interpret causal effect estimates, we initially assume that the effect of each of the risk factors on the outcome is not mediated by another of the risk factors: We could intervene on each risk factor independently of all the  other risk factors, and an intervention on one risk factor will not influence any other risk factor. We refer to such risk factors as "causally independent." A causal directed acyclic graph corresponding to these assumptions with 3 genetic variants and 2 risk factors is presented in Figure 3A. We later relax the assumption of causal independence and allow causal effects between the risk factors, as in Figure 3B. In this paper, we assume that all associations are linear.

Individual-level data: 2-stage least squares method
If individual-level data are available on the genetic variants, the risk factors, and the outcome, causal effects of the risk factors on the outcome can be estimated using a 2-stage least squares (2SLS) approach (14). The risk factors are regressed on the genetic variants in a multivariate linear regression (first stage; a multivariate multiple regression, since there are multiple dependent variables and multiple independent variables), and then the outcome is regressed linearly on the fitted values of each of the risk factors (second stage; a univariate multiple regression, since there is 1 dependent variable and multiple independent variables). An alternative model for the genetic association with the risk factor could be proposed (such as one including interaction terms), but a model which is additive and linear in the variants is used here for comparability with the summarized data methods considered in the next section. Although a sequential regression approach gives the correct point estimates, the use of 2SLS software (such as the ivreg2 command in Stata (StataCorp LP, College Station, Texas) (15)) is recommended for estimation in practice to derive correct standard errors (16). Estimates from the method are valid even if the genetic variants are in linkage disequilibrium.

Summarized data: likelihood-based method
If individual-level data are not available but rather we have summarized (aggregated) data on the beta coefficients and standard errors for the associations between the genetic variants and the risk factors and outcome from separate univariate regressions, then the causal effects of the risk factors on the outcome can be estimated using a likelihood-based method (17). For example, if there are 2 risk factors X 1 and X 2 , each of which has no causal effect on the other, a multivariate normal distribution can be assumed for the beta coefficients representing the genetic associations with each of the risk factors X 1 and X 2 and the outcome Y from univariate linear regressions. Specifically, we assume that the estimate of association of genetic variant j, j = 1, ..., J, with X 1 is X 1j with standard error σ X1j , and similarly with X 2 (X 2j , standard error σ X2j ) and with Y (Y j , standard error σ Yj ): Estimates of the causal effects of X 1 and X 2 on Y (β 1 and β 2 ) can be obtained by numerical maximization of this likelihood function or by Bayesian methods (18). If there are K risk factors, data on K + 1 beta coefficients and corresponding standard errors are required for each genetic variant (X 1j , ..., X Kj , Y j ); there are K(J + 1) parameters in the model, and equation 1 shows a (K + 1)-variate normal distribution. If the outcome is binary and the beta coefficients for the genetic association with the outcome represent log relative risks or log odds ratios, then the causal effect estimates will represent log relative risks or log odds ratios, respectively. The model for the genetic associations with the outcome is linear in contributions from the genetic associations with the risk factors.
The parameters ρ 12 , ρ 1Y , and ρ 2Y represent the correlations between the beta coefficients. These will be nonzero if the beta coefficients are derived from the same data. Although these correlations can only be estimated from individual-level data, they should be approximately equal to the observational correlations between the variables X 1 , X 2 , and Y. It is advisable to conduct a sensitivity analysis to assess the impact of these parameter values on the causal estimates. If data on the associations with the risk factors and outcome are obtained from separate data sources, the relevant correlations will be zero.
Because the likelihood function comprises contributions from each variant, it is necessary that the information on the causal parameters provided by each variant be independent. Therefore, the genetic variants used in a summarized data analysis must be uncorrelated (not in linkage disequilibrium); otherwise confidence intervals estimated by the method will be too narrow (17). If the genetic variants are in linkage disequilibrium and the correlations between the variants are known, then these correlations can be used in a modified likelihood-based model: The correlations between the genetic variants are the same as the correlations between the beta coefficients corresponding to the genetic variants. If all of the variants are correlated, then instead of a separate (K +1)variate normal distribution for each of the J genetic variants, we can employ a J(K + 1)-variate normal distribution in equation 1 for all of the variants together.
Summarized data: regression-based method A further method which has been proposed for the analysis of summarized data is a linear regression-based approach, which gives estimates for each of the risk factors separately (19). This is performed in 2 stages. First, the beta coefficients for the genetic association with the outcome are regressed on the beta coefficients for the competing risk factors. Then the residuals from the first regression are regressed on the beta coefficients for the risk factor of interest.
For example, if there are 2 risk factors X 1 and X 2 and we want to estimate the effect of X 1 on Y, the 2 stages are: 2. Regress the residuals Y 1 ; Y 2 ; :::; Y J on the beta coefficients X 11 , X 12 , ..., X 1J . The regression-based estimate is the regression coefficient for X 1 .
The intuitive rationale is that these residuals represent any causal effects that are not explained by the alternative risk factors but are potentially explained by the risk factor of interest. However, it is an ad hoc approach which has no clear theoretical basis and which ignores the uncertainty in the beta coefficients (20). The causal nature of the associations of various lipid fractions, including low-density lipoprotein cholesterol (LDL-C), high-density lipoprotein cholesterol (HDL-C), and triglycerides, with the risk of coronary heart disease (CHD) is an issue with important consequences for disease prevention and drug development strategies. Observational studies have shown Figure 3. Causal directed acyclic graph illustrating multivariable Mendelian randomization in associations between variants G 1 , G 2 , and G 3 , risk factors X 1 and X 2 , and outcome Y. Confounders U 1 and U 2 are assumed to be unknown. A) Risk factors are causally independent (no causal effects between X 1 and X 2 ); B) risk factors are causally dependent (X 1 has a causal effect on X 2 ). associations of LDL-C with increased CHD risk, associations of HDL-C with decreased CHD risk, and a null association of triglycerides with CHD risk upon adjustment for a number of risk factors, including HDL-C and non-HDL-C concentrations, systolic blood pressure, and body mass index (21). However, a causal interpretation of these results may be misleading due to unmeasured confounding and the possibility that some of the covariates adjusted for lie on causal pathways, making their inclusion in a regression model inappropriate. Efforts to elucidate causal relationships using genetic variants in a Mendelian randomization approach have indicated that LDL-C plays a causal role in increasing the risk of CHD (22) and have suggested a null causal effect of HDL-C on CHD risk (23). However, the latter estimate had wide confidence intervals, because only a few variantsthose not associated with other lipid fractions-were included in the analysis. The inability to find variants associated with triglycerides and not associated with LDL-C or HDL-C has precluded reliable Mendelian randomization investigations for triglycerides.
Here we address the question of the causal effects of LDL-C, HDL-C, and triglycerides on CHD risk by multivariable Mendelian randomization using published data. Waterworth et al. (24) reported genetic associations from univariate regression analyses of 28 genetic variants with log-transformed LDL-C, HDL-C, and triglyceride concentrations and with the log odds of CHD. Details on the variants and the β coefficients for the associations are given in Web Table 1. We combine these beta coefficients in a multivariable Mendelian randomization analysis using the likelihoodbased and regression-based methods. Estimates using the likelihood-based method were obtained in a Bayesian framework using WinBUGS (http://www.mrc-bsu.cam.ac.uk/bugs). Technical details on the analysis and the software code used are provided in Web Appendix 1. A sensitivity analysis for the values of the correlation parameters (ρ 1Y , ρ 2Y , ...) is given in Web Table 2. Initially, we do not account for linkage disequilibrium between the genetic variants so that the analysis methods can be more directly compared.

SIMULATION STUDY
In order to assess the statistical properties of the analysis methods used, we perform a simulation study. The setup corresponds to the example above.
We generate data for 30,000 individuals indexed by i on 3 risk factors (X 1 , X 2 , X 3 ) and an outcome (Y) from the following data-generating model: g ij ∼ Binomialð2, 0:3Þ independently for each j ¼ 1; :::; 28: u 1i ; u 2i ; u 3i ∼ N ð0; 1Þ independently: ϵ X1i ; ϵ X2i ; ϵ X3i ; ϵ Yi ∼ N ð0; 1Þ independently: We set the genetic association parameters α G1j , α G2j , and α G3j for j = 1, ..., 28 to take the values shown in Web Table 1 for log-transformed LDL-C, HDL-C, and triglycerides, respectively, to be similar to the applied example. The instrumental variables g ij are drawn from binomial distributions, representing independent single-nucleotide polymorphisms with minor allele frequencies of 0.3. The causal effects β 1 , β 2 , and β 3 are set to 0.3, 0, and −0.1, respectively. The variables U 1 , U 2 , U 3 represent confounders, leading to correlations between X 1 , X 2 , X 3 and Y. The parameters β U1 , β U2 , β U3 are each fixed at 0.3 throughout, and the parameters α U1 , α U2 , α U3 are varied to take the value 0.3 or −0.3, leading to 8 different scenarios. The mean R 2 values, representing the proportion of variation in each risk factor explained by the 28 instrumental variables together, for X 1 , X 2 , and X 3 are 0.6%, 0.5%, and 3.2%, respectively, corresponding to mean F statistics of 6.6, 5.2, and 35.4. Estimates for 1,000 data sets generated in each of the 8 scenarios considered were derived using the 2SLS, likelihoodbased, and regression-based methods. The Monte Carlo standard errors for the mean estimates were approximately 0.003, and for the power they were approximately 1%. The likelihood-based method was applied in a Bayesian framework using WinBUGS; technical details on the analyses are provided in Web Appendix 2. Table 1 shows, for each scenario, the mean estimate, the mean standard error, the standard deviation of the estimates, and the statistical power to detect a nonzero effect at a nominal 5% significance level. For β 2 = 0, the expected power is 5%. We see that the mean estimates from the 2SLS and likelihood-based methods are close to the true values, with some deviation depending on the direction of confounding. This may represent the effect of weak instrument bias, corresponding to the low F statistics above for X 1 and X 2 (25). The efficiencies of the 2SLS and likelihood-based methods are similar, despite the reliance of the likelihood-based method on only summarized data.
In contrast, estimates from the regression-based method are biased, although they appear to give approximately valid inferences for the presence of a causal effect under the null. However, the power, especially the power to estimate β 3 , is much lower than that from the other methods. We therefore recommend the likelihood-based method for use in practice when summarized data are available.

Causal relationships between risk factors
In order to investigate the performance of the methods when there are causal relationships between the risk factors, we repeated the simulation but replaced the first line with The additional terms α X2 x 2i and α X3 x 3i represent causal effects of X 2 and X 3 (which were evaluated first) on X 1 . We set α U1 , α U2 , and α U3 equal to 0.3 (the first scenario considered above) and took 9 values of the parameters α X2 and α X3 ( Table 2). All other parameters were taken as in the original simulation study. Table 2 shows the mean estimates of the causal parameters derived from each of the methods. Aside from the  Abbreviations: SD, standard deviation; SE, standard error. a Three analytical methods (2-stage least squares, likelihood-based, and regression-based) were used to estimate the causal effects of X 1 on Y (β 1 = 0.3), X 2 on Y (β 2 = 0), and X 3 on Y (β 3 = -0.1). b Empirical power to detect a causal effect at a nominal 5% significance level.
regression-based method, which produces widely varying results, we see that the estimates do not change substantially as the parameters vary. This indicates that the 2SLS and likelihood-based methods estimate the direct causal effect of each risk factor on the outcome, not including paths operating via the other risk factors. This can lead to misleading conclusions about the total effects of the variables. For example, when α X2 = 0 and α X3 = 0.5, the total causal effect of X 3 on Y is β 3 þ α X1 β 1 ¼ À0:1 þ 0:5 × 0:3 ¼ 0:05 (including the path operating via X 1 ). The mean estimates from the 2SLS and likelihood-based methods are in the opposite direction of the true total effect. The differences between the estimated values of β 1 , β 2 , and β 3 in Table 2 and their true values can be attributed to weak instrument bias. Weak instrument scenarios will be common in multivariable Mendelian randomization, as it is necessary to use multiple instrumental variables to estimate the different causal effects. If the genetic associations with the risk factors and with the outcome are measured in the same data set, this will lead to bias in the direction of the observational association (25), whereas if the genetic associations with the risk factors and with the outcome come from different sources (known as a 2-sample instrumental-variable analysis), the bias will be in the direction of the null (26). To demonstrate this, we repeat the simulation study outlined in Web Appendix 3 with fewer nonweak instrumental variables. The mean estimates of the causal effect parameters are very close to their true values (Web Table 3). In this simulation context, we also explore the impact of interactions between genetic variants in their effects on the risk factors. The likelihoodbased method is robust to these misspecifications of the analysis model (Web Table 4). We also investigate a modification of the 2SLS method referred to as "sequential adjustment" by Holmes et al. (27), in which the causal effects of each of the risk factors are estimated in turn, and alternative risk factors are adjusted for as if they are confounders. Web Tables 5 and  6 indicate that substantial bias in the sequential adjustment method is evident even under the null, and its direction depends on the unknown confounders.
A nonzero causal estimate from a multivariable Mendelian randomization approach when there are causal relationships between the risk factors implies that the variable is an independent causal risk factor, in the sense that an intervention on the variable keeping the other risk factors constant (the controlled direct effect) would affect the outcome. However, the magnitude of the causal estimate may not represent the total causal effect of the variable on the outcome.

DISCUSSION
In this paper, we have introduced multivariable Mendelian randomization, an important and practically relevant extension of the Mendelian randomization paradigm for estimation of causal effects using genetic variants associated with more than 1 risk factor. For a valid analysis, the variants must satisfy a set of assumptions that are similar to those for an instrumental variable in conventional Mendelian randomization but are modified to take account of the multiple risk factors. A multivariable Mendelian randomization analysis may be beneficial where genetic variants are associated with several related risk factors, such as in the example with lipid fractions. It permits causal evaluation of a risk factor even if no variants are uniquely associated with it, as for triglycerides.
There are several limitations to this approach, many of which are shared with conventional Mendelian randomization (28,29). The specific association of a genetic variant with a single risk factor may be a reasonable assumption if the function of the genetic region where the variant is located is well-characterized. The assumption of an exclusive association between genetic variants and a set of risk factors is unlikely unless the risk factors have strong biological associations. However, if they are strongly associated, an assumption that the risk factors are causally independent is less plausible. Weak instrument bias, a phenomenon by which instrumental-variable estimates using variants not strongly associated with the risk factor of interest are biased, may be substantial if large numbers of genetic variants are used (30), as may be necessary in a multivariable Mendelian randomization a Three analytical methods (2-stage least squares, likelihood-based, and regression-based) were used to estimate direct causal effects of X 1 on Y (β 1 = 0.3), X 2 on Y (β 2 = 0), and X 3 on Y (β 3 = -0.1).
experiment. The effects of the risk factors on the outcome are assumed to be linear. While some researchers do not view this as a crucial assumption, citing the interpretation of an instrumental-variable estimate as an average causal effect (6,31), others have shown that departures from linearity can affect the findings of an instrumental-variable analysis (32). Our preference is to take a less literal view of causal estimates and to emphasize the outcome of a Mendelian randomization analysis as reflecting testing of a causal effect, rather than necessarily estimation of a causal parameter. From this perspective, while the linearity assumption is important, it is less important than the other instrumentalvariable assumptions.
Although multivariable Mendelian randomization is able to allow for genetic variants with "measured" pleiotropic associations, under the assumptions discussed in this paper, it is unable to deal with unmeasured or unknown pleiotropy. If an apparent causal finding is dependent on the association of a small number of variants with the outcome, then the result may plausibly be due to pleiotropic variants rather than being a true causal effect. However, if several variants in different genetic regions all demonstrate consistent associations with the outcome, then it is perhaps unlikely that all of these associations reflect pleiotropic mechanisms (33). In the case of our applied example, we constructed a lipid risk score for each variant by multiplying the genetic associations with each lipid fraction by the estimate of the lipid fraction's causal effect on CHD risk; details are given in Web Appendix 4. Web Figure 2 displays the lipid risk score plotted against the log odds ratio of CHD risk for each variant. Aside from variant rs2304130, it seems that the estimated causal effect of the lipid fractions on CHD risk is consistent across variants, and so unmeasured pleiotropy is unlikely to explain the causal effects.
We performed the applied analysis for the causal effects of lipid fractions on CHD risk in this paper at face value, assuming that the instrumental-variable assumptions were satisfied. In reality, the assumptions that there were only 3 lipid categories and that the effects of the genetic variants were restricted to these lipid fractions are oversimplifications. Some lipid fractions (e.g., intermediate-density lipoprotein cholesterol) were omitted from the analysis, and the variability of particle size within the categories was ignored (34). An assumption that the causal effect of triglycerides on CHD risk is independent of the effects of LDL-C and HDL-C may not be satisfied, particularly as evidenced by the attenuation of the observational association of triglycerides with CHD risk upon adjustment for HDL-C and non-HDL-C (21) and studies of the apolipoprotein A5 gene (APOA5) (35). Our simulations above have shown that in the case of causal effects between risk factors, estimates represent the direct causal effect of each risk factor on the outcome by a pathway that is not operating via the other risk factors. This may not equal the total causal effect of the risk factor, but it provides important evidence on the independent causal effect of the risk factor. Finally, our estimated causal odds ratio for a 30% decrease in LDL-C was surprisingly large in comparison with not only estimates of the effect of statin usage, which also reduces LDL-C levels by approximately 30% (22), but also a Mendelian randomization analysis that included variants solely associated with LDL-C (17). A nuanced interpretation of Mendelian randomization estimates, and of multivariable Mendelian randomization estimates in particular, is required in the light of the uncertainty of the underlying assumptions in any applied analysis. These aspects are considered in more detail elsewhere (20).
In conclusion, these findings provide some evidence of a causal effect of triglyceride-related pathways on CHD risk independent of the effects of LDL-C and HDL-C, but the weight of evidence attributed to the findings is a matter of interpretation depending on the degree of validity attributed to the instrumental-variable assumptions.
Note added in proof: While our manuscript was in press, we discovered a simple modification of the discussed regressionbased method that uses available statistical software to produce estimates with much better theoretical and statistical properties. See Burgess et al. (36) for further information.