A comparison of robust Mendelian randomization methods using summary data

1Mendelian randomization (MR) uses genetic variants as instrumental variables (IV) to determine whether an observational association between a modifiable exposure (often also called the intermediate variable under study or risk factor) and an outcome is consistent with a causal effect (Davey Smith & Ebrahim, 2003; Smith & Ebrahim, 2004). This approach is less vulnerable to traditional problems of epidemiological studies such as confounding and reverse causality. With the increasing availability of genome‐wide association studies that find robust associations between genetic variants and exposures of interest (Welter et al., 2014; Zheng et al., 2017), the potential of this approach is rapidly evolving. A genetic variant is a valid IV if (a) it is associated with the exposure, (b) it has no direct effect on the outcome, and (c) there are no associations between the variant and any potential confounders.

2There has been much discussion on the potentials and limitations of MR, as the IV assumptions cannot be fully tested (Davey Smith & Ebrahim, 2003; Glymour, Tchetgen Tchetgen, & Robins, 2012; VanderWeele, Tchetgen, Cornelis, & Kraft, 2014). Violation of the IV assumptions can lead to invalid conclusions in applied investigations. In practice, the exclusion restriction assumption that the proposed instruments (genetic variants) should not have a direct effect on the outcome of interest is debatable, particularly if the biological roles of the genetic variants are insufficiently understood (Glymour et al., 2012; von Hinke, Smith, Lawlor, Propper, & Windmeijer, 2016).

3Some genetic variants are associated with multiple traits (Sivakumaran et al., 2011; Solovieff, Cotsapas, Lee, Purcell, & Smoller, 2013). This is referred to as pleiotropy. There are two types of pleiotropy. Vertical pleiotropy occurs when a variant is directly associated with the exposure and another trait on the same biological pathway. This does not lead to violation of the IV assumptions provided the only causal pathway from the genetic variant to the outcome passes via the exposure. Horizontal pleiotropy occurs when the second trait is on a different biological pathway, and so there may exist different causal pathways from the variant to the outcome. This would violate the exclusion restriction assumption. To solve the problems that arise due to horizontal pleiotropy, several robust methods for MR have been developed that can provide reliable inferences when some genetic variants violate the IV assumptions, or when genetic variants violate the IV assumptions in a particular way. To our knowledge, a comprehensive review and simulation study to compare the statistical performance of these different methods has not been performed.

