Robust inference in summary data Mendelian randomization via the zero modal pleiotropy assumption

Abstract Background Mendelian randomization (MR) is being increasingly used to strengthen causal inference in observational studies. Availability of summary data of genetic associations for a variety of phenotypes from large genome-wide association studies (GWAS) allows straightforward application of MR using summary data methods, typically in a two-sample design. In addition to the conventional inverse variance weighting (IVW) method, recently developed summary data MR methods, such as the MR-Egger and weighted median approaches, allow a relaxation of the instrumental variable assumptions. Methods Here, a new method - the mode-based estimate (MBE) - is proposed to obtain a single causal effect estimate from multiple genetic instruments. The MBE is consistent when the largest number of similar (identical in infinite samples) individual-instrument causal effect estimates comes from valid instruments, even if the majority of instruments are invalid. We evaluate the performance of the method in simulations designed to mimic the two-sample summary data setting, and demonstrate its use by investigating the causal effect of plasma lipid fractions and urate levels on coronary heart disease risk. Results The MBE presented less bias and lower type-I error rates than other methods under the null in many situations. Its power to detect a causal effect was smaller compared with the IVW and weighted median methods, but was larger than that of MR-Egger regression, with sample size requirements typically smaller than those available from GWAS consortia. Conclusions The MBE relaxes the instrumental variable assumptions, and should be used in combination with other approaches in sensitivity analyses.


Introduction
Using germline genetic variants as instrumental variables of modifiable exposure phenotypes can strengthen causal inference in observational studies by applying the principles of Mendelian randomization (MR). 1,2 This method has already been used to address causality in several exposure-outcome combinations and has become a common feature in the recent epidemiological literature. 3 Causal inference using MR relies on the instrumental variable assumptions, which require that the genetic variant is: (i) associated with the exposure; (ii) independent of confounders of the exposure-outcome association; and (iii) independent of the outcome after conditioning on the exposure and all exposure-outcome confounders.
Recent MR methods allow performing MR with multiple genetic instruments, typically single nucleotide polymorphisms (SNPs), using summary data estimates from genome-wide association studies (GWAS). 4 Given the increasing number of publicly available summary statistics from large GWAS consortia, summary data MR methods enable many causal hypotheses to be rapidly interrogated without the administrative burden and cooperation required to perform equivalent individual-level data analyses. 5,6 However, using many instruments in an MR analysis increases the probability of including at least one invalid instrument, which could easily bias the estimate. For example, the inverse variance weighting (IVW) method requires that either all variants are valid instruments or that there is balanced horizontal pleiotropy (i.e. horizontal pleiotropic effects of individual instruments sum to zero) and that such pleiotropic effects are independent of instrument strength across all variants (i.e. the Instrument Strength Independent of Direct Effects -InSIDE -assumption). 4,7 More recently, other summary data MR methods that allow relaxion (but not elimination) of the instrumental variable assumptions regarding horizontal pleiotropy have been proposed. 8,9 In this paper, we describe a new summary data MR method -the mode-based estimate (MBE). We clarify when this will be a consistent estimate of the causal effect, compare it with established summary data MR methods using simulations and illustrate its application using real data examples.

