Software application profile: mrrobust—a tool for performing two-sample summary Mendelian randomization analyses

Abstract Motivation In recent years, Mendelian randomization analysis using summary data from genome-wide association studies has become a popular approach for investigating causal relationships in epidemiology. The mrrobust Stata package implements several of the recently developed methods. Implementation mrrobust is freely available as a Stata package. General features The package includes inverse variance weighted estimation, as well as a range of median, modal and MR-Egger estimation methods. Using mrrobust, plots can be constructed visualizing each estimate either individually or simultaneously. The package also provides statistics such as IGX2, which are useful in assessing attenuation bias in causal estimates. Availability The software is freely available from GitHub [https://raw.github.com/remlapmot/mrrobust/master/].


Introduction
Mendelian randomization 1 has developed into a popular approach to examining causal relationships in epidemiology. 2,3 By employing genetic variants as instrumental variables (IVs) it is possible to limit bias from confounding, provided variants satisfy the assumptions of IV analysis. 1,4 For a genetic variant to serve as a suitable instrument, three assumptions must hold: (i) it must be associated with the exposure of interest; (ii) there must be no confounders of the instrument and outcome; and (iii) the instrument must not affect the outcome except via the exposure of interest. 5 Candidate variants are usually identified through large genome-wide association studies (GWASs). 6 However, IV analyses using single variants rarely have sufficient power to test hypotheses of interest. 6,7 One approach to increase the statistical power of Mendelian randomization studies is to use multiple genetic variants as instruments within a twosample summary framework. 8,9 Two-sample Mendelian randomization estimates the effect of the exposure using instrument-exposure and instrument-outcome associations from different samples, often through methods originally developed for meta-analysis. 8,9 This is particularly useful as MR estimators, such as MR Egger and median based regression, are robust to certain forms of violation of the third instrumental variable assumption. 8,10,11 Violations of this assumption can occur through directional pleiotropy, where a genetic variant affects the study outcome through pathways that are not mediated via the exposure. Such developments have contributed to the increasing popularity of twosample summary MR. 5 This paper introduces the mrrobust Stata package as a tool to help researchers implement two-sample MR analyses, and can be viewed as the Stata counterpart to toolkits such as the MR-Base web application, and the MendelianRandomization and TwoSampleMR R packages. 12,13 Whereas it is possible to conduct individual-level IV analyses in Stata using modules such as IVREG2, 14 twosample summary MR has previously required bespoke code to implement. The mrrobust package addresses this limitation, providing a suite of popular two-sample MR methods and sensitivity analyses. Before continuing, we briefly outline the three primary estimation methods included in the mrrobust package, using the notation of Bowden et al. 10,15 Methods

Inverse variance weighting (IVW)
To perform IVW, a weighted averageb IVW is calculated using the set of ratio estimatesb J for each individual variant J ¼ 1; 2; . . . ; j: 9 Ratio estimates are obtained for each variant by dividing the instrument-outcome association by the corresponding instrument-exposure association. Such association estimates are obtained by fitting simple linear regression models of the outcome and exposure upon the genetic variant, primarily by conducting a GWAS. Letĉ j and r 2 Yj denote the instrument-outcome association and variance, respectively, for the j th variant. The IVW estimate is then defined as: This corresponds to the estimate one would obtain from a weighted linear regression of the set of instrumentoutcome associations upon the set of instrument-exposure associations, constraining the intercept at the origin. 9 One drawback of the IVW approach is that causal effect estimates can be biased in cases where one or more variants exhibit directional pleiotropy. 9

MR-Egger regression
MR-Egger regression is valid under weaker assumptions than IVW, as it can provide unbiased causal effect estimates even if the variants have pleiotropic effects. In this case, the set of instrument-outcome associations is regressed upon the set of instrument-exposure associations, weighting the regression using precision of the instrumentoutcome associations, as in the IVW case. 8 However, MR-Egger does not constrain the intercept at the origin, and the intercept represents an estimate of the average directional pleiotropic effect across the set of variants. The slope of the model provides an unbiased estimate of the causal effect. 8,10 If there is little evidence of systematic differences between the IVW and MR-Egger, then the IVW should be preferred. The IVW is more efficient, but potentially less robust, and in such cases the IVW estimate is often the most appropriate estimate to adopt due to the greater precision of IVW estimates in comparison with other approaches. 10 If there are differences between the IVW and MR-Egger estimates, this may be due to pleiotropy or to heterogeneous treatment effects.
The utility of MR Egger regression hinges upon two core assumptions. First, the INstrument Strength Independent of Direct Effect (InSIDE) assumption requires the effects of single nucleotide polymorphisms (SNPs) on the exposure and their pleiotropic effects on the outcome to be independent. If the InSIDE assumption holds, estimates for variants with stronger instrument-exposure associations ðĉ j Þ will be closer to the true causal effect parameter than variants with weaker associations. 8 Second, the NO Measurement Error (NOME) assumption requires no measurement error to be present in the instrumentexposure associations, and therefore that the variance of the instrument-exposure association r 2 Xj ¼ 0. In cases where NOME is strictly satisfied, estimatesĉ j will be equal to c j and the variance of the ratio estimate for each var- . We further note that the NOME assumption applies to other two-sample MR approaches and is not therefore a unique feature of the MR Egger approach.
In cases where the NOME assumption is violated, individual variants will suffer from weak instrument bias, leading to attenuation of MR Egger estimates towards the null. This can occur if the SNPs were not genome-wide significant (p ¼ 5 Â 10 À8 ) or were selected from small GWAS. One novel approach to assessing the strength of the NOME assumption is to evaluate the I 2 GX statistic, interpreted as the relative degree of attenuation bias in the MR Egger regression in the interval (0, 1). 10 Thus, for example, an I 2 GX value of 0.7 represents an estimated relative bias of 30% towards the null. Further details regarding calculation of the I 2 GX statistic are presented in the Supplementary material, available at IJE online.

Weighted median
The weighted median approach is an adaptation of the simple median estimator for two-sample summary MR. 15 For a total number of variants J ¼ 2k þ 1, the simple median approach selects the middle ratio estimateb kþ1 , from ordered ratio estimatesb 1 ;b 2 ; . . .b j : 15 In cases where the total number of variants is even, the median is interpolated as 1 2b k þb kþ1 . As the simple median approach is inefficient, particularly in cases with variable precision in the set of ratio estimates, it is preferable to incorporate weights in a fashion similar to the IVW and MR Egger approaches. Let s j ¼ P j k¼1 w k be the sum of weights for the set of variants 1; 2; . . . j, standardized so the sum of weights s J ¼1. The weighted median estimator is the median of a distribution having estimateb j as its p j ¼ 100 s j À wj 2 À Á th percentile. 15 For the range of percentile values, we perform a linear extrapolation between neighbouring ratio estimates. An important assumption of the median summary MR approaches is that more than 50% of the genetic variants do not exhibit directional pleiotropy. In the simple median case, this threshold refers to the number of variants, whereas in the weighted median case, the 50% threshold is with respect to the weights of the non-pleiotropic variants. 15

Additional estimators
As two-sample MR represents a developing area of genetic epidemiology, novel approaches to causal effect estimation are incorporated into the mrrobust package through frequent updates. One such method is the mode-based estimator put forward by Hartwig et al. 16 Details on the implementation of this approach with accompanying examples can be found in the Supplementary material, available at IJE online.

Visualizing MR estimates
One useful approach to presenting the results of MR analyses is to produce a scatterplot, with the x and y axes representing the instrument-exposure and instrument-outcome associations, respectively, for each variant. If one were to draw a hypothetical regression line leading from the origin to each variant, the slope of the line would represent a ratio estimate of the causal effect using the single variant as an instrument, that is dividing the instrument-outcome association by the instrument-exposure association (defined as b j above). The precision of the instrument-outcome association estimate for each variant is illustrated using vertical error bars, whereas horizontal error bars pertaining to the instrument-exposure association may be omitted for clarity. As the IVW, MR-Egger, median and modal approaches essentially meta-analyse the set of ratio estimates, it is possible to include regression lines highlighting effect estimates of each approach for comparison. For such regression lines, positive and negative slopes are indicative of a positive or negative effect, respectively, whereas a slope of zero represents the absence of an observed association.

Implementation
The mrrobust package uses functions from moremata, 17 addplot 18 and the heterogi 19 command. For versions of Stata 13 and higher, it can be installed using the .net install command from [https://raw.github.com/remlapmot/mrro bust/master/]. For older versions of Stata, a zip archive of the files is freely available for download at: [https://github. com/remlapmot/mrrobust].
The package facilitates two-sample summary MR analyses with key features including: • IVW and MR-Egger regression approaches, including fixed effects MR-Egger regression, standard error correction and weighting options; • unweighted, weighted and penalized weighted median IV estimators, providing pleiotropy robust estimates in cases where fewer than 50% of the genetic instruments are valid; • modal estimation following Hartwig et al., 16 including weighted and unweighted variations; • presentation of heterogeneity statistics, statistics such as I 2 GX for use in assessing attenuation bias, 10 and Simulation Extrapolation (SIMEX) correction following Bowden et al; 10 • plotting tools to visualize IVW, MR-Egger and weighted median estimators, as well as density plotting with respect to implementing the modal estimator; • and illustrative examples and documentation using data from Do et al. 20 Applied examples: adiposity and height as predictors of serum glucose levels To illustrate key features of the mrrobust package, we perform two analyses investigating potential relationships between adiposity, height and serum glucose. Adiposity was selected owing to the vast body of evidence supporting a positive association with serum glucose levels, [21][22][23][24] whereas height was based upon limited evidence of association. [25][26][27] Glucose was selected as an outcome with respect to its hypothesized role in the development of type 2 diabetes. 21,27 Datasets were obtained from the MR-Base web application and pruned for linkage disequilibrium before conducting the analyses. 13 Applied example I: adiposity and serum glucose Though the relationship between adiposity and glucose has received much attention in the literature, such studies are predominantly observational and therefore may be subject to bias from confounding. This provides motivation for considering Mendelian randomization techniques which are able to control for such unobserved confounding. In the initial analysis, we select adiposity as an exposure measured using standardized body mass index (BMI), obtaining estimates of its associations with genotypes and their respective standard errors from Locke et al. 28 For the outcome, we consider log transformed measures of serum glucose logðmMÞ using effect estimates and standard errors from Shin et al. 29 The summary data used for this analysis are provided in the Supplementary material, available at IJE online. Adopting a GWAS significance Pvalue threshold of 5 Â 10 À8 ; a total of 79 independent SNPs were identified in both samples. We confirmed the linkage equilibrium (LD) between the SNPs using a clumping algorithm, a clumping distance of 10000 kb and an LD R 2 of 0.001. This resulted in a total of 79 SNPs for use as instrumental variables, details of which are presented in the Supplementary material, available at IJE online.
Using mrrobust, we conducted IVW, MR-Egger and weighted median regression approaches using the above summary data. The code for our analysis is in the Supplementary material, available at IJE online. For IVW and MR Egger approaches, the regression was weighted using the variance of the instrument-outcome association. The set of summary MR estimates are presented in Table 1A.
We find strong evidence of a positive association between BMI and serum glucose, using both IVW and weighted median methods. Considering the MR Egger case, a substantial average directional pleiotropic effect was not detected, and the lack of significance with respect to the effect estimate can be attributed to a lack of statistical power. An I 2 GX value of 0.88 was reported, which can be interpreted as a relative bias in the MR-Egger estimate of 12% towards the null. The estimates are shown in Figure 1A, constructed using the mreggerplot command which generates a scatterplot of the instrument-exposure and instrument-outcome associations for each variant. This shows the set of estimates to be in agreement, with the plot being constructed as previously described.
Applied example II: height and serum glucose As a further example, we consider the effect of standardized height (metres) upon serum glucose, using summary data from Wood et al. 30 and outcome summary data on log transformed serum glucose from Shin et al. 29 The summary data used for this analysis are provided in the Supplementary material, available at IJE online. We assess the SNPS for LD using criteria from the previous example and identify 367 SNPs as suitable instruments for the analysis, details of which are presented in the Supplementary material as above. The set of summary MR estimates are presented in Table 1B. From Table 1B, we find no evidence against the null hypothesis of no association between height and serum glucose levels using IVW, weighted median and MR Egger regression. Considering the MR Egger case, there appeared to be no evidence of directional pleiotropy, with an I 2 GX value of 0.90 indicating a relative bias of 10% towards the null. As in the previous example, a plot of the MR estimates can be generated using the mreggerplot command as shown in Figure 1B. In this scenario, the estimates appear in agreement, indicating a lack of evidence for a substantial directional pleiotropic effect.

Discussion
The mrrobust package is a freely available Stata package, containing a number of summary MR estimation methods which can be used to estimate causal effects. In the applied example, the mrrobust package was able to provide a series of estimates, finding evidence of a positive association between BMI and serum glucose and no evidence of association between height and serum glucose. One possible conclusion that can be drawn from these results is that previously reported associations between height and glucose are driven by confounding factors. 31,32 It is important, however, to consider the extent to which Mendelian randomization is appropriate for a given analysis and, by extension, situations in which mrrobust is suitable.
In the first instance, Mendelian randomization studies only produce unbiased estimates when genetic instruments satisfy the assumptions of each estimator (e.g. IVW, MR-Egger or weighted median). In two-sample analyses, genetic instruments should be associated with the exposure of interest at genome-wide levels of significance (satisfying the first instrumental variable assumption), and pruned for LD to limit the overlap between SNPs. The IVW estimator also requires that genetic variants should not have directional pleiotropic effects. The MR Egger and median estimators are robust to directional pleiotropy if the effects of the exposure are constant. MR Egger regression requires the InSIDE assumption, whereas median methods assume that the number of valid instruments is greater than 50%. For MR-Egger estimation where the value of I 2 GX is low, it is possible to use SIMEX to correct for regression attenuation towards the null. This is implemented using the mreggersimex command.
In this paper, we have presented the mrrobust Stata package as an accessible toolkit for performing summary MR and instrumental variable analysis using many instruments. It contains a range of summary MR approaches, and should make examining causal relationships using Mendelian randomization more accessible for genetic epidemiologists.

Supplementary Data
Within the Supplementary material, available at IJE online, we include example code and Stata output for each of the analyses performed within this paper, as well as the summary data obtained from the MRBase GWAS catalogue. We also include a brief summary of the I 2 GX statistic, as well as guidance on implementing and interpreting the modal estimator.