4To focus our simulation study and compare the most relevant robust methods for applied practice, we concentrate on methods that satisfy two criteria. First, the method requires only summary data on estimates (beta‐coefficients and standard errors) of genetic variant–exposure and genetic variant–outcome associations. We exclude methods that require individual participant data (Guo, Kang, TonyCai, & Small, 2018; Jiang et al., 2017; Kang, Zhang, Cai, & Small, 2016; Tchetgen Tchetgen, Sun, & Walter, 2017), and those that require data on additional variants not associated with the exposure (DiPrete, Burik, & Koellinger, 2018; O'Connor & Price, 2018). This is because the sharing of individual participant data is often impractical, so that many empirical researchers only have access to summary data, and for fairness, to ensure that all methods are using the same information to make inferences. Second, the method must be performed using standard statistical software packages. We exclude methods requiring convergence checks that cannot be easily automated for a simulation study (Berzuini, Guo, Burgess, & Bernardinelli, 2018) or are computationally infeasible for large numbers of variants in a reasonable running time (Burgess, Zuber, Gkatzionis, & Foley, 2018).

5In this article, we review nine robust methods for MR from a theoretical perspective, and evaluate their performance in a simulation study set in a two‐sample summary data setting. The methods differ in how they estimate a causal effect of the exposure on the outcome, as well as in the assumptions required for consistent estimation. We consider the weighted median, mode‐based estimation (MBE), MR‐Pleiotropy Residual Sum and Outlier (MR‐PRESSO), MR‐Robust, MR‐Lasso, MR‐Egger, contamination mixture, MR‐Mix, and MR‐RAPS methods. Some methods take a summarized measure of the variant‐specific causal estimates as the overall causal effect estimate (weighted median, and MBE), whereas others remove or downweight outliers (MR‐PRESSO, MR‐Lasso, and MR‐Robust), or attempt to model the distribution of the estimates from invalid IVs (MR‐Egger, contamination mixture, MR‐Mix, and MR‐RAPS). We also consider the performance of the methods in an empirical example to evaluate the causal effect of body mass index (BMI) on coronary artery disease risk.

6This paper is organized as follows. First, we give an overview of the robust methods and compare their theoretical properties. Then, we introduce the simulation framework and applied example to compare their properties in practice. Finally, we discuss the implications of this study for applied practice.

7We consider a model as previously described by Palmer, Thompson, Tobin, Sheehan, and Burton (2008) and Bowden et al. (2017) for J genetic variants G1,G2,…,GJ that are independent in their distributions, a modifiable exposure X, an outcome variable Y, and a confounder U. We assume that all relationships between variables are linear and homogeneous without effect modification, meaning that the same causal effect is estimated by any valid IV (Didelez & Sheehan, 2007). A visual representation of the model is shown in Figure 1.

8We assume that summary data are available on genetic associations with the exposure (beta‐coefficient βˆXj and standard error σXj) and with the outcome (beta‐coefficient βˆYj and standard error σYj) for each variant Gj.

9The causal effect of the exposure on the outcome can be estimated using a single genetic variant Gj by the ratio method The ratio estimate θˆRj is a consistent estimate of the causal effect if variant Gj satisfies the IV assumptions (Didelez & Sheehan, 2007). If the uncertainty in the genetic association with the exposure is low, then the standard error of the ratio estimate σRj can be approximated as (Thomas, Lawlor, & Thompson, 2007)

10The individual ratio estimates can be combined to obtain a single more efficient estimate. The optimally efficient combination of the ratio estimates is referred to as the inverse‐variance weighted (IVW) estimate (Burgess, Butterworth, & Thompson, 2013): The IVW estimate is equal to the estimate from the two‐stage least squares method that is performed using individual participant data (Burgess, Dudbridge, & Thompson, 2016). It is a weighted mean of the ratio estimates, where the weights are the inverse‐variances of the ratio estimates. The IVW estimate can also be obtained by weighted regression of the genetic associations with the outcome on the genetic associations with the exposure

11However, the IVW method has a 0% breakdown point, meaning that if only one genetic variant is not a valid IV, then the estimator is typically biased (Bowden, Davey Smith, Haycock, & Burgess, 2016). Bias will be present unless the pleiotropic effects of genetic variants average to zero (balanced pleiotropy) and the pleiotropic effects are independent of the genetic variant–exposure associations (see MR‐Egger method below; Bowden et al., 2017). With the increasing number of variants used in MR investigations, it is increasingly unlikely that all variants are valid IVs. Hence, it is crucial to consider robust estimation methods despite their lower statistical efficiency (i.e., lower power to detect a causal effect).

12We proceed to introduce the different robust methods we consider in this study in three categories: consensus methods, outlier‐robust methods, and modelling methods. A summary table comparing the methods is presented as Table 1.

13A consensus method is one that takes its causal estimate as a summary measure of the distribution of the ratio estimates. The most straightforward consensus method is the median method. Rather than taking a weighted mean of the ratio estimates as in the IVW method, we take the median of the ratio estimates. The median estimator is consistent (i.e., unbiased in large samples) even if up to 50% of the variants are invalid (Bowden et al., 2016). We consider a weighted version of the median method, where the median is taken from a distribution of the ratio estimates in which genetic variants with more precise ratio estimates receive more weight. Here, an unbiased estimate will be obtained if up to 50% of the weight comes from variants that are valid IVs. We refer to this as the “majority valid” assumption.

14A related assumption is the “plurality valid” assumption (Guo et al., 2018). In large samples, while ratio estimates for all valid IVs should equal the true causal effect, ratio estimates for invalid IVs will take different values. The “plurality valid” assumption is that, out of all the different values taken by ratio estimates in large samples (we term these the ratio estimands), the true causal effect is the value taken for the largest number of genetic variants (i.e., the modal ratio estimand). For example, the plurality assumption would be satisfied if only 40% of the genetic variants are valid instruments, provided that out of the remaining 60% invalid instruments, no larger group with the same ratio estimand exists. This assumption is also referred to as the Zero Modal Pleiotropy Assumption (ZEMPA; Hartwig, Davey Smith, & Bowden, 2017).

15This assumption is exploited by MBE method (Hartwig et al., 2017). As no two ratio estimates will be identical in finite samples, it is not possible to take the mode of the ratio estimates directly. In the MBE method, a normal density is drawn for each genetic variant centered at its ratio estimate. The spread of this density depends on a bandwidth parameter, and (for the weighted version of the MBE method) the precision of the ratio estimate. A smoothed density function is then constructed by summing these normal densities. The maximum of this distribution is the causal estimate.

16As these consensus methods take the median or mode of the ratio estimate distribution as the causal estimate, they are naturally robust to outliers, as the median and mode of a distribution are unaffected by the magnitude of extreme values. However, they are still influenced by outliers, as these variants still contribute to determining the location of the median or mode of a distribution. These methods can also be sensitive to changes in the ratio estimates for variants that contribute to the median or mode, and to the addition and removal of variants from the analysis. Additionally, the methods may not be as efficient as those that base their estimates on all the genetic variants.

17Next, we present three outlier‐robust methods. These methods either downweight or remove genetic variants from the analysis that have outlying ratio estimates. They provide consistent estimates under the same assumptions as the IVW method for the set of genetic variants that are not identified as outliers.

18In MR‐PRESSO method (Verbanck, Chen, Neale, & Do, 2018), the IVW method is implemented by regression using all the genetic variants, and the residual sum of squares (RSS) is calculated from the regression equation. The RSS is a heterogeneity measure for the ratio estimates. Then, the IVW method is performed omitting each genetic variant from the analysis in turn. If the RSS decreases substantially compared to a simulated expected distribution, then that variant is removed from the analysis. This procedure is repeated until no further variants are removed from the analysis. The causal estimate is then obtained by the IVW method using the remaining genetic variants.

19In MR‐Robust, the IVW method is performed by regression, except that instead of using ordinary least squares regression, MM‐estimation is used combined with Tukey's biweight loss function (Burgess, Bowden, Dudbridge, & Thompson, 2016). MM‐estimation provides robustness against influential points and Tukey's loss function provides robustness against outliers. Tukey's loss function is a truncated quadratic function, meaning that there is a limit in the degree to which an outlier contributes to the analysis (Mosteller & Tukey, 1977). This contrasts with the quadratic loss function used in ordinary least squares regression, which is unbounded, meaning that a single outlier can have an unlimited effect on the IVW estimate.

20In MR‐Lasso, the IVW regression model is augmented by adding an intercept term for each genetic variant (Burgess, Bowden, et al., 2016). The IVW estimate is the value of θ that minimizes In MR‐Lasso, we minimize where λ is a tuning parameter. As the regression equation contains more parameters than there are genetic variants, a lasso penalty term is added for identification (Windmeijer, Farbmacher, Davies, & DaveySmith, 2016). The intercept term θ0j represents the direct (pleiotropic) effect on the outcome, and should be zero for a valid IV, but will be non‐zero for an invalid IV. The causal estimate is then obtained by the IVW method using the genetic variants that had θ0j=0 in Equation (6). A heterogeneity criterion is used to determine the value of λ. Increasing λ means that more of the pleiotropy parameters equal zero and so the corresponding variants are included in the analysis; we increase λ step‐by‐step until one step before there is more heterogeneity in the ratio estimates for variants included in the analysis than expected by chance alone.

21The MR‐PRESSO and MR‐Lasso methods remove variants from the analysis, whereas MR‐Robust downweights variants. These methods will be valuable when there is a small number of genetic variants with heterogeneous ratio estimates, as they will be removed from the analysis or heavily downweighted, and so will not influence the overall estimate. In such a case, these methods are likely to be efficient, as they are based on the IVW method. The methods are less likely to be valuable when there is a larger number of genetic variants that are pleiotropic, particularly if the pleiotropic effects are small in magnitude, and when the average pleiotropic effect of non‐outliers is not zero.

22Finally, we present four methods that attempt to model the distribution of estimates from invalid IVs or make a specific assumption about the way in which the IV assumptions are violated. The MR‐Egger method is performed similarly to the IVW method, except that the regression model contains an intercept term θ0: This differs from the MR‐Lasso method, as there is only one intercept term, which represents the average pleiotropic effect. The MR‐Egger method gives consistent estimates of the causal effect under the Instrument Strength Independent of Direct Effect (InSIDE) assumption, which states that pleiotropic effects of genetic variants must be uncorrelated with genetic variant–exposure association. As the regression model is no longer symmetric to changes in the signs of the genetic association estimates (which result from switching the reference and effect alleles), we first reorientate the genetic associations before performing the regression by fixing all genetic associations with the exposure to be positive, and correspondingly changing the signs of the genetic associations with the outcome if necessary. The intercept in MR‐Egger also provides a test of the IV assumptions. The intercept will differ from zero when either the average pleiotropic effect is not zero, or the InSIDE assumption is violated. These two conditions (average pleiotropy of zero and InSIDE assumption satisfied) are precisely the conditions required for the IVW estimate to be unbiased.

23The contamination mixture method assumes that only some of the genetic variants are valid IVs (Burgess, Foley, Allara, Staley, & Howson, 2020). We construct a likelihood function from the ratio estimates. If a variant is a valid instrument, then its ratio estimate is assumed to be normally distributed about the true causal effect θ with variance σRj2. If a variant is not a valid instrument, then its ratio estimate is assumed to be normally distributed about zero with variance ψ2+σRj2, where ψ2 represents the variance of the estimands from invalid IVs. This parameter is specified by the analyst. We then maximize the likelihood over different values of the causal effect θ and different configurations of valid and invalid IVs. Maximization is performed in linear time by first constructing a profile likelihood as a function of θ, and then maximizing this function with respect to θ. The value of θ that maximizes the profile likelihood is the causal estimate.

24The MR‐Mix method (Qi & Chatterjee, 2020) is similar to the contamination mixture method, except that rather than dividing the genetic variants into valid and invalid IVs, the method divides variants into four categories: (a) variants that directly influence the exposure only (valid instruments), and (b) variants that influence the exposure and outcome, (c) that influence the outcome only, and (d) that neither influence the exposure or outcome (invalid instruments). This allows for more flexibility in modelling genetic variants, although potentially leads to more uncertainty in assigning genetic variants to categories.

25The MR‐Robust Adjusted Profile Score (RAPS; Zhao, Wang, Bowden, & Small, 2018) method models the pleiotropic effects of genetic variants directly using a random‐effects distribution. The pleiotropic effects are assumed to be normally distributed about zero with unknown variance. Estimates are obtained using a profile‐likelihood function for the causal effect and the variance of the pleiotropic effect distribution. To provide further robustness to outliers, either Tukey's biweight loss function or Huber's loss function (Mosteller & Tukey, 1977) can be used.

26Modelling methods are likely to be valuable when the modelling assumptions are correct, but not when the assumptions are incorrect. For example, the MR‐Egger method requires the InSIDE assumption to be satisfied to give a consistent estimate. The MR‐RAPS method is likely to perform well when pleiotropic effects truly are normally distributed about zero, but less well when they are not. The MR‐Mix method is likely to require large numbers of genetic variants to correct classify variants into the different categories. The contamination mixture method is less likely to be affected by modelling assumptions as it does not make such strict assumptions, but it is likely to be sensitive to specification of the variance parameter.

27To compare the performance of these methods in a realistic setting, we perform a simulation study. Full details of the simulation study are given in the Supporting Information Material.

28For each participant i, we simulate data on J genetic variants Gi1,Gi2,…,GiJ, a modifiable exposure Xi, an outcome variable Yi, and a confounder Ui (assumed unknown). The confounder is a linear function of the genetic variants and an independent error term εiU. The effect of variant j on the confounder is represented by coefficient ϕj (this is zero for a valid IV). The exposure is linear in the genetic variants, the confounder and an independent error term εiX. The effect of variant j on the exposure is represented by coefficient γj. The outcome is linear in the genetic variants, exposure, confounders, and an independent error term εiY. The effect of variant j on the outcome is represented by coefficient αj (again, this is zero for a valid IV). The effect of the exposure on the outcome is represented by θ. The genetic variants are modelled as single nucleotide polymorphisms (SNPs), with a varying minor allele frequency mafj, and take values 0, 1, or 2. The minor allele frequencies are drawn from an uniform distribution. The error terms εiU, εiX, and εiY each follow an independent normal distribution with mean 0 and unit variance.

29We can represent the model mathematically as

30In brief, we consider three scenarios: balanced pleiotropy, InSIDE satisfied—invalid IVs have direct effects on the outcome generated from a normal distribution centered at zero (for invalid instruments αj∼N(0,0.15), ϕj=0);directional pleiotropy, InSIDE satisfied—invalid IVs have direct effects on the outcome generated from a normal distribution centered away from zero (for invalid instruments αj∼N(0.1,0.075), ϕj=0);directional pleiotropy, InSIDE violated—invalid IVs have direct effects on the outcome generated from a normal distribution centered away from zero, and indirect effects on the outcome via the confounder (for invalid instruments αj∼N(0.1,0.075), ϕj∼U(0,0.1)).

31balanced pleiotropy, InSIDE satisfied—invalid IVs have direct effects on the outcome generated from a normal distribution centered at zero (for invalid instruments αj∼N(0,0.15), ϕj=0);

32directional pleiotropy, InSIDE satisfied—invalid IVs have direct effects on the outcome generated from a normal distribution centered away from zero (for invalid instruments αj∼N(0.1,0.075), ϕj=0);

33directional pleiotropy, InSIDE violated—invalid IVs have direct effects on the outcome generated from a normal distribution centered away from zero, and indirect effects on the outcome via the confounder (for invalid instruments αj∼N(0.1,0.075), ϕj∼U(0,0.1)).

34We simulated data on J=10, 30, and 100 genetic variants. A portion of the genetic variants were invalid IVs (30%, 50%, and 70%), and the direct effects of the variants explain 10% of the variance in the exposure. Summary genetic associations were calculated for the exposure and the outcome on nonoverlapping sets of individuals, each consisting of 10,000 individuals (Haycock et al., 2016). This situation is often referred to as two‐sample summary data MR (Pierce & Burgess, 2013). We considered situations with a null causal effect (θ=0) and a positive causal effect (θ=0.2). In total, 10,000 data sets were generated in each scenario.

35Methods can be compared by many metrics, including bias, empirical power, and standard deviation of estimates. We use mean squared error, which is the sum of bias squared plus variance, as the main criterion for comparing methods, as this provides a compromise between bias and precision. However, the relative importance of each metric will depend on the specific features of the application.

36We also compare the methods in an empirical example considering the effect of BMI on CAD risk. Since BMI is influenced by several biological mechanisms (Monnereau, Vogelezang, Kruithof, Jaddoe, & Felix, 2016), it is likely that the exclusion restriction is not satisfied for all associated genetic variants. Hence it is necessary to use robust methods to analyse these data. Additionally, we consider methods that detect outliers (MR‐Presso, MR‐Robust, MR‐Lasso, contamination mixture, MR‐Mix, and MR‐RAPS), and compare whether the same outliers are detected in each of these methods.

37We take 97 genome‐wide significant variants associated with BMI from the GIANT consortium (Locke et al., 2015). Associations with BMI are estimated in up to 339,224 participants from this consortium. Associations with coronary artery disease risk are estimated in up to 60,801 CAD cases and 123,504 controls from the CARDIoGRAMplusC4D Consortium (Nikpay et al., 2015). Association estimates for CAD were available for 94 of these variants.

38The scatter plot of the genetic associations with BMI and CAD risk is shown in Figure 2. While most variants seem to suggest a harmful effect of increased BMI on CAD risk, there is apparent heterogeneity in the IV estimates from each genetic variant individually, as evidenced by Cochran's Q test (Q‐statistic = 235.7, p <.001). Even after removing the five outliers as judged by the MR‐PRESSO method, which makes use of the heterogeneity statistic to identify outliers, we still reject the null hypothesis of that the regression model (including an intercept) fits the regression model with no additional variability than would be expected by chance (Q‐statistic = 125.9, p =.005). This suggests that some of the variants violate the IV assumptions.

39The results of the simulation study are presented in Table 2 (10 variants), Table 3 (30 variants), and Table 4 (100 variants). For each scenario, we present the mean, median, and standard deviation of estimates across simulations, and the empirical Type 1 error rate (for a null causal effect) or empirical power (for a positive causal effect) at a 95% confidence level. The empirical Type 1 error rate and empirical power are calculated as the proportion of simulated data sets where zero was not included in the 95% confidence interval. The mean squared error across simulations for the different methods with a null causal effect is presented in Figure 3 (Scenario 2), and Figure 4 (Scenario 3) for 30 variants. The corresponding plots for 10 variants (Figures S1 and S2) and 100 variants (Figures S3 and S4) were broadly similar.

40Overall, judging by mean squared error, the contamination mixture method performed best with 30% and 50% invalid variants. In some scenarios, other methods had lower mean squared error with 70% invalid variants. However, with some isolated exceptions, all the methods performed badly with 70% invalid instruments. Coverage for the contamination mixture method was around 10% or less when there were up to 50% invalid variants. This was also true for the MR‐Robust method, although that method had slightly lower power to detect a causal effect in some scenarios. Several other methods performed well in particular scenarios.

41Among consensus methods, estimates from the MBE method were less biased than those from the weighted median method, with lower Type 1 errors. The weighted median method had slightly higher power to detect a causal effect, although comparisons of power lose much of their value when a method has inflated Type 1 error rates. Performance of the MBE method improved as the number of variants increased. Among outlier‐robust methods, bias was greater for the MR‐Robust than the MR‐Lasso method. The MR‐Lasso method generally had the lower mean squared error when the invalidity was 50% or 70%, but MR‐Robust had the lower Type 1 error rates. Performance of the MR‐Robust method was better when there were at least 30 genetic variants. MR‐PRESSO had biased estimates with inflated Type 1 error rates even with 30% invalid variants, and performed particularly badly as the number of variants increased.

42The modelling methods performed well in some scenarios, but less well in others. This is unsurprising, as in some scenarios, consistency assumptions for the methods were satisfied, and in others they were not. The MR‐Egger method performed well in terms of Type 1 error rate in Scenarios 1 and 2, where the InSIDE assumption was satisfied. Estimates from the method were generally imprecise with low power. However, power in the MR‐Egger method depends on the genetic associations with the exposure varying substantially between variants, which was not the case in the simulation study (Burgess & Thompson, 2017). The contamination mixture method performed well with 30% and 50% valid instruments, with low bias and Type 1 error rates at or below 8% with 10 variants, 10% with 30 variants, and 11% with 100 variants. The MR‐Mix method performed badly throughout, with highly inflated Type 1 error rates in almost all scenarios with less than 100 instruments and comparatively low power to detect a causal effect. It performed slightly better with more genetic variants, although its performance was still worse than other methods. However, the method performed much better in a simulation comparison of methods performed by the authors of the MR‐Mix method (Qi & Chatterjee, 2019), in which the data‐generating model was more similar to the model assumed by the MR‐Mix method. The MR‐RAPS method performed well in Scenario 1, where its consistency assumption was satisfied, but less well in other scenarios with inflated Type 1 error rates. Its performance also worsened as more variants were included in the analysis.

43Results from the empirical example are shown in Table 5. All methods agree that there is a positive effect of BMI on CAD risk, except for the MR‐Mix method which gives a wide confidence interval that includes the null. The narrowest confidence intervals are for the outlier‐robust methods (MR‐Lasso, MR‐Robust, MR‐PRESSO), followed by the modelling methods except MR‐Mix and MR‐Egger (contamination mixture, MR‐RAPS), then the consensus methods (weighted median, MBE), and finally MR‐Egger and MR‐Mix.

44While the methods that detect outliers varied in terms of how lenient or strictly they identified outliers, they agreed on the order of outliers (Table S3). The MR‐Robust method was the most lenient, downweighting two variants as outliers. Each subsequent method in order of strictness identified all previously identified variants as outliers. MR‐PRESSO excluded the two variants identified by MR‐Robust plus an additional three variants. MR‐RAPS identified these five plus an additional two variants. MR‐Lasso identified an additional three variants, 10 in total. The contamination mixture method identified an additional 14 variants, 24 in total. MR‐Mix identified an additional 21 variants, 45 in total. This suggests that any difference between results from outlier‐robust methods are likely due to the strictness of outlier detection, rather than due to intrinsic differences in how the different methods select outliers. In several methods, the threshold at which outliers are detected can be varied by the analyst (e.g., by varying the penalization parameter λ in MR‐Lasso, or the significance threshold in MR‐PRESSO). In practice, rather than performing different outlier‐robust methods, it may be better to concentrate on one method, but vary this threshold. In our example, some of the variants that were the most pleiotropic in terms of their associations with other measured risk factors were only removed from the analysis by the MR‐Mix method (Table S3).

45In this paper, we have provided a review of robust methods for MR, focusing on methods that can be performed using summary data and implemented using standard statistical software. We have divided methods into three categories: consensus methods, outlier‐robust methods, and modelling methods. Methods were compared in three ways: by their theoretical properties, including the assumptions required for the method to give a consistent estimate, in an extensive simulation study, and in an empirical investigation.

46While the use of robust methods for MR analyses with multiple genetic variants is highly recommended, it is not practical or desirable to perform and report results from every single robust method that has been proposed. Guidance is therefore needed as to which robust methods should be performed in practice. As an example, if an investigator performed the MR‐PRESSO, MR‐Robust, and MR‐Lasso methods, they would have assessed robustness of the result to outliers, but they would not have not assessed other potential violations of the IV assumptions. The categorization of methods proposed here is not the only possible division of methods, but we hope it is practically useful. For instance, the contamination mixture and MR‐Mix methods make the same “plurality valid” assumption as the MBE method, and so could have been placed in the same category.

47The similarity and ubiquity of the “outlier‐robust” and “majority/plurality valid” assumptions should encourage investigators to consider methods that make alternative assumptions, such as the MR‐Egger method. While the InSIDE assumption is often not plausible (Burgess & Thompson, 2017), the MR‐Egger method and the intercept test have value in providing a different route to testing the validity of an MR study. Another potential choice is the constrained IV method, which uses information on measured confounders to construct a composite IV that is not associated with these confounders (Jiang et al., 2017). This method was not considered in the simulation study, as it requires additional data on confounders and individual participant data. Further methods development is needed to develop robust methods for summary data that make different consistency assumptions.

48We encourage researchers to perform robust methods from different categories, and that make varied consistency assumptions. For example, an investigator could perform the weighted median method (majority valid assumption), the contamination mixture method (plurality valid assumption), and the MR‐Egger method (InSIDE assumption). If there are a few clear outliers in the data, then an outlier‐robust method such as MR‐PRESSO (best used with few very distinct outliers) or MR‐Robust could also be performed. While we are hesitant to make a definitive recommendation as each method has its own strengths and weaknesses, this set of methods would be a reasonable compromise between performing too few methods and not adequately assessing the IV assumptions, and performing so many methods that clarity is obscured. Another danger of the use of large numbers of methods is the possibility to cherry‐pick results, either by an investigator seeking to present their results in a more positive light, or a reader picking the one method that gives a different result (such as the MR‐Mix method in our empirical example).

49One important limitation of these methods is the assumption that all valid IVs estimate the same causal effect. Particularly for complex exposures such as BMI, it is possible that different genetic variants have different ratio estimates not because they are invalid IVs, but because there are different ways of intervening on BMI that lead to different effects on the outcome. This can be remedied somewhat in methods based on the IVW method by using a random‐effects model (Bowden et al., 2017), or in the contamination mixture method, where causal effects evidenced by different sets of variants will lead to a multimodal likelihood function, and potentially a confidence interval that consists of more than one region.

50In summary, while robust methods for MR do not provide a perfect solution to violations of the IV assumptions, they are able to detect such violations and help investigators make more reliable causal inferences. Investigators should perform a range of robust methods that operate in different ways and make different assumptions to assess the robustness of findings from a MR investigation.