Methods
In order to motivate the summary data methods discussed in this paper, we assume the following data-generating model linking genetic variant G j (j ¼ 1; . . . ; L), a continuous exposure X and outcome Y for subject i: Here, b Xj and b Yj ¼ bb Xj þ a j represent G j 's true association with the exposure and outcome, respectively. bb Xj is the effect of G j on Y through X, where b is the causal effect of X on Y we wish to estimate. The term a j represents the association between G j and Y not through the exposure of interest, due to horizontal pleiotropy. The error terms k Xij and k Yij will generally be correlated when collected on the same individuals. However, we will mainly focus on the two-sample setting where the error terms are independent, because independent samples are used to fit models (1) and (2). For simplicity, we will also assume that all L genetic variants are mutually independent of one another.
Letb Xj andb Yj represent the SNP-exposure and SNPoutcome association estimates for variant j, respectively, and let r 2 Xj and r 2 Yj represent the variance ofb Xj andb Yj , respectively. The ratio estimate 10,11 for the causal effect b using variant j alone is equal to: the standard error of which (r Rj ) can be obtained using the delta method 12 as follows: The standard error in (4) can be simplified to r Y j =jb Xj j when the variance of the SNP-exposure association r 2 Xj is small enough to be considered 'ignorable', or equivalently thatb Xj ¼ b Xj . This is referred to as the NO Measurement Error (NOME) assumption. 13 The ratio estimateb Rj is a crude measure of causal effect, but has a major advantage over more sophisticated methods in that it can be calculated using summary data estimates for b Xj and b Y j alone. These estimates can then be used to furnish a summary data MR analysis using the framework of a meta-analysis.
Under models (1) and (2), variant j is a valid instrument when a j ¼ 0 and invalid when a j 6 ¼ 0. When a j 6 ¼ 0, e. a bias term). In the Supplementary Methods (available as Supplementary data at IJE online), we briefly review three such summary data methods -IVW, 4 MR-Egger regression 8 and weighted median 9 -and discuss the conditions under which each method returns a consistent causal effect estimate (i.e. estimate converges in probability to the true value as the sample size increases).

