title: Pleiotropy robust methods for multivariable Mendelian randomization
authors: Andrew J. Grant Stephen Burgess
text: 
  	
  
  1]Andrew J. GrantCorresponding author. Email address: andrew.grant@mrc-bsu.cam.ac.uk
  1,2]Stephen Burgess
  [1]MRC Biostatistics Unit, University of Cambridge, Cambridge, UK
  [2]Cardiovascular Epidemiology Unit, University of Cambridge, Cambridge, UK
  
  Pleiotropy robust methods for multivariable Mendelian randomization
      [
      
  ===================================================================
  
  
  
  
  	Mendelian randomization is a powerful tool for inferring the presence, or otherwise, of causal effects from observational data. However, the nature of genetic variants is such that pleiotropy remains a barrier to valid causal effect estimation. There are many options in the literature for pleiotropy robust methods when studying the effects of a single risk factor on an outcome. However, there are few pleiotropy robust methods in the multivariable setting, that is, when there are multiple risk factors of interest. In this paper we introduce three methods which build on common approaches in the univariable setting: MVMR-Robust; MVMR-Median; and MVMR-Lasso. We discuss the properties of each of these methods and examine their performance in comparison to existing approaches in a simulation study. MVMR-Robust is shown to outperform existing outlier robust approaches when there are low levels of pleiotropy. MVMR-Lasso provides the best estimation in terms of mean squared error for moderate to high levels of pleiotropy, and can provide valid inference in a three sample setting. MVMR-Median performs well in terms of estimation across all scenarios considered, and provides valid inference up to a moderate level of pleiotropy. We demonstrate the methods in an applied example looking at the effects of intelligence, education and household income on the risk of Alzheimer's disease.
  
  
  
  
  § INTRODUCTION
  
  Mendelian randomization is a technique for estimating the causal effect of a risk factor on an outcome using observational data <cit.>. It uses genetic variants as instrumental variables and can provide valid causal effect estimation in the presence of unmeasured confounding. Three  assumptions are required in order that a genetic variant is a valid instrument: it must be associated with the risk factor of interest; it must not be associated with any confounder of the risk factor-outcome relationship; and it must be independent of the outcome conditional on the risk factor and confounders <cit.>.
  
  Genetic variants are good candidates for instrumental variables: they are naturally independent of many environmental factors which are common sources of confounding, and mitigate the potential for reverse causation. Furthermore, methods for Mendelian randomization have been developed which allow for combining many instruments in a single analysis, and which can also be used when only summary statistics of the associations between the genetic variants and traits are available <cit.>. These features allow practitioners to harness publicly available summary data from genome wide association studies (GWAS). Two-sample approaches, where the genetic variant-risk factor and genetic variant-outcome associations are estimated in different samples, open up vast combinations of risk factor-outcome relationships to be studied <cit.>. The major limitation in Mendelian randomization analyses is therefore the potential presence of pleiotropy, which is when genetic variants associate with traits other than the risk factor of interest. If any such trait provides an alternative causal pathway to the outcome not via the risk factor, then the corresponding genetic variants are invalid instruments and causal effect estimates may be biased.
  
  Multivariable Mendelian randomization fits multiple risk factors in a single model <cit.>. One motivation for its use is to account for pleiotropy in a univariable analysis via a set of measured covariates. It can be an important sensitivity analysis if there are known biological pathways linking the genetic variants and the outcome. Another motivation is if there are a number of correlated traits with shared genetic predictors which are all hypothesized to have potential causal effects on the outcome. A multivariable model can distinguish between the direct effects of the risk factors on the outcome and the total effects inclusive of mediators <cit.>. A genetic variant is a valid instrument for multivariable Mendelian randomization if: it is associated with at least one risk factor; it is independent of any confounder of each risk factor-outcome relationship; and it is independent of the outcome conditional on all risk factors and confounders. Causal pathways from a genetic variant to the outcome that do not pass via one or more of the risk factors are referred to as unmeasured pleiotropy (in contrast to measured pleiotropy, where such pathways are entirely account for via the set of risk factors). For the purposes of this paper, we use the word pleiotropy to mean unmeasured pleiotropy.
  
  There are a number methods in the literature for univariable Mendelian randomization (that is, when there is a single risk factor) which are robust to pleiotropy <cit.>. Each method provides valid estimation of the causal effect under different sets of assumptions. Although these assumptions are, generally, untestable, an applied analysis will typically employ a range of methods. Consistency of results across various methods which rely on different assumptions gives strength of evidence to the findings <cit.>. There are, however, few methods for pleiotropy robust multivariable Mendelian randomization. Valid estimation of causal effects, therefore, typically relies on the assumption that all causal pathways between the genetic variants and the outcome are accounted for via the measured risk factors.
  
  In this paper we propose a number of novel approaches to multivariable Mendelian randomization which provide robustness to different forms of pleiotropy. The methods are developed for use with summary level data, and so access to individual level data is not required. We examine the performance of the methods under various pleiotropic settings in a simulation study. We then demonstrate the methods in an applied analysis looking at the effects of intelligence, years of education and household income on the risk of Alzheimer's disease.
  
  
  
  § MODELLING ASSUMPTIONS
  
  
  
   §.§ Data generating model
  
  We assume the following model, which is similar to a multivariable version of the one set out by <cit.> in the single risk factor case. For individual i, let Y_i be the outcome, X_i1, …, X_iK be K risk factors, G_i1, …, G_ip be p genetic variants and U_i represent confounders of the risk factor-outcome relationships. The data generating model is:
  
      X_ik   = β_X0k + ∑_j=1^pβ_Xjk G_ij + γ_Xk U_i + v_Xik ,    k = 1, …, K 
      
      	Y_i   = θ_0 + ∑_k=1^Kθ_k X_ik + ∑_j=1^pα_j G_ij + γ_Y U_i + v_Yi,
  
  where v_Xik and v_Yi are independent error terms with mean zero. Note that the v_Xik are not necessarily independent of each other, and so the risk factors may be correlated via the correlation between these error terms as well as their common association with U_i. We assume that the genetic variants are independent of each other and independent of U_i. We further assume that p > K and the p × K matrix with (j, k)th element β_Xjk is of full column rank.
  
  
  
   §.§ Instrument validity and pleiotropy
  
  The relationships between a single genetic variant, the risk factors, confounders and outcome in model (<ref>)–(<ref>) are represented by the directed acyclic graph in Figure <ref>. For the jth genetic variant, pleiotropy is caused by the α_j term. Since the model allows for no direct association between the genetic variant and U_i, G_j is a valid instrument if at least one of β_Xj1, …, β_XjK are non-zero and α_j = 0. Note that although the model suggests α_j represents direct effects of the genetic variant on the outcome, it may also represent an association via an unmeasured trait.
  
  If α_j = 0 for all j, then there is no pleiotropy and all genetic variants are valid instruments. When not all α_j's are zero, we consider two patterns of pleiotropy. The first is referred to as balanced pleiotropy, which is where the α_j's are distributed with mean zero. The second is referred to as directional pleiotropy, which is where the α_j's are distributed with mean not equal to zero. In each case we assume that the α_j's are independent of each other. It may also be that most of the α_j are equal to zero but a relatively small number of them are non-zero and possibly large in magnitude. We will refer to these non-zero α_j's as outliers.
  
  
  
  
  
   §.§ Summary level data
  
  We denote by β̂_Xjk and β̂_Yj the estimates obtained by regressing the kth risk factor and outcome, respectively, on the jth genetic variant. We have that
  
      β̂_Xjk   = β_Xjk + ε_Xjk
      β̂_Yj   = α_j + ∑_k=1^Kβ_Xjkθ_k + ε_Yj ,
  
  where var( ε_Xjk) = σ_Xjk^2 and var( ε_Yj) = σ_Yj^2. In two sample Mendelian randomization the genetic variant-risk factor and genetic variant-outcome associations are estimated in separate samples and so ε_Xjk and ε_Yj are independent for all j. If β̂_Xjk and β̂_Xjl are obtained from separate samples, then ε_Xjk and ε_Xjl are independent. Otherwise, the correlation between ε_Xjk and ε_Xjl depends on the correlation between the kth and lth  risk factors. Finally, although in practice they are estimated from data, it commonly assumed that σ_Xjk^2 and σ_Yj^2 are known without error for all j,k.
  
  Although the model assumes that the risk factors and outcome are continuous, categorical traits are possible, and in fact common in practice. In this case, the relevant genetic variant-trait associations are estimated by logistic regression (or ordinal logistic regression, for ordinal variables with more than two categories) and represent the change in log odds ratio of the trait per extra effect allele in the genetic variant.
  
  
  
   §.§ Genetic variant orientation
  
  Each genetic variant can be coded in two ways, depending on which allele is chosen as the effect allele. The choice of effect allele is arbitrary, but will change the sign of the genetic variant-trait associations. Some Mendelian randomization methods may give different results depending on the orientation of the genetic variants. For example, in the single risk factor case, the inverse-variance weighted method <cit.> is not affected by genetic variant orientation, but methods which model pleiotropic effects, such as the MR-Egger method <cit.> and the lasso-based approach of <cit.>, are.
  
  In univariable Mendelian randomization, it is conventional to orientate the genetic variants such that an additional copy of the effect allele has a positive association with the risk factor. In the multiple risk factor case, however, this may be done in multiple ways, since forcing positive associations with respect to one risk factor may change the sign of the associations with respect to the others. <cit.> suggest orientating the genetic variants such that each one has a positive association with the primary risk factor of interest. If there is no single primary risk factor of interest, or as an additional sensitivity analysis, the impact of the orientation may be assessed by repeating the analysis multiple times, re-orientating the genetic variants with respect to each risk factor.
  
  
  
   §.§ The InSIDE assumption
  
  An assumption that is often required for pleiotropy robust Mendelian randomization is the InSIDE assumption (instrument strength independent of direct effects) <cit.>. There are two forms of the InSIDE assumption: a `population' level version is that, for all j, α_j is independent of each of β_Xj1, …, β_XjK; a `finite sample' version is that, for each j, the correlation between the sample estimates of α_j and each of β_Xj1, …, β_XjK, for the given set of data, is equal to zero. The latter will rarely be true in practice, since there will typically be residual correlation due to random variation. If the former is true, then this sample correlation will tend to 0 as the number of instruments increases.
  
  
  
  § METHODS
   
  We proceed to recall existing methods for multivariable Mendelian randomization in Sections <ref> and <ref> before introducing new approaches in Sections <ref>–<ref>.
  
  
  
   §.§ The inverse-variance weighted method
   
  The multivariable inverse-variance weighted (MVMR-IVW) method <cit.> fits the multiple linear regression model
  
      β̂_Yj = ∑_k=1^Kβ̂_Xjkθ_k + ε_j,
  
  j = 1, …, p, where ε_j is normally distributed with mean zero and variance σ_Yj^2. The estimator is obtained using weighted least squares estimation and is thus given by
  
      _θ_1, …, θ_K∑_j=1^p1/σ_Yj^2( β̂_Yj - ∑_k=1^Kβ̂_Xjkθ_k)^2.
  
  If all genetic variants are valid instruments, θ̂_IVW is a consistent estimator of θ. If not all genetic variants are valid instruments, the estimator remains consistent if pleiotropy is balanced and InSIDE is met. Thus, it is sensitive to outliers and directional pleiotropy.
  
  
  
   §.§ Existing pleiotropy robust methods
   
  If some of the genetic variants are invalid and pleiotropy is directional, the causal effect can still be consistently estimated using the MVMR-Egger method <cit.>. This method fits an intercept term in (<ref>) to account for pleiotropy. That is, we obtain the estimator from
  
      _θ_0, θ_1, …, θ_K∑_j=1^p1/σ_Yj^2( β̂_Yj - θ_0 - ∑_k=1^Kβ̂_Xjkθ_k)^2.
  
  Although the MVMR-Egger estimator is robust to invalidity of instruments, even when all instruments are invalid, it relies on the InSIDE assumption for consistent estimation. Furthermore, it results in lower precision. A final drawback is that it may produce different results depending on the orientation of the genetic variants.
  
  The MR-PRESSO method <cit.> has been proposed to handle the case where pleiotropy is balanced but there are outliers. Broadly speaking, the method performs a test based on a heterogeneity measure to identify outliers, which are then removed from the analysis. Although <cit.> describe the method for the single risk factor case, the authors have also produced a multivariable version, which is a straightforward extension. Specifically, the method computes inverse-variance weighted estimates by leaving out one genetic variant at a time. Letting θ̂_1,-j, …, θ̂_K,-j be the estimates obtained after leaving out the jth genetic variant, it then computes the following quantity, termed the global observed residual sum of squares:
  
      RSS_obs = ∑_j=1^p1/σ_Yj^2( β̂_Yj - ∑_k=1^Kβ̂_Xjkθ̂_k, -j)^2.
  
  This is compared with an expected residual sum of squares, which is computed multiple (M) times:
  
      RSS_exp^m = ∑_j=1^p1/σ_Yj^2( β̂_Yj^(m) - ∑_k=1^Kβ̂_Xjk^(m)θ̂_k, -j)^2 ,
  
  where β̂_Xjk^(m), j=1, …, p, k=1, …, K, are drawn from the normal distribution with mean β̂_Xjk and variance σ_Xjk^2, β̂_Yj^(m), j=1, …, p, are drawn the normal distribution with mean ∑_k=1^Kβ̂_Xjkθ̂_k, -j and variance σ_Yj^2, and m=1, …, M. Finally, for each j, an empirical p-value is computed as
  
      1/M∑_m=1^M1_> RSS_obs( RSS_exp^m),
  
  where 1_A( x ) is the indicator function. If the jth empirical p-value, multiplied by the number of variants (in order to apply a Bonferroni correction), is greater than the chosen significance level (for example, 0.05), then the respective genetic variant is identified as an outlier. If there are no outliers identified, the estimate obtained is the same as MVMR-IVW. If true outliers are identified and removed, it is expected to reduce the bias and be more efficient than MVMR-IVW. However, the method is not expected to perform well when there is directional pleiotropy, or there is a large number of invalid instruments.
  
  
  
   §.§ Robust regression
   
  
  A natural extension to MVMR-IVW is to use robust regression methods, for example MM-estimation. These methods provide robustness to observations which “contaminate” the data, such as outliers and influential observations (that is, those for which a small change in observed value results in a large change in parameter estimate). A method for performing robust regression in univariable Mendelian randomization is described in <cit.>, which uses MM-estimation along with Tukey's bisquare objective function. It is straightforward to extend this approach to the multivariable model: MM-estimation as described by <cit.> is done in a multivariable setting, and it can be implemented using existing software.
  
  This method of robust regression provides robustness to outliers by effectively capping residuals of a certain magnitude. The approach is thus expected to be robust to pleiotropy when there are a relatively small number of invalid instruments. In this case it should be unbiased and more efficient than MVMR-IVW. However, it may not perform well if there are a relatively large number of invalid instruments.
  
  
  
   §.§ Median based estimation
   
  An alternative approach to robust regression is to use least absolute deviations regression. That is, we estimate θ by
  
      _θ_1, …, θ_K∑_j=1^p1/σ_Yj^2| β̂_Yj - ∑_k=1^Kβ̂_Xjkθ_k|.
  
  Least absolute deviations regression is a special case of quantile regression which estimates the 50th percentile. Thus, (<ref>) is easily computed using techniques developed for quantile regression <cit.>. Since β̂_Yj and β̂_Xjk are continuous, (<ref>) has a unique solution with probability one.
  
  Similar to robust regression, least absolute deviations regression is less affected by outliers than least squares regression. It is not robust to influential observations, as robust regression is. However, it may be expected to perform better when the distribution of the β_Yj's are not symmetric. That is, it also provides robustness to directional pleiotropy. When K = 1, the estimator obtained using least absolute deviations regression is equivalent to the weighted median estimator for univariable Mendelian randomization proposed by <cit.> with weights given by |β̂_Xj1| / σ_Yj^2 (note that, strictly speaking, it is equivalent to the weighted empirical distribution method described in the supplementary material to that paper). The least absolute deviations regression approach can thus be thought of as a natural extension of median-based methods to the multivariable setting. We therefore refer to the method as MVMR-Median.
  
  A disadvantage of least absolute deviations regression is that we lose the asymptotic theory of least squares estimation which leads to easy to compute and accurate standard errors for use, for example, in inference. Confidence intervals are typically produced using a rank inversion technique, or via resampling methods (see, for example, <cit.>). Here we take advantage of the fact that we know the distribution of the genetic variant-trait associations, and implement a parametric bootstrap procedure, as follows. For each genetic variant, a genetic variant-outcome association is drawn from the normal distribution with mean β̂_Yj and variance σ_Yj^2, and genetic variant-risk factor associations are drawn from the multivariate normal distribution with mean ( β̂_Xj1, …, β̂_XjK) ' and covariance matrix diag( σ_Xj1^2, …, σ_XjK^2 ). The estimated standard error is the standard deviation of the estimates computed from multiple replications of this sampling. This approach does not take in to account correlation between the risk factors, however the simulation results presented in Section <ref> and the supplementary material show it still performs well in the correlated risk factor case.
  
  
  
   §.§ Regularization methods
   
  Under the assumption that some of the α_j's are zero and some are not, regularization methods for univariable Mendelian randomization have been proposed which include an intercept term for each genetic variant in the least squares equations (<ref>) and then apply lasso-type penalization to these terms. The penalization tends to shrink the intercept terms corresponding to the valid instruments toward zero. It thus accounts for the pleiotropy caused by invalid instruments, without the loss of power and need for the InSIDE assumption of Egger regression. The approach was first proposed by <cit.> in the individual level setting, and followed up by <cit.>. <cit.> developed a regularization approach using summary level data.
  
  In the multivariable setting we propose using
  
      _θ_01, …, θ_0p, θ_1, …, θ_K∑_j=1^p1/σ_Yj^2( β̂_Yj - θ_0j - ∑_k=1^Kβ̂_Xjkθ_k)^2 + λ∑_j=1^p| θ_0j|,
  
  for some tuning parameter λ > 0. This is not a standard lasso problem, since not all regression parameters are being penalized. However, the parameter estimates can be easily computed using the algorithm given in Section S.1 of the supplementary material, which uses only standard regression and lasso procedures. The tuning parameter controls the level of sparsity. The larger the value, the fewer genetic variants will be identified as invalid, and the estimate will approach the MVMR-IVW estimate. A data driven approach to choosing the tuning parameter is to use the heterogeneity stopping rule described by <cit.>.
  
  The lasso penalty will shrink some θ_0j's exactly to zero, thus identifying the corresponding genetic variants as being valid instruments. A post-lasso estimator takes the genetic variants identified as valid and fits a standard MVMR-IVW model using only these variants. Post-lasso estimators have been advocated by, for example, <cit.> and <cit.>, in order to avoid bias caused by the shrinkage of parameter estimates. The lasso algorithm is thus effectively used as a model selection technique.
  
  A limitation of regularization techniques generally is the inability to compute accurate standard errors. We can compute standard errors for the post-lasso estimator using a random effects model <cit.> in the post-selection regression. However, this ignores the uncertainty associated with the model selection event. As a result, the standard errors are likely to be too small, and the type I error rate inflated. We examine the effect of this in the simulation study presented in Section <ref>. A way around this is to use a three sample approach: here, a set of genetic variant-trait associations is used for performing the MVMR-Lasso procedure which are taken from a sample (or samples) which are independent of those from which the genetic variant-trait associations used for the post-lasso estimator are taken. In this way, the model selection and estimation procedures are independent and the correct type I error rate will be retained <cit.>. Although this restricts the potential scope for analyses, since multiple independent samples of genetic variant associations with the traits of interest are required, there are still a number of risk factor-outcome combinations which can be studied given the wide variety of GWAS results which are publicly available. Another promising development in the univariable setting is the use of a selective inference approach, which aims to derive a conditional distribution of the estimator given the model selection event <cit.>.
  
  One final point to note is that the solution to (<ref>) may be different depending on the orientation of the genetic variants. Following the convention used when performing MVMR-Egger, we propose orientating the genetic variants such that the genetic variant associations with the primary risk factor of interest are all positive.
  
  
  
  § SIMULATIONS
   
  We conducted a simulation study to compare the performance of the methods described in the previous section under scenarios with different amounts and types of pleiotropy. We simulated from model (<ref>)–(<ref>) with the intercepts set to zero, p=100 genetic variants, K=4 risk factors, n=20 000, γ_Xj = 1 / K, γ_Y = 1, β_Xjk∼Uniform( 0, 0.1 ), G_ij∼Binomial( 2, 0.3 ),
  
      U_i = ∑_j=1^pδ_j G_ij + w_i
  
  and v_Xi1, …, v_XiK, v_Yi, w_i∼ N ( 0, 1 ), independently. These parameter values give R^2 statistics (that is, the proportion of the variance in each risk factor explained by the genetic variants) of approximately 12%. Two sets of values for the causal effects were considered: in the first, θ_1=0.2, θ_2=0.1, θ_3=0.3, θ_4=0.4; in the second,  θ_1=0, θ_2=-0.1, θ_3=0.1, θ_4=0.2. Three scenarios were considered with different patterns of pleiotropy. For each scenario either 10%, 30% or 50% of genetic variants were invalid.
  
  	
    * Balanced pleiotropy and InSIDE assumption met: All δ_j's were set to zero, the proportion of α_j's set to zero was either 0.9, 0.7 or 0.5 and non-zero α_j's were generated from the N ( 0, 0.2^2) distribution.
  	
    * Directional pleiotropy and InSIDE assumption met: All δ_j's were set to zero, the proportion of α_j's set to zero was either 0.9, 0.7 or 0.5 and non-zero α_j's were generated from the N ( 0.1, 0.2^2) distribution.
  	
    * Directional pleiotropy and InSIDE assumption violated: All α_j's were set to zero, the proportion of δ_j's set to zero was either 0.9, 0.7 or 0.5 and non-zero δ_j's were generated from the Uniform( 0, 0.1 ) distribution (for detail on how the δ_j parameters cause the InSIDE assumption to be violated see, for example, <cit.>).
  
  
  For each scenario, level of pleiotropy, and set of θ_k values, the simulations were replicated 1 000 times. For each replication, the genetic variant-trait association estimates and their standard errors were computed from the individual level data using simple linear regression with an intercept. The causal effects were then estimated using the methods described in Section <ref>. The parameter of interest which we report on was the causal effect of the first risk factor on the outcome (that is, for the first set of θ_k values, there is a true causal effect, and for the second set of values there is no causal effect). The mean, standard deviation of estimates, mean standard error and power / type I error rate, at the 0.05 significance level, are shown in Tables <ref>–<ref>. The log of the mean squared errors across all scenarios are shown in Figure <ref>. Note that here the MVMR-Lasso method refers to the two sample post-lasso estimator (that is, with the estimate computed from the same samples that the instruments were selected in).
  
  All methods performed well in terms of bias when there was balanced pleiotropy. The MVMR-IVW and MVMR-Egger methods were biased when pleiotropy was directional, increasing as the proportion of pleiotropy increased. These methods were also less precise than all other methods, with the largest standard deviations of estimates, and were very low powered. In theory, MVMR-Egger should be robust to directional pleiotropy when InSIDE is met. However, there was a fair amount of bias in these scenarios. An explanation for this is that the bias is due to weak instruments, which this method is particularly susceptible to. This was further examined in the supplementary simulations (see the discussion at the end of this section).
  
  MVMR-Robust outperformed MVMR-PRESSO in all scenarios with lower bias, more precision and correct type I errors rates. MVMR-PRESSO had low bias at the lower level of pleiotropy, but did not perform well with moderate or high amounts. MVMR-Lasso was generally the most precise estimate: it had similar mean squared error to MVMR-Robust at 10% pleiotropy, but retained its performance in this regard at the higher levels of pleiotropy also. Similarly, it had comparable power to MVMR-Robust at 10% pleiotropy, but did much better as the proportion of pleiotropy increased. As expected, MVMR-Lasso had inflated type I error rates. MVMR-Median had comparable bias to MVMR-Lasso across all levels of pleiotropy. It was less precise and lower powered than MVMR-Lasso, but had type I error rates closer to the significance level. In terms of mean squared error, MVMR-Median was bettered uniformly across all scenarios only by MVMR-Lasso. As a further analysis, in Figure S1 in the supplementary material, we compare the mean squared errors with those from the estimates of the causal effect of the fourth risk factor, that is, where the true causal effect is θ_4 = 0.4. The average proportion of variation in the outcome explained by the fourth risk factor is approximately 19%, compared with approximately 11% for the first risk factor. However, with the exception of MVMR-Egger, the results are almost identical across all scenarios.
  
  The simulations were repeated for the cases where there were fewer instruments (p = 20, with the distribution of the β_Xjk parameters adjusted to retain similar R^2 values), where the risk factors were correlated by setting cor( v_Xik, v_Xil) = 0.5 for all k ≠ l, and where the genetic variant-trait associations were all estimated from the same sample (one sample Mendelian randomization). The results are shown in Tables <ref>–<ref> and Figures <ref>–<ref>. In each case, the results followed a similar pattern as before in terms of comparative performance of the different methods. In the fewer instruments case, there was slightly higher mean squared error across the board, which would be expected with fewer instruments, but the differences were not great. The results also support the assertion that the bias from MVMR-Egger shown in Tables <ref>–<ref> is due to weak instruments, since in this case the genetic variant-risk factor associations were adjusted to control the R^2 values, but the average F statistics were higher. The results from the correlated risk factor case suggested that the parametric bootstrap procedure (which effectively ignores risk factor correlation) still performs well.
  
  
  
  
  
  
  
  
  
  § APPLIED EXAMPLE: THE CAUSAL EFFECT OF INTELLIGENCE, EDUCATION AND HOUSEHOLD INCOME ON ALZHEIMER'S DISEASE
  
  In this section we consider an applied example looking at the causal effects of intelligence, years of education and household income on Alzheimer's disease. The effects of intelligence and years of education on health outcomes have been studied by <cit.> and <cit.>. A multivariable approach is important in this case since intelligence and years of education are highly correlated. <cit.> used univariable Mendelian randomization with intelligence and years of education, separately, as risk factors, and Alzheimer's disease as the outcome. The results suggest that both risk factors have a protective effect on Alzheimer's disease. However, when both risk factors are included in a multivariable model, using MVMR-IVW, the effect of years of education, independent of intelligence, shifts toward the null. The implication is that years of education only has a causal effect on the odds ratio of Alzheimer's disease via its effect on intelligence. Here, we reconsider this example and include household income as an extra risk factor.
  
  Genetic variant associations with intelligence and years of education are taken from the GWAS of <cit.> and <cit.>, respectively. By clumping the combined list of genetic variants which associate with each risk factor at the genome wide significance level, <cit.> arrived at a list of 219 independent genetic variants to be used as instruments in multivariable Mendelian randomization analyses. We obtained the associations between these genetic variants and household income from the UK Biobank (sourced from <http://www.nealelab.is/uk-biobank/>). Note that household income is an ordinal categorical variable, and so the genetic variant associations represent the increase in log odds of being in a higher income category per extra effect allele. Genetic variant associations with Alzheimer's disease were obtained from the GWAS of <cit.>. In total, 213 of the genetic variants used by <cit.> were available in both of the household income and Alzheimer's disease datasets, and we used these as instruments in our analysis. Note that the genetic variant associations with both intelligence and years of education were all in the same direction, and so they were orientated in our analysis to be all positive with respect to these traits.
  
  Figure <ref> shows a plot of the residuals vs fitted values after fitting the MVMR-IVW model to the data. The vertical error bars indicate ±σ_Yj for each genetic variant. The plot provides a way of visualising heterogeneity in the multivariable setting, similar to the scatterplots of β̂_Xj against β̂_Yj commonly used in the univariable case. Although there is little evidence of directional pleiotropy, there may be some outliers. Figure <ref> shows scatterplots of each pair of genetic variant-risk factor associations. There appears to be reasonably strong correlation between the genetic variant associations with years of education and household income, and low to moderate correlation between the other two pairs of associations.
  
  We used the univariable inverse-variance weighted method to estimate the causal effect of each risk factor separately on the odds ratio of Alzheimer's disease. We then applied each of the multivariable methods discussed in Section <ref> with all three risk factors included. The results are shown in Figure <ref>(a)–(c).
  
  The univariable analyses suggest that intelligence and years of education both have a protective effect on Alzheimer's disease, in line with the results of <cit.>. The estimated log causal odds ratio of Alzheimer's disease per one standard deviation increase in intelligence is -4.20 (95% CI -0.57 to -0.27), and per one standard deviation increase in years of education is -0.59 (95% CI -0.83 to -0.36). The estimated log causal odds ratio of Alzheimer's disease per unit increase in log odds ratio for a higher household income bracket is -0.60 (95% CI -0.89 to -0.31). Using the MVMR-IVW model, the estimates of the log causal odds ratio from both years of education and household income attenuated to the null, with 95% confidence intervals overlapping zero. The multivariable model however still suggests a protective effect from intelligence, with an estimated odds ratio of -0.47 (95% CI -0.86 to -0.07).
  
  The pleiotropy robust multivariable methods gave results which were broadly consistent with the MVMR-IVW results (see Table <ref>). The MVMR-Egger method suggested a null causal effect of intelligence on Alzheimer's disease, however the point estimate is still in the same direction and all other methods were in line with MVMR-IVW. Note that the MR-PRESSO outlier test did not detect any outliers, but the MVMR-Lasso method identified 15 genetic variants as pleiotropic, which were removed before computing the post-lasso estimator. These genetic variants are indicated in Figure <ref>. Interestingly, the post-lasso estimate for the causal effect of years of education on Alzheimer's disease was positive, whereas all other estimates of this effect were negative. However, the confidence interval still included the null.
  
  The consistency of the findings give strength to the assertion that intelligence has a causally protective effect on Alzheimer's disease, independent of years of education and household income. There are two potential explanations for the findings relating to years of education and household income. One is that these risk factors affect Alzheimer's disease via their association with intelligence. That is, intelligence is a mediator of the effect of years of education and/or household income. The other is that these risk factors sit on pleiotropic pathways between the genetic variants and Alzheimer's disease, not necessarily passing through intelligence. It is also possible that both explanations are true (that is, that there is both pleiotropic effects and mediation via intelligence).
  
  
  
  
  
  
  
  
  
  
  
  § DISCUSSION
  
  In this paper we have presented methods for performing multivariable Mendelian randomization which are robust to pleiotropy. Existing methods either allow for invalidity at the cost of low precision and the InSIDE assumption (MVMR-Egger), or were developed for the case where pleiotropy is balanced and there are a relatively small number of outliers (MVMR-PRESSO). We have considered methods which can handle higher proportions of invalidity and directional pleiotropy.
  
  When there is evidence of relatively few invalid instruments, MVMR-Robust was shown to outperform MVMR-PRESSO in all scenarios considered. MVMR-Lasso, another method which aims to identify and downweight outliers, performed best overall in terms of mean squared error, even when half of the genetic variants were invalid instruments and pleiotropy was directional. Although type I error rates were inflated, this can be mitigated when a three sample approach is possible. MVMR-Median was shown to perform almost as well as MVMR-Lasso in terms of mean squared error, and retained correct type I error rates at higher levels of pleiotropy. As demonstrated in the applied example, a plot of the residuals versus fitted values from the MVMR-IVW method can be used to visualise potential outliers and pleiotropy, and to help determine the most appropriate choice of robust method.
  
  The work has some limitations in the modelling assumptions made, in particular of the linearity and homogeneity (that is, no effect modification) of the effects of the risk factors on the outcome. Furthermore, although we can handle non-continuous traits via the use of logistic regression to produce summary statistics, this may cause bias in the causal effect estimates due to the non-collapsibility of the odds ratio. Nonetheless, violations of these assumptions tend to attenuate causal effect estimates toward the null <cit.>.
  
  Another limitation is the assumption of no measurement error of the genetic variant-risk factor associations, equivalent to assuming σ_Xjk = 0 for all j, k. This is a common assumption in Mendelian randomization analyses, and is justified by the very large sample sizes that these associations are typically estimated in, in contrast to the genetic variant-outcome associations which may be estimated using a relatively small number of cases vs controls. Provided the genetic variants strongly predict each risk factor, conditional on all the other risk factors, this assumption will have little influence on the analysis. Otherwise, the results may be subject to weak instrument bias. In practice, weak instrument bias is mitigated by selecting a set of genetic variants which associate with the risk factors according to some threshold related to, for example, an F-statistic or p-value. <cit.> provide a test to help to diagnose (conditionally) weak instruments and propose an approach to adjust standard multivariable Mendelian randomization estimators to account for the σ_Xjk's being non-zero. Although in univariable Mendelian randomization measurement error will bias causal effect estimates toward the null <cit.>, this will not necessarily be the case in the multivariable setting. Assessing the impact of measurement error in the multiple risk factor case, and how to account for this, is an active area of research.
  
  In summary, the methods we have presented provide new ways for performing Mendelian randomization with multiple risk factors which are robust to different forms of pleiotropy. Each has advantages when applied to specific scenarios. Together with MVMR-Egger, these methods provide a suite of sensitivity analyses for multivariable Mendelian randomization.
  
  
  
  § SOFTWARE
  
  R code for performing the methods described in this paper, and for reproducing the simulation results, can be found at <https://github.com/aj-grant/robust-mvmr>. Existing R packages that are used to implement the various methods include: MendelianRandomization <cit.>; MR-PRESSO <cit.>; robustbase <cit.>; quantreg <cit.>; and glmnet <cit.>.
  
  
  
  § ACKNOWLEDGMENTS
  
  The authors thank Neil Davies for supplying some of the data used in the applied example, and Max Mandl for supplying a correction to the multivariable version of the mr_presso R function.
  
  
  
  § FUNDING
  
  Andrew J. Grant and Stephen Burgess are supported by a Sir Henry Dale Fellowship jointly funded by the Wellcome Trust and the Royal Society (grant number 204623/Z/16/Z). This research was funded by the NIHR Cambridge Biomedical Research Centre. The views expressed are those of the authors and not necessarily those of the NHS, the NIHR or the Department of Health and Social Care.
  
  
  unsrtnat_ag
  
  	
  
  
  
  § SUPPLEMENTARY MATERIAL
  
  
  
  
  
  
  
  
  
  
  
   §.§ Algorithm for performing the regularization approach
   
  Let β̂_X be the p × K matrix with (j, k)th element β̂_Xjk, β̂_Y be the vector of length p with jth element β̂_Yj, S be the p × p diagonal matrix with (j, j)th element σ_Yj^-2 and θ_0 be the vector of length p with jth element θ_0j. We denote by P_β̂_X = S^1/2β̂_X( β̂_X'Sβ̂_X)^-1β̂_X'S^1/2 the projection onto the column space of S^1/2β̂_X and I_p the identity matrix of dimension p. We can solve (<ref>), for a given value of λ, using the following procedure.
  	
  		
    * Let
  		
      θ̂_0λ = _θ_0‖( I_p - P_β̂_X) S^1/2( β̂_Y - θ_0) ‖^2 + λ∑_j=1^p| θ_0j| ,
  
  		where ‖·‖ denotes the ℓ_2 norm.
  		
    * let
  		
      θ̂_λ = ( β̂_X'Sβ̂_X)^-1( β̂_X'Sβ̂_Y) .
  
  	
  The kth element of θ̂_λ is the estimate of θ_k for given λ. Note that Step 1 is now a standard lasso, with responses ( I_p - P_β̂_X) S^1/2β̂_Y and design matrix ( I_p - P_β̂_X) S^1/2, and can be computed with standard software.
  
  
  
   §.§ Supplementary simulation results
   
  In this section we present the results of the supplementary simulation studies described in Section <ref>. Tables <ref>–<ref> and Figure <ref> show the results for the case where p=20 and β_Xjk∼Uniform( 0, 0.22 ). Tables <ref>–<ref> and Figure <ref> show the results for the case where the risk factors are correlated with cor( v_Xik, v_Xil) = 0.5 for all k ≠ l. Tables <ref>–<ref> and Figure <ref> show the results for the case where the genetic variant-trait associations were all estimated from the same sample (one sample Mendelian randomization).
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  