The MBE
In this paper we propose a new causal effect estimatorthe MBE -that offers robustness to horizontal pleiotropy in a different manner to that of the IVW, MR-Egger or weighted median methods. Its ability to consistently estimate the true causal effect relies on the following fundamental assumption termed the ZEro Modal Pleiotropy Assumption (ZEMPA): across all instruments, the most frequent value (i.e., the mode) of b j is 0.
In order to formalize this, let k 2 f1; 2; . . . ; Lg represent the number of unique values of b j among the L variants. If all b j terms are identical then k ¼ 1, but if all are unique then k ¼ L. Now, let n 1 ; n 2 ; . . . ; n k represent the number of instruments that have the same non-zero value of b j , where n 1 represents those with the smallest non-zero identical value of b j and n k represents those with the largest non-zero identical value. Finally, let n 0 represent the number of valid instruments whose b j terms are identically zero. We then have that n 0 þ n 1 þ . . . þ n k ¼ L. ZEMPA implies that n 0 is larger than any other n l for l in 1; 2; . . . ; k (i.e., n 0 > maxðn 1 ; . . . ; n k Þ). For a weighted version of the MBE, that is an MBE derived by allowing the weight given to each ratio estimate to vary, ZEMPA implies that the weights associated with the valid instruments are the largest among all k subsets of instruments (ie. w 0 > maxðw 1 ; . . . ; w k Þ, where w l is the weight contributed by the lth subset of instruments using our previous subset definition based on b j . The breakdown level (i.e. the maximum proportion of information that can come from invalid instruments before the method is inconsistent) of the simple (i.e. unweighted) MBE ranges from 100 L=2þ1 L % to 100 LÀ2 L À Á %. The lower limit corresponds to the situation where there are some valid instruments, but all invalid instruments estimate the same (biased) causal effect parameter (i.e. k ¼ 2) implying that ZEMPA is satisfied (i.e. n 0 > maxðn 1 ; . . . ; n k Þ) if up to, but not including, half of the instruments are invalid. The upper limit corresponds to the situation where all invalid instruments estimate different causal effect parameters (i.e. n 1 ¼ n 2 ¼ . . . ¼ n k ¼ 1), implying that ZEMPA would be satisfied if just two variants were valid (n 0 ¼ 2Þ and the remainder (L À 2) were invalid. Given that maxðn 1 ; . . . ; n k Þ is often unknown and is likely to vary depending on the set of genetic instruments and the outcome variable, the true breakdown level of the MBE in any given applied investigation is difficult to determine.
For example, in Figure 1A, six out of eight instruments are invalid (so n 0 ¼ 2), but all non-zero b j s are unique, implying that k ¼ L À 1 ¼ 7 and n 1 ¼ n 2 ¼ . . . ¼ n 7 ¼ 1. In this situation, ZEMPA is satisfied and the simple MBE is a consistent estimate of the causal effect b. However, when the largest number of identical estimates comes from invalid instruments (i.e. n 0 < n l for some l; ZEMPA violated), then the simple MBE will be inconsistent for b (i.e. asymptotically biased). This is illustrated in Figure 1B, which shows causal effect estimates from six invalid and two valid variants (n 0 ¼ 2). Since three variants have precisely the same horizontal pleiotropic effect in this example (n 2 ¼ 3), ZEMPA is violated.
The breakdown level of the weighted MBE can be similarly defined as ranging from 50% (exclusive) to 100% (exclusive). In other words, the weighted MBE is biased if w 0 < w l for some l. Of note, the limits are open intervals because the weights are real numbers, unlike number of instruments (in the case of the simple MBE), which is a natural number. However, as L increases, then the lower and upper limits of the breakdown level of the simple MBE also tend to 50% and 100%, respectively.

Implementing the MBE
To calculate the MBE, we propose using the mode of the smoothed empirical density function of allb Rj s as the causal effect estimate. This strategy is straightforward to implement, easily deals with sampling variation in asymptotically identicalb Rj s and allows different weights to be given to different instruments. We refer to the mode of the unweighted and inverse-variance weighted empirical density function as the simple and weighted MBEs, respectively. The standardized weights for the weighted MBE can be computed as follows: For the simple MBE, Consider the normal kernel density function of theb Rj s: where h is the smoothing bandwidth parameter. 14 where sdðb RJ Þ and madb RJ are the standard deviation and median absolute deviation from the median of the L b Rj s, respectively. An intuitive explanation of the MBE based on an analogy with histograms is provided in the Supplementary Methods (available as Supplementary data at IJE online).

Simulation model
The simulations were performed using the following model to generate individual i's exposure X i , outcome Y i and confounder U i , based on their underlying genetic data vector (G i1 ; . . . ; G iL ): where: Z U , Z X and Z Y represent the additive allele scores of L independent SNPs on U, X and Y, modulated by the parameters d Uj ; d Xj ; d Y j (j ¼ 1,. . . L). b denotes the true causal effect of X on Y that we wish to estimate. The underlying genetic variables (G ij ) were generated independently by sampling from a Binomial (2, p) distribution with p itself drawn from a Uniform(0.1,0.9) distribution, to mimic biallelic SNPs in Hardy-Weinberg equilibrium. The resulting allele scores were then divided by their sample standard deviations r ZU ; r ZX ; r ZY ð ), to set variances to one. The direct effects of U on X and Y are denoted by h X and h Y , respectively. h X and h Y are set to positive values in all simulations, so as to always induce positive confounding. Error terms e Ui ; e Xi ; e Y i were independently generated from a normal distribution, with mean ¼ 0 and variances r 2 eU , r 2 eX and r 2 eY , respectively, whose values were chosen to set the variances of U, X and Y to one.
Constraining the variances in this way enables easy interpretation of the parameters in models (8)-(10). For example, b ¼ 0.1 implies that one standard deviation increment in X causes a 0.1 standard deviation increment in Y, and that the causal effect of X on Y explains 0.1 2 ¼ 1% of Y variance. A summary data interpretation of our simulation model is provided in the Supplementary Methods (available as Supplementary data at IJE online).

Simulation scenarios
Although the consistency property of an estimator provides a formal justification of the approach, it is equally important to understand how well it works in practice for realistically sized datasets in comparison with other methods. Therefore, we evaluated our proposed estimator in four different simulation scenarios. In all simulations, the number of variants . . . ; 30 is the number of invalid instruments.
Simulations 1 and 2 were aimed at evaluating the performance of the MBE under the causal null (b ¼ 0) in the two-sample setting. Datasets of 100 000 individuals were simulated and divided in half at random, and each was used to estimate either SNP-exposure or SNP-outcome associations. Simulations 3 and 4 were aimed at evaluating weak instrument bias in the twosample and single-sample settings; sample sizes used to estimate instrument-exposure (N X ) and instrumentoutcome (N Y ) associations were allowed to vary, as described below.
Simulation 1. In this scenario, d Uj was 0 for all instruments, implying that there is no InSIDE-violating horizontal pleiotropy. InSIDE-respecting horizontal pleiotropic effects d Y j were drawn from a Uniform(0.01, 0.2) distribution for the q invalid instruments or were set to 0 for valid instruments. Given that b ¼ 0, power can be interpreted as the type-I error rate.
Simulation 2. InSIDE-violating horizontal pleiotropy was induced by setting d Y j ¼ 0 for all instruments, whereas d Uj values were drawn from a Uniform(0.01, 0.2) distribution for the q invalid instruments.
Simulation 4. This simulation evaluated the performance of the estimators under the causal null when SNP-exposure and SNP-outcome associations are estimated in partially (50%) or fully (100%) overlapping samples (the latter being equivalent to the single sample setting). It was implemented as for simulation 3, except b ¼ 0 and N X ¼ N Y 2 {1 000, 5 000, 10 000}. We used smaller sample sizes to purposely increase the bias due to sample overlap, thus facilitating comparisons between methods.
Applied examples: plasma lipid fractions and urate levels and coronary heart disease risk Do and colleagues 16 performed a two-sample MR analysis to evaluate the causal effect of low-density lipoprotein cholesterol (LDL-C), high-density lipoprotein cholesterol (HDL-C) and triglycerides on coronary heart disease (CHD) risk, using a total of 185 genetics variants. Summary association results were obtained from the Global Lipids Genetics Consortium 17 and the Coronary Artery Disease Genome-Wide Replication and Meta-Analysis Consortium, 18 and were downloaded from Do and colleagues' supplementary material (standard errors were estimated based on the regression coefficients and P-values). Genetic variants were classified as instruments for each lipid fraction using a statistical criterion (P < 1 Â 10 À8 ), resulting in 73 instruments for LDL-C, 85 for HDL-C and 31 for triglycerides.
White and colleagues 19 performed a similar analysis, but with plasma urate levels rather than lipid fractions. 31 variants associated with urate levels (P < 5 Â 10 À7 ) were used as genetic instruments, and the required summary statistics were obtained from the GWAS catalogue [https:// www.ebi.ac.uk/gwas/].

Statistical analyses
In all simulation scenarios, causal effect estimates were obtained using established MR methods (multiplicative random effects IVW, 7 multiplicative random effects MR-Egger regression 7 and weighted median, all implemented using inverse-variance weights calculated under NOME), as well as the simple and the weighted MBEs. Each version of the MBE was evaluated using weights calculated with and without making the NOME assumption, thus yielding four MBEs. Each of these four methods was evaluated for two values of the tuning parameter u 2 f1; 0:5g, totalling eight versions of the MBE method. Parametric bootstrap was used to estimate the standard errors of the MBE using the median absolute deviation from the median (multiplied by 1.4826 for asymptotically normal consistency) of the bootstrap distribution of causal effect estimates. These were used to derive symmetrical confidence intervals.
In each scenario, coverage, power and average causal effect estimates, standard errors, F GX À1 F GX and I 2 GX statistics (which quantify the magnitude of violation of the NOME assumption in IVW and MR-Egger regression estimates, respectively 7,13 ) were obtained across 10 000 simulated datasets. Power was defined as the proportion of times that 95% confidence intervals excluded zero, and coverage as the proportion of times that 95% confidence intervals included the true causal effect.
MR methods were also applied to estimate the causal effect of plasma lipid fractions and urate levels on CHD risk. The magnitude of regression dilution bias in IVW and MR-Egger regression was assessed by the F GX À1 F GX and I 2 GX statistics, respectively. Cochran's Q test was used to test for the presence of horizontal pleiotropy (under the assumption that this is the only source of heterogeneity betweenb Rj s other than chance). 20 All simulations and analyses were performed using R 3.3.1 [www.r-project. org]. R code for implementing the MBE is provided in Supplementary Methods (available as Supplementary data at IJE online).

Performance under the causal null in the twosample context
The results of simulation 1 -where directional horizontal pleiotropy (if any) occurs only under the InSIDE assumption -are shown in Table 1. When all instruments were valid, all methods were unbiased with type-I error rates 5%. As expected, MR-Egger regression (which is consistent if InSIDE holds) was the least biased method in this scenario, especially when many instruments were invalid. The four MBEs in Table 1 were less biased and less precise than the IVW and the weighted median methods. The simple MBE was more biased than the weighted MBE (noticeable especially when the proportion of invalid instruments was high). Using weights derived under the NOME assumption increased bias and false rejection rates. Setting u ¼ 0.5 (i.e. setting the bandwidth to half of the default value) reduced both bias and precision (Supplementary Table 1, available as Supplementary data at IJE online).
When InSIDE is violated (Table 2), again the MBEs were less biased than IVW and weighted median methods. In this case, however, they were also less biased than MR-Egger regression estimates, which is known to be highly sensitive to InSIDE violation. 8 The exception was for large proportions (i.e. ! 80%) of invalid instruments, where MR-Egger estimates were the least biased. This is because the degree of InSIDE violation, as quantified by the inverse-variance weighted Pearson correlation between instrument strength and horizontal pleiotropic effects, 8 is smaller in those situations (Supplementary Table 2, available as Supplementary data at IJE online). Moreover, in this scenario, the simple MBE was generally less biased than the weighted counterparts, and setting u ¼ 0.5 had a smaller effect when compared with simulation 1 (and indeed only clear for the simple MBE-Supplementary Table 3). The NOME assumption again increased bias and false rejection rates.
Power to detect a causal effect in the two-sample context Table 3 displays the results for simulation 3 (no invalid instruments). The IVW method was the most powered to detect a causal effect, followed by the weighted median method, the weighted MBE, the simple MBE and MR-Egger regression. Assuming NOME reduced the bias towards the null in the weighted MBEs and improved power. Setting u ¼ 0.5 had no consistent effect on bias, but substantially reduced power (Supplementary Table 4, available as Supplementary data at IJE online).

Performance under the causal null in overlapping samples
Supplementary      precision of the MBE was very low, suggesting that the method may be prohibitively underpowered in small samples, thus being best suited for the two-sample setting using precise summary association results. Gains in precision by making the NOME assumption were more noticeable than in the other simulations with larger sample sizes.

Causal effect of plasma lipid fractions and urate levels on CHD risk
We used real datasets of summary association results to further explore the influence of the u parameter on the MBE. First, we visually explored the distribution of ratio estimates (Figure 2). In the case of LDL-C (panel A), most of the distribution was above zero, and increasing the stringency of u did not reveal substantial multimodality, although there were some pronounced density peaks at the left of the main distribution (which corresponds to the true causal effect under the ZEMPA assumption), which may result in attenuation of the causal effect estimate. However, setting u ¼ 0.25 resulted in some small peaks in the main distribution which may suggest over-stringency, so we used u ¼ 0.5 in the MR analysis. For HDL-C (panel B), the bulk of the distribution was centred close to zero, and setting u ¼ 0.25 revealed some peaks at the left of the main distribution, suggesting that horizontal pleiotropy could lead to an apparent protective effect. Since setting u ¼ 0.5 was sufficient to substantially reduce the density at the tails, this was used in the MR analysis. Regarding triglycerides (panel C), the main distribution was above zero ; sample size of the dataset used to estimate instrument-exposure associations; N Y ; sample size of the dataset used to estimate instrument-outcome associations; IVW, inverse-variance weighting; SE, estimated standard error; NOME, NO Measurement Error; MBE, mode-based estimate. a u ¼ 1. and the plot suggested that there may be negative horizontal pleiotropy, leading to an underestimation of the causal effect (u ¼ 0.25 was used in MR analysis). Finally, in the case of urate levels (panel D), by decreasing u it became increasingly evident that the distribution was bi-modal, which could only be clearly distinguished by setting u ¼ 0.25 (which was used in MR analysis) because the main peaks were similar to one another. Comparing the two distributions, the main one was the closest to zero, suggesting that horizontal pleiotropy is biasing the causal effect estimate upwards.
Results of the MR analysis are shown in Table 4. The smallest values of F GX À1 F GX and I 2 GX were 0.996 and 0.993, respectively, suggesting that IVW and MR-Egger regression estimates were not materially affected by regression dilution bias. P-values of the Cochran's Q test ranged from 0.0003 (urate) to 1.7 Â 10 À21 (HDL-C), thus providing strong statistical evidence for heterogeneity between the ratio estimates. Nevertheless, results for LDL-C and triglycerides consistently suggested risk-increasing causal effects. In the case of HDL-C, the IVW method suggested a protective effect, with one standard deviation increase in HDL-C being associated with a 0.254 (95% CI: 0.115; 0.393) decrease in CHD ln(odds). However, the other methods did not confirm this result, suggesting that it was due to negative horizontal pleiotropy (as suggested by visually inspecting the distribution of ratio estimates). Finally, the IVW method suggested a 0.163 (95% CI: 0.027; 0.298) increase in CHD ln(odds) per standard deviation increase in urate levels. Other methods did not confirm this finding, suggesting that it could be a result of positive horizontal pleiotropy (as the empirical density plot suggested).

Discussion
We have proposed a new MR method -the MBE -for causal effect estimation using summary data of multiple genetic instruments. Its performance was evaluated in a simulation study and its application illustrated in real data examples. An overview of the summary data MR methods that we evaluated (as well as the simple median) is provided in Table 5.
Consistent causal effect estimation using the MBE requires that ZEMPA holds. ZEMPA is an assumption that relates to the underlying bias parameters (the b j ) that contribute to the ratio estimand b j ¼ b þ b j identified by the jth genetic instrument. If ZEMPA is satisfied, then the MBE yields a consistent estimate for the causal effect. However, due to imprecision in theb j 's in finite samples, in practice the MBE may be contaminated by some invalid invariants even if ZEMPA holds. This can be seen in our simulations, where ZEMPA is only violated when all instruments are invalid, but nevertheless there is bias in the MBE when some of the instruments are valid. In practice, the MBE also depends on the magnitude of the bias, with invalid genetic instruments identifying causal effect parameters that are close to the true causal effect being more likely to contaminate the MBE estimate. However, this also means that genetic instruments that would introduce strong bias are less likely to contaminate the MBE.
In our simulations, we evaluated eight different versions of the MBE. Decreasing the tuning parameter u reduced bias (at the cost of reduced precision) when horizontal pleiotropy did not violate the InSIDE assumption. However, when InSIDE was violated, a similar behaviour could only be clearly seen for the simple MBE. Choosing the value of the tuning parameter u is a bias-variance trade-off and depends on how stringent the smoothing bandwidth needs to be and how stringent it can be before being prohibitively imprecise. In our applied example, we identified the stringency required through a graphical examination, and verified that the MBEs were powered enough to detect a causal effect between HDL-C and triglycerides on CHD risk. Moreover, in the case of urate levels, the weighted MBE was similarly precise to the IVW and weighted median methods. This suggests that it may be feasible to set u to stringent values in practice, especially when there are multiple instruments selected based on genome-wide significance. Evaluating a range of u values through a graphical examination may be useful to investigate how susceptible the MBE is to contamination from invalid instruments.
Assuming NOME increased bias and reduced the coverage of the 95% confidence intervals in the presence of invalid instruments, but reduced regression dilution bias and improved power in the two-sample setting. However, such gains were relatively small and virtually disappeared in simulations with larger sample sizes. Moreover, the results in the applied example were virtually identical whether or not NOME was assumed. These findings suggest that the NOME assumption is not necessary (and might be even unwarranted) when deriving weights for the MBE.
Although the simple MBE was less precise than the weighted MBE, it was less prone to bias due to violations of the InSIDE assumption. However, it was more prone to bias when InSIDE held. Indeed, a similar pattern has been previously shown for the simple and weighted median. 9 This suggests that comparing both methods would be a useful sensitivity analysis in practice, although care must be taken since the simple MBE may in some cases (as in our real data example with urate levels) be prohibitively imprecise. Importantly, all the recommendations above are general, and we strongly encourage researchers to consider study-specific factors when deciding upon these aspects. One way of doing so is to perform simulations that reflect the study-specific context and compare different thresholds and filters in a range of different scenarios, keeping observable parameters (e.g. sample size) constant. Such simulations would also be useful to identify how strong the violations of the assumptions must be in order to obtain the observed results, which may be a useful sensitivity analysis that will either strengthen or weaken causal inference.
In our simulations, the 95% confidence intervals of the MBE computed using the normal approximation presented over-coverage (i.e. coverage larger than 95%). This may be due to the MBE being less influenced by outlying instruments (which is indeed the basis of the method), which correspond to the most imprecise ones when all instruments are valid. Therefore, the causal effect estimate fluctuates less around the true causal effect b (i.e. is less influenced by sampling variation). This may also explain the less pronounced overcoverage in the weighted median. We compared the normal approximation with the percentile method (Supplementary  Table 7, available as Supplementary data at IJE online), but over-coverage in the latter was even greater when there were no or few invalid instruments. Moreover, after a certain proportion of invalid instruments (around 50%), coverage of the percentile method reduced markedly, whereas this occurred gradually in the normal approximation method. We therefore proposed the latter method to compute confidence intervals, but there might be better alternatives.
Another aspect of the MBE method (and of the weighted median) that requires further research is regression dilution bias in the two-sample setting. Understanding how regression dilution bias operates in IVW and MR-Egger contributed to developing correction methods, 13 thus reinforcing the importance of research in this area regarding the MBE and the weighted median.
Although this is the first description of using the MBE as a causal effect estimate in MR, other closely related methods have already been published. For example, Guo et al. 21 have recently described a method based on bivariate comparisons of all pairs of instruments, which classify instruments as estimating or not estimating the same causal effect. The largest identified set of concordant instruments can then be used to estimate the causal effect using, for example, the IVW method. Therefore, Guo et al.'s approach also relies on the assumption that the most common causal effect estimate is a consistent estimate of the true causal effect (i.e. ZEMPA). In fact, both our approach and Guo et al.'s can be viewed as methods that fully exploit the power of the consistency criterion defined originally by Kang et al, 22 who used it to propose a LASSObased variable selection procedure to detects and adjusts for horizontally pleiotropic variants. However, Guo et al.'s method and the MBE (which was developed independently from their work) are very different in their implementation. Ours is designed to be simple to understand and implement, does not require selecting instruments, and is easy to extend to any weighting scheme one desires. Moreover, plotting the empirical density function using different bandwidths may be a useful tool to visually explore the distribution of theb Rj s, and provides an intuitive way to select the optimal bandwidth value. In separate work we conduct a thorough review Guo et al.'s method after translating it to the two-sample context, and suggest some simple modifications to improve its performance. 23 It is also important to consider that there are other strategies to compute the mode of continuous data. In preliminary simulations, the modified Silverman's rule was both generally more robust against horizontal pleiotropy than the original Silverman's rule 24 and more powered to detect a causal effect. Therefore, we opted for the modified rule. However, many other kernels and bandwidth selection rules could be used, as well as strategies that are not based on the smoothed empirical density function, such as the simple and robust parametric estimators, 15 Grenander's estimators 25 and the half-sample mode method. 14 Further research is required to translate these mode estimators into the summary data MR context and compare their performance under different scenarios.
We propose the MBE as an additional MR method that should be used in combination with other approaches in a sensitivity analysis framework. Using several methods that make different assumptions, rather than a single method, is a useful strategy to assess the robustness of the results against violations of the instrumental variable assumptions. 26,27 Further developments in this area (including some aspects of the MBE itself) will contribute to expanding the arsenal of tools available to applied researchers to interrogate causal hypotheses with observational data.