Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
83,709 characters · 17 sections · 81 citation commands
Unbiased estimation and asymptotically valid inference in multivariable Mendelian randomization with many weak instrumental variables
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 \fi
\if01 \fi
{\it Keywords:} Causal Inference, Genome-Wide Association Studies, Inverse-Variance Weighting, Mendelian Randomization, Weak Instrumental Variables.
\spacingset{1.05}
A genome-wide association study (GWAS) refers to the identification of genetic variants statistically associated with complex traits or diseases across the whole genome using large population cohorts visscher201710. GWAS typically examines associations between single-nucleotide polymorphisms (SNPs) and a trait but can also handle other genetic variants such as insertion and deletions (indels) and structural variations (SVs) gresham2008comparing. The first example of a successful GWAS was the 2005 GWAS which revealed two genetic variants significantly associated with age-related macular degeneration klein2005complement. To date, over 5,000 human GWAS have investigated approximately 2,000 diseases and traits and have identified more than 400,000 genetic associations wijmenga2018importance. This groundbreaking work has uncovered numerous compelling associations with human complex traits and diseases, shedding light on the disease mechanisms and enhancing clinical care and personalized medicine tam2019benefits.
Mendelian randomization (MR) is an epidemiological method that utilizes genetic variants as instrumental variables (IVs) to infer whether an exposure causally influences an outcome burgess2021mendelian. Since the genotypes of individuals are randomly inherited from their parents and generally do not change during their lifetime, genetic variants are considered to be independent of underlying confounders and hence can be used as IVs to eliminate confounding bias. Early MR studies progressed slowly because individual-level data simultaneously including genotypes and phenotypes were rarely available ebrahim2008mendelian. Recently, many large GWAS have been published and the corresponding summary statistics are available in databases such as the GWAS Catalog macarthur2017new (\url{https://www.ebi.ac.uk/gwas/}), dbGaP (\url{https://www.ncbi.nlm.nih.gov/gap/}), and UK biobank (UKBB, sudlow2015uk) (\url{https://www.ukbiobank.ac.uk/}). The accuracy of causal effect estimation is improved and valuable insights into the causal relationships between common risk factors and diseases are uncovered by utilizing MR with GWAS summary data wang2022mendelian.
The inverse-variance weighted (IVW) method is the most popular approach used to perform MR with GWAS summary data. A causal effect estimate yielded by the IVW method is supposedly unbiased if three so-called valid IV conditions are satisfied: (IV1) the genetic variants are strongly associated with the exposure; (IV2) the genetic variants are associated with the outcome only through the exposure; and (IV3) the genetic variants are independent of confounders bowden2015mendelian. The directed acyclic graph (DAG) of valid IV conditions is shown in panel (a) in Figure (ref). Due to the complexity of genetic architecture, conditions IV2 and IV3 are often difficult to validate zhu2020mendelian. Meanwhile, it is challenging to quantify the instrument strength and define a universal criterion for concluding that an IV satisfies condition IV1, although the F statistic can be utilized as a rough measure of weak instrument bias burgess2011avoiding. Thus, quantifying and eliminating the bias of IVW estimate in MR analysis will lead to valid causal inference and help to understand disease etiology.
A genetic variant is termed a pleiotropic variant or pleiotropy if it simultaneously affects multiple traits through different pathways. There are two types of pleiotropy: vertical and horizontal pleiotropy, where the former refers to the genetic variant associated with one trait through the mediation of another trait (as described in panel (a) in Figure (ref)), while the latter refers to the genetic variant independently associated with both traits (as illustrated in panel (b) in Figure (ref)). IVs with evidence of horizontal pleiotropy should be removed before applying IVW because it violates either the (IV2) condition or the (IV3) condition; otherwise, a biased causal effect estimate is likely obtained. In the literature, there are three strategies to remove the effect of horizontal pleiotropy: 1) Identifying and excluding horizontally pleiotropic IVs by using hypothesis tests, such as the MR pleiotropy residual sum and outlier (MR-PRESSO, verbanck2018detection) and the iterative MR and pleiotropy (IMRP, zhu2021iterative); 2) Eliminating the effect of horizontal pleiotropy by applying robust tools; e.g., the MR-Egger bowden2015mendelian, MR-Median bowden2016consistent, and MR-Lasso/MR-Robust rees2019robust; 3) Automatically separating vertical pleiotropy from horizontal pleiotropy through a mixture mode, among which the representative methods include MR-Mix qi2019mendelian and MR contamination mixture (MR-ConMix, burgess2020robust).
It has been gradually realized that horizontal pleiotropy can be divided into uncorrelated horizontal pleiotropy (UHP) and correlated horizontal pleiotropy (CHP). UHP violates the (IV2) condition and usually refers to a genetic variant that is directly associated with the outcome. In contrast, CHP violates the IV3 condition and may occur when a genetic variant indirectly affects the outcome through the mediation of unspecified exposures. The DAG of UHP and CHP is shown in panel (b) in Figure (ref). morrison2020mendelian proposed causal analysis using summary effect (CAUSE), which is the first MR approach accounting for UHP and CHP simultaneously. cheng2022mr proposed MR-Corr to detect CHP by a Bayesian mixture model and cheng2022mendelian extended MR-Corr to MR-CUE (MR with CHP Unraveling shared Etiology and confounding) to detect the UHP and CHP simultaneously. Alternatively, xue2021constrained proposed the constrained maximum likelihood-based MR (cML-MR) method that identifies UHP and CHP through Bayesian information criterion (BIC, schwarz1978estimating). In addition, yuan2022likelihood derived MR with automated instrument determination (MRAID) to address UHP and CHP, which allows vertical pleiotropy to be in high linkage disequilibrium (LD).
A significant disadvantage of most existing approaches is that they assume both UHP and CHP to have similar properties as outliers in the traditional regression approach. However, there is substantial evidence that most traits have shared moderate or high genetic correlations bulik2015atlas, violating this technical assumption required by most existing approaches. Consequently, it is challenging to remove the effect of horizontally pleiotropic variants by considering only one exposure in MR analysis. Multivariable MR, which simultaneously estimates the causal effects of multiple exposures on an outcome, is compelling in resolving this problem burgess2015multivariable. Multivariable MR recognizes the bias caused by horizontal pleiotropy as an omitted-variable bias (OVB), which will disappear automatically if all the omitted exposures are specified in the multivariable MR model. The DAG of multivariable MR is exhibited in panel (c) in Figure (ref). So far, the multivariable versions of the IVW method, MR-Egger, MR-Median, and MR-Lasso/MR-Robust have been developed burgess2015multivariable,rees2017extending,grant2021pleiotropy. sanderson2019examination showed that the multivariable MR is able to unbiasedly estimate the causal effects of a target exposure when the other exposure is confounder, collider, or mediator of this exposure.
Weak instrument bias arises when the majority of IVs are weakly associated with the exposures, therefore violating the (IV1) condition and making conventional MR methods unreliable burgess2011avoiding. It is widely recognized that a common trait is often polygenic affected by hundreds or even thousands of independent variants/genes with small effect sizes. With the increasing sample sizes of GWAS, more and more trait-associated variants are being identified. Thus, the weak instrument bias is likely to become a considerable problem in future MR studies. burgess2011avoiding and sanderson2021testing suggested using the F and conditional F statistics to measure the weak instrument bias in MR and multivariable MR, respectively. burgess2016bias illustrated that the weak instrument bias also depends on the sample overlap in two-sample MR. sadreev2021navigating examined the impact of sample overlap and winner's curse when weak instrument bias exists and observed that the weak instrument bias grew dramatically in the presence of winner's curse. For the univariate MR model with no sample overlap, zhao2020statistical proposed the robust adjusted profile score to estimate the causal effect unbiasedly, while ye2021debiased provided the debiased IVW (DIVW) method to remove the weak instrument bias of IVW estimate. Overall, these aforementioned methods have neither provided a comprehensive theoretical analysis of weak instrument bias nor a general solution to remove the weak instrument bias in both univariable MR and multivariable MR.
As the first contribution of this paper, we theoretically characterize the bias in multivariable IVW causal estimate. Specifically, we demonstrate that the bias of IVW causal estimate is the product of weak instrument and estimation error biases. Meanwhile, we demonstrate that the estimation error bias is a linear combination of measurement error yi2017statistical and confounder biases, and the sample overlaps among multiple GWAS cohorts trade off the proportions of these two biases. With moderate conditions on the MR model, we theoretically illustrate how the number of IVs, sample sizes of GWAS studies, and sample overlap among GWAS cohorts influence the asymptotic behavior of multivariable IVW estimate. These theoretical findings are summarized in Theorem (ref) that to our best knowledge is the first comprehensive investigation of multivariable IVW estimate.
As the second contribution of this paper, we demonstrate our novel multivariable MR approach, MR using Bias-corrected Estimating Equations (MRBEE), can estimate causal effects unbiasedly in the presence of many weak IVs. Under moderate conditions, we investigate the asymptotic behaviors of IVW and MRBEE, revealing that MRBEE is superior to IVW in terms of strongly asymptotic unbiasedness. In particular, only when an estimate is strongly asymptotically unbiased, the inference made based on this estimate is asymptotically valid. Simulations show that only MRBEE can provide unbiased causal effect estimates in the presence of many weak IVs. Applied to data from the UK Biobank, MRBEE can successfully remove the weak instrument and estimation error biases and therefore make valid causal inferences.
This paper is arranged as follows. In section 2, we study the asymptotic behavior of multivariate IVW estimate. In section 3, we introduce MRBEE and examine its asymptotic properties. In section 4, simulations are conducted to compare MRBEE with the existing methods. In section 5, we apply MRBEE to estimate the causal effects of exposures on cardiovascular disease. Discussion is presented in section 6 and proofs of the related theorems are shown in Appendix. R package MRBEE (\url{https://github.com/noahlorinczcomi/MRBEE}) and supplementary materials are available online.
In this section, we introduce the notations, the model of the multivariable MR, and the bias of the multivariable IVW. Since univariable MR is a special case of multivariable MR, MR refers to the multivariable MR unless otherwise specified.
For a vector $\boldsymbol a=(a_j)_{p\times 1}$, $||\boldsymbol a||_q=(\sum_{j=1}^p|a_j|^q)^{1/q}$ with $q\in[0,\infty]$. For a symmetric matrix $\mathbf A=(A_{ij})_{p\times p}$, $\lambda_{\max}(\mathbf A)$ and $\lambda_{\min }(\mathbf A)$ is its the maximum and minimum eigenvalues, $\mathbf A^+$ is its Moore–Penrose inverse; and $||\mathbf A||_q=\max\{||\mathbf A\boldsymbol a||_q,\ ||\boldsymbol a||_q=1\}$. Let diag($\boldsymbol\alpha$) be the diagonalizing operator of vector $\boldsymbol\alpha$ and $\mathbf A\odot\mathbf B$ be the Hadamard product of matrices $\mathbf A$ and $\mathbf B$. For a set $\mathcal A$, $|\mathcal A|$ is the number of elements in $\mathcal A$. Notations $O(\cdot)$ and $o(\cdot)$ are the infinitely large and small quantities, while $O_P(\cdot)$ and $o_P(\cdot)$ mean that such relationships hold in probability.
The central aim of MR is to estimate causal effects between exposures and an outcome unbiasedly. Let $\boldsymbol g_i=(g_{i1},\dots,g_{im})^\top$ be an ($m\times 1$) genotype value vector of $m$ genetic variants, $\boldsymbol x_i=(x_{i1},\dots,x_{ip})^\top$ be an ($p\times 1$) vector representing $p$ exposures, and $y_i$ be an outcome. Here, $m$ is the number of specified IVs, which is usually the number of independent loci with $p$-values reaching the genome-wide significant level. Let $\mathbf B=(\boldsymbol\beta_1,\dots,\boldsymbol\beta_m)^\top$ be an ($m\times p$) matrix of genetic effects on exposures with $\boldsymbol\beta_j=(\beta_{j1},\dots,\beta_{jp})^\top$ being an $(p\times 1)$ vector, and $\boldsymbol\theta=(\theta_{1},\dots,\theta_p)^\top$ be an $(p\times1)$ vector of causal effects of the $p$ exposures on the outcome. The MR model is
where $\boldsymbol u_i$ and $v_i$ are the noise terms. Substituting for $\boldsymbol x_i$ in ((ref)), we obtain the equation
where $\boldsymbol\alpha=\mathbf B\boldsymbol\theta$. In the literature, ((ref)) - ((ref)) have been named as the structural form, first-stage, and reduced form, respectively stock2002survey.
In this paper, we assume that the total number of exposures $p$ is fixed and the causal effect $\boldsymbol\theta$ is fixed and bounded. The genetic variant $g_{ij}$ is standardized so that E$(g_{ij})=0$ and var$(g_{ij})=1$, and all IVs are in linkage equilibrium (LE), i.e., cov$(g_{ij},g_{ik})=0$ for $j\neq k$. The genetic effect $\boldsymbol\beta_j$ is random with zero-mean, covariance matrix $\bm\Sigma_{\beta\beta}$, and cumulative covariance matrix $\bm\Psi_{\beta\beta}$
The covariance matrix $\bm\Sigma_{\beta\beta}$ will vanish as $m$ increase, but the cumulative covariance matrix $\bm\Psi_{\beta\beta}$ is still a constant matrix, representing the total genetic covariance contributed from the $m$ IVs. Next, the noise terms $\boldsymbol u_i$ and $v_j$ have zero-means and joint covariance matrix \[ \bm\Sigma_{u\times v}=cov((\boldsymbol u_i^\top,v_j)^\top)=
\] Thus, the exposure vector $\boldsymbol x_i$ and outcome $y_i$ have zero-means and joint covariance matrix \[ \bm\Sigma_{x\times y}=cov((\boldsymbol x_i^\top,y_j)^\top)=
\] where $\mathbf\Sigma_{xx}=\bm\Psi_{\beta\beta}+\bm\Sigma_{uu}$, $\boldsymbol\sigma_{xy}=\bm\Psi_{\beta\beta}\boldsymbol\theta+\mathbf\Sigma_{uu}\boldsymbol\theta+\boldsymbol\sigma_{uv}$, and $\sigma_{yy}=\boldsymbol\theta^\top\bm\Psi_{\beta\beta}\boldsymbol\theta+\boldsymbol\theta^\top\bm\Sigma_{uu}\boldsymbol\theta+2\boldsymbol\theta^\top\boldsymbol\sigma_{uv}+\sigma_{vv}.$ Note that $\boldsymbol\sigma_{uv}\neq\mathbf 0$ means the confounders simultaneously affect $\boldsymbol x_i$ and $y_i$.
In genetics, the genetic effect $\beta_{js}$ can be treated as a random variable with mean 0 and variance $\psi_{\beta_s\beta_s}/m$, where $\psi_{\beta_s\beta_s}$ is the IV-heritability, i.e., the variance explained by additive effects of specified instrumental variants of the $s$th exposure bulik2015ld. Since a complex trait is often polygenic with a contribution from thousands of independent variants, the number of causal variants can be regarded as a number approaching infinity. Subject to this principle, a random effect model can describe the variation of these effects more simply and essentially. Although the fixed effect model has also been adopted by some works to study the asymptotic properties of the corresponding MR approaches zhao2020statistical,ye2021debiased, the random effect model is still the most commonly used at genome-wide level studies. Moreover, even if all causal variants were identified, the random effect model should still be more efficient than the fixed effect model to characterize the statistical property of the MR model diggle2002analysis.
The existing univariable MR methods, such as CAUSE morrison2020mendelian and MR-CUE cheng2022mendelian, can successfully remove the effects of CHP only when a small fraction of IVs have CHP, which is easy to violate because common traits may share a large fraction of pleiotropic variants bulik2015atlas. In contrast, multivariable MR resolves the pleiotropic variant problem by specifying all the relevant exposures in the model ((ref)), as the multivariable regression can automatically account for the pleiotropic variants shared by these exposures. This is one of the greatest advantages of multivariable MR over univariable MR. Hence, we assume that all the exposures can be included in the multivariable MR; therefore, the CHP effect is ignorable in ((ref)). On the other hand, some IVs may still present strong UHP effects. To account for potential UHP, we propose using an iterative procedure to remove these invalid IVs, which is similar to detecting outliers in MR analysis verbanck2018detection,zhu2020mendelian. Thus, the bias introduced by UHP and CHP can be greatly alleviated, as demonstrated in our proposed MRBEE.
With the rapid development of GWAS, large GWASs of exposures and disease outcomes have been conducted and their summary statistics including effect sizes, SEs, and variant information are publicly available for download sudlow2015uk,macarthur2017new. Thus, many recently developed MR methods are often designed based on GWAS summary statistics, as is in this paper.
With GWAS summary statistics, MR is mainly based on the linear regression model
where $\hat\alpha_j$ and $\hat{\boldsymbol\beta}_j$ are estimated from outcome and exposure GWAS for $j$th IV, and $\varepsilon_j$ represents the residual of this regression model. Let $\mathbf y^{[0]}=(y^{[0]}_1,\dots,y^{[0]}_{n_0})^\top$ be the sample vector from outcome GWAS, $\boldsymbol x^{[1]}=(x^{[1]}_1,\dots,x^{[1]}_{n_1})^\top,\dots,\boldsymbol x^{[p]}=(x^{[p]}_1,\dots,x^{[p]}_{n_p})^\top$ be the sample vectors of the 1st$,\dots,p$th exposure GWAS cohorts, and $\mathbf G^{[0]}=(g^{[0]}_{ij})_{n_0\times m},\dots,\mathbf G^{[p]}=(g^{[p]}_{ij})_{n_p\times m}$ be the sample matrices of $m$ genetic variants of the outcome and 1st$,\dots,p$th exposure GWAS cohorts. The sample size of the $s$th cohort is $n_s$, the overlapping sample size between the $s$th and the $k$th cohorts is $n_{sk}$, and the minimum sample size is $n_{\rm min}=\min\{n_0,\dots,n_p\}$.
The GWAS summary data are generated as follows. Suppose that $\boldsymbol y^{[0]}$, $\{\boldsymbol x^{[s]}\}$, and $\{\mathbf G^{[s]}\}$ are centered, and the $m$ genetic variants are in LE, i.e. ${\rm E}(\mathbf G^{[s]\top}\mathbf G^{[s]}/n_s)=\mathbf I_m$ for $j=0,1,\dots,p$. This orthogonality enables the following genetic effects to be estimated separately
where the corresponding variance estimates are given by
Then the GWAS summary data are formed by $\hat{\boldsymbol\alpha}=(\hat\alpha_1,\dots,\hat\alpha_m)^\top$, $\hat{\boldsymbol\beta}_j=(\hat\beta_{j1},\dots,\hat\beta_{jp})^\top$, $\hat{\mathbf B}=(\hat{\boldsymbol\beta}_1,\dots,\hat{\boldsymbol\beta}_m)^\top$, the related SE estimates, the p-values, and sample sizes $n_0,n_1,\dots,n_p$, and SNPs information.
The IVW method is equivalent to a weighted regression which estimates $\boldsymbol\theta$ by
where $\mathbf V=\text{diag}(1/\text{var}(\hat\alpha_1),\dots,1/\text{var}(\hat\alpha_m))$. In practice, we often standardize $\hat\alpha_j$ and $\hat{\beta}_{js}$ by $\hat{\alpha}_j/\text{se}(\hat{\alpha}_j)$ and $\hat{\beta}_{js}/\text{se}(\hat{\beta}_{js})$ to remove the minor allele frequency effect zhu2022genome. Therefore, ${\text{var}}(\hat\alpha_j)=1$ for all $j$ and ((ref)) reduces to
Here, we qualitatively show that the IVW estimate is biased due to the estimation errors of $\hat{\boldsymbol\alpha}$ and $\hat{\mathbf B}$, i.e., ${\boldsymbol w}_\alpha=\hat{\boldsymbol\alpha}-\boldsymbol\alpha$ and ${\mathbf W}_\beta=\hat{\mathbf B}-\mathbf B$, and meanwhile, the weak IVs can inflate th estimation error bias. Specifically, consider the estimating equation and Hessian matrix of $\hat{\boldsymbol\theta}_{\rm IVW}$:
That is, $\boldsymbol{S}_{\rm IVW}(\boldsymbol\theta)$ is the score function of ((ref)) and $\hat{\boldsymbol\theta}_{\rm IVW}$ is estimated by solving $\boldsymbol{S}_{\rm IVW}(\hat{\boldsymbol\theta}_{\rm IVW})=\mathbf 0$, and $\textbf{H}_{\rm IVW}$ is the second order derivative matrix of ((ref)). In particular, since the third derivative of the quadratic loss function ((ref)) is zero, we have $\hat{\boldsymbol\theta}_{\rm IVW}-\boldsymbol\theta=-\textbf{H}_{\rm IVW}^{-1}\boldsymbol{S}_{\rm IVW}(\boldsymbol\theta)$. As a result, the expectation of the bias of $\hat{\boldsymbol\theta}_{\rm IVW}$ is approximately:
where ${\boldsymbol w}_{\beta_j}$ is the $j$th row of ${\mathbf W}_\beta$, $w_{\alpha_j}$ is the $j$th element of ${\boldsymbol w}_\alpha$, and \[ cov(({\boldsymbol w}_{\beta_j}^\top,w_{\alpha_j})^\top)=\bm\Sigma_{W_\beta \times w_\alpha}=
, \] Intuitively, the bias of $\hat{\boldsymbol\theta}_{\rm IVW}$ has a product structure “weak instrument bias $\times$ estimation error bias". We call $\{\bm\Sigma_{W_\beta W_\beta}\boldsymbol\theta-\boldsymbol\sigma_{W_\beta w_\alpha}\}$ the estimation error bias because it comes from the covariance matrix of estimation errors $\bm\Sigma_{W_\beta \times w_\alpha}$. We term $\{\mathbf\Sigma_{\beta\beta}+\bm\Sigma_{W_\beta W_\beta}\}$ the weak instrument bias because the bias of $\hat{\boldsymbol\theta}_{\rm IVW}$ is inflated if the covariance matrix of effect sizes $\mathbf\Sigma_{\beta\beta}$ is not considerably larger than the covariance matrix of estimation errors $\bm\Sigma_{W_\beta W_\beta}$, which often happens if the majority of IVs used to infer the causal effect have weak effects.
In this subsection, we investigate the asymptotic behavior of the IVW estimate as the number of IVs $m$ and the minimum sample size $n_{\min }$ go to infinity. To facilitate the theoretical derivation, we specify the following three definitions and four regularity conditions.
A sub-Gaussian variable is one of the basic concepts in modern statistics vershynin2018high. It generalizes the scope of ordinary Gaussian variables to include all bounded discrete and common continuous variables. The well-conditioned covariance matrix is another important concept bickel2008regularized. A well-conditioned covariance matrix will ensure that the related statistical optimization is nondegenerate. In addition, we define the strongly asymptotic unbiasedness to distinguish the consistent estimate whose bias square vanishes with an equal and a smaller rate than its variance, respectively. If an estimate is consistent but its bias square and variance vanish at the same rate, the classic confidence interval cannot cover the true parameter with a probability of 0.95, thus leading to invalid statistical inference. This problem widely exists in all fields of statistics, especially, in nonparametric statistics and high-dimensional statistics, and many novel methods are derived to reduce the bias such that the bias square vanishes faster than the variance hall1992effect,van2014asymptotically,jankova2018semiparametric,calonico2018effect.
Conditions (C1)-(C4) restrict that all variables involved in this paper are sub-Gaussian distributed. In practice, $g_{ij}$ is standardized from a binomial variable with status 0, 1, and 2. Hence, it is supposedly a bounded sub-Gaussian variable as long as its minor allele frequency is not rare. Besides, we assume $\sqrt m\boldsymbol \beta_j$ to be sub-Gaussian with a well-conditioned covariance matrix $\bm\Psi_{\beta\beta}$ because the covariance explained by each variant $\bm\Sigma_{\beta\beta}$ decreases as the number of instrumental variants $m$ increases.
Theorem (ref) demonstrates the asymptotic normal distribution of the estimation errors, based on which we are able to obtain
where the $(j,s)$th element of $\mathbf\Delta_{xx}$ is $n_{js}/(n_jn_s)$ and the $j$th element of $\boldsymbol\delta_{xy}$ is $n_{j0}/(n_0n_j)$. As a result, the expectations of $\boldsymbol S_{\rm IVW}(\boldsymbol\theta)$ and $\mathbf H_{\rm IVW}$ are given by
By expressing $\boldsymbol\sigma_{xy}=\mathbf\Sigma_{xx}\boldsymbol\theta+\boldsymbol\sigma_{uv}$, we obtain an alternative expectation of $\boldsymbol S_{\rm IVW}(\boldsymbol\theta))$:
From this expectation, it is clear that there are two sources of the estimation error bias: $\{(\mathbf\Delta_{xx}-\boldsymbol\delta_{xy}\mathbf 1^\top)\odot\bm\Sigma_{xx}\}\boldsymbol\theta$ comes from the measurement error, while $\{\boldsymbol\delta_{xy}\odot\boldsymbol\sigma_{uv}\}$ is caused by the confounder. Here, we call $\{(\mathbf\Delta_{xx}-\boldsymbol\delta_{xy}\mathbf 1^\top)\odot\bm\Sigma_{xx}\}\boldsymbol\theta$ the measurement error bias because it has the same statistical impact, i.e., shrinking the coefficient estimate toward zero, as in measurement error analysis yi2017statistical. In contrast, we term $\{\boldsymbol\delta_{xy}\odot\boldsymbol\sigma_{uv}\}$ the confounder bias because $\boldsymbol\sigma_{uv}\neq\mathbf0$ implies that there are underlying confounders simultaneously affecting both $\boldsymbol x_i$ and $y_i$. In addition, the overlapping fraction vector $\boldsymbol\delta_{xy}$ trades off these two sources of biases. Generally, the measurement error bias is dominant when the elements of $\boldsymbol\delta_{xy}$ are small, while the confounder bias dominates when the elements of $\boldsymbol\delta_{xy}$ are large, and there may exist a special sample overlap such that $\boldsymbol\delta_{xy}\odot\boldsymbol\sigma_{uv}=\{(\mathbf\Delta_{xx}-\boldsymbol\delta_{xy}\mathbf 1^\top)\odot\bm\Sigma_{xx}\}\boldsymbol\theta$. In univariable MR, this special fraction is $n_{01}/n_0=\sigma_{xx}\theta/\sigma_{xy}$, which guarantees that E$(S_{\rm IVW}(\theta))=0$ and E$(\hat\theta_{\rm IVW})=\theta$. This theoretical result explains why in the empirical studies (e.g., Figures 1 and 2 in sadreev2021navigating), $\hat{\theta}_{\rm IVW}$ has a negative bias when $n_{01}/n_0$ is small, positive bias when $n_{01}/n_0$ is large, and is unbiased at this specific point.
Theorem (ref) is one of two main theorems in this paper and points out four scenarios. First, if $m$ goes to infinity with a lower rate than $\sqrt{n}_{\rm min}$, $\hat{\boldsymbol\theta}_{\rm IVW}$ is strongly asymptotically unbiased. In other words, $\hat{\boldsymbol\theta}_{\rm IVW}$ is able to reliably infer causality only when the sample size of GWAS data is quadratically larger than the number of IVs. On the other hand, the asymptotic covariance matrix of $\hat{\boldsymbol\theta}_{\rm IVW}$ is the inverse of the cumulative covariance matrix $\mathbf\Psi_{\beta\beta}=\sum_{j=1}^m\text{cov}(\boldsymbol\beta_j)$, therefore, it is optimal to include as many associated variants as possible in order to have $\mathbf\Psi_{\beta\beta}$ large enough. In contrast, using a few top significant variants to perform MR analysis is not recommended.
Second, if $m$ tends to infinity with the same rate as $\sqrt{n}_{\rm min}$, $\sqrt n_{\min }(\hat{\boldsymbol\theta}_{\rm IVW}-\boldsymbol\theta)$ converges to an asymptotic normal distribution with a non-zero asymptotic bias $\{-c_0\bm\Psi_{\beta\beta}^{-1}(\bm\Psi_{W_\beta W_\beta}\boldsymbol\theta-\boldsymbol\psi_{W_\beta w_\alpha})\}$. In this asymptotic bias, $\{-c_0(\bm\Psi_{W_\beta W_\beta}\boldsymbol\theta-\boldsymbol\psi_{W_\beta w_\alpha})\}$ is caused by $\boldsymbol S_{\rm IVW}(\boldsymbol\theta)$ and $\bm\Psi_{\beta\beta}^{-1}$ is resulted by $\mathbf H_{\rm IVW}^{-1}$. Since the asymptotic bias and asymptotic covariance matrix are of the same order in this scenario, the inference made is invalid although the bias of $\hat{\boldsymbol\theta}_{\rm IVW}$ is infinitesimal. Scenario $(iii)$ is more serious than $(ii)$ because the bias of $\hat{\boldsymbol\theta}_{\rm IVW}$ will not vanish even when $\sqrt n_{\rm min}$ goes to infinity. In the fourth scenario, $\hat{\boldsymbol\theta}_{\rm IVW}$ converges to a term irrelevant to $\boldsymbol\theta$. Scenarios (ii) - (iv) indicate that the IVW method is unlikely to make valid causal inference unless the sample sizes are quadratically larger than the number of IVs.
It is crucial to understand the asymptotic behaviors of $\hat{\boldsymbol\theta}_{\rm IVW}$ since the IVW method serves as the foundation for practically all MR techniques. Specifically, IMRP and MR-PRESSO use hypothesis tests to identify invalid IVs and then apply the IVW method to estimate causal effects based on valid IVs only. MR-Robust and MR-Median replace the quadratic loss function used in IVW by a robust loss function and absolute loss function, respectively. Although there have been literature studying the bias of $\hat{\boldsymbol\theta}_{\rm IVW}$ empirically burgess2011avoiding,burgess2016bias, they could not explain what causes the bias and how it behaves asymptotically. In contrast, Theorem (ref) points out the asymptotic properties of $\hat{\boldsymbol\theta}_{\rm IVW}$, representing a significant advance in understanding the IVW method and its extensions.
According to ((ref)), it is possible to remove the bias of $\boldsymbol S_{\rm IVW}(\boldsymbol\theta)$ by subtracting the measurement error bias $\{\bm\Sigma_{W_\beta W_\beta}\boldsymbol\theta-\boldsymbol\sigma_{W_\beta w_\alpha}\}$. Motivated by this principle, we propose MRBEE that estimates the causal effect estimates by solving the new unbiased estimating equation. In this section, we introduce the estimation of MRBEE, investigate its asymptotic properties, and discuss three implementation issues including the estimations of the bias-correction terms, the estimation of sandwich formula of causal effect estimate, and the detection of potential pleiotropy.
There are many methods that can remove the measurement error bias, including maximum likelihood estimation, unbiased estimating functions, and simulation-extrapolation (SIMEX) methods; see, e.g., yi2017statistical. MRBEE is a subtraction correction method belonging to the class of unbiased estimating function methods. Specifically, MRBEE estimates $\boldsymbol\theta$ by solving the following unbiased estimating equation:
where $\boldsymbol{S}_{\rm IVW}(\boldsymbol\theta)=-\hat{\mathbf B}^\top(\hat{\boldsymbol\alpha}-\hat{\mathbf B}\boldsymbol\theta)/m$. The solution $\hat{\boldsymbol\theta}_{\rm BEE}$ such that $\boldsymbol{S}_{\rm BEE}(\hat{\boldsymbol\theta}_{\rm BEE})=\mathbf0$ is
In practice, $\hat{\boldsymbol\theta}_{\rm BEE}$ is unreliable when the minimum eigenvalue of $\hat{\mathbf B}^\top\hat{\mathbf B}/m-\bm\Sigma_{W_\beta W_\beta}$ is negative, which is also a common problem for subtraction correction methods. In this case, we recommend first adjusting the negative eigenvalues to be 0 and then using the generalized inverse of this semi-positive matrix to yield $\hat{\boldsymbol\theta}_{\rm BEE}$.
Theorem (ref) indicates the following three scenarios. First, if $m/n\to0$,$\sqrt n_{\rm min}(\hat{\boldsymbol\theta}_{\rm BEE}-\boldsymbol\theta)$ converges to a normal distribution with a zero mean and the covariance matrix being exactly the same as $\hat{\boldsymbol\theta}_{\rm IVW}$. In other words, $\hat{\boldsymbol\theta}_{\rm BEE}$ not only enjoys the strongly asymptotic unbiasedness but also loses no efficiency in comparison to $\hat{\boldsymbol\theta}_{\rm IVW}$. Second, if $m/n_{\min }\to c_0\in(0,\infty)$, there is an additional covariance matrix $c_0\bm\Psi_{\beta\beta}^{-1}\bm\Psi_{\rm BC}\bm\Psi_{\beta\beta}^{-1}$ in the asymptotic normal distribution, where $\bm\Psi_{\rm BC}$ is introduced by the bias-correction terms: \[ \bm\Psi_{\rm BC}=\lim_{n_{\rm min}\to\infty}\text{var}\bigg[\frac{n_{\rm min}}{\sqrt m}\bigg((\mathbf W_\beta^\top\mathbf W_\beta-m\mathbf\Sigma_{W_\beta W_\beta})\boldsymbol\theta-(\mathbf W_\beta^\top\boldsymbol w_\alpha-m\boldsymbol\sigma_{W_\beta w_\alpha})\bigg)\bigg]. \] In this scenario, $\hat{\boldsymbol\theta}_{\rm BEE}$ is again strongly asymptotically unbiased with a convergence rate $\sqrt n_{\rm min}$, while $\hat{\boldsymbol\theta}_{\rm IVW}$ suffers from a bias not vanishing asymptotically. In the third scenario, $\hat{\boldsymbol\theta}_{\rm BEE}$ is still strongly asymptotically unbiased with a convergence rate $\sqrt{n_{\rm min}^2/m}$, and the asymptotic distribution is dominated by the bias correction term. In contrast, $\hat{\boldsymbol\theta}_{\rm IVW}$ converges to a term irrelevant to $\boldsymbol\theta$. Note that $\hat{\boldsymbol\theta}_{\rm IVW}$ is not consistent unless $m/n\to0$ and the inference made by $\hat{\boldsymbol\theta}_{\rm IVW}$ is unreliable unless $m/\sqrt n_{\rm min}\to0$. Therefore, MRBEE is superior to IVW in terms of both unbiasedness and asymptotic validity.
Most previous works of MR introduced their methods from the perspective of empirical applications and have not discussed the asymptotic properties; see, e.g., bowden2015mendelian,bowden2016consistent,verbanck2018detection,morrison2020mendelian. Some works zhao2020statistical,ye2021debiased described the asymptotic behaviors of the causal effect estimates yielded by their univariate MR methods, but the convergence rates and related conditions were not straightforward. For example, zhao2020statistical showed that the convergence rate of their causal effect estimate is $O(V_1/\sqrt{V_2})$ where $V_1$ and $V_2$ are two $m$-concentrations, which may mislead that this estimate has a $O(\sqrt m)$ convergence rate. From Theorem (ref), it is easy to see that $\hat{\boldsymbol\theta}_{\rm BEE}$ is strongly asymptotically unbiased, the asymptotic covariance matrix is $\psi_\theta\bm\Psi_{\beta\beta}^{-1}$, $\psi_\theta\bm\Psi_{\beta\beta}^{-1}+c_0\bm\Psi_{\beta\beta}^{-1}\bm\Psi_{\rm BC}\bm\Psi_{\beta\beta}^{-1}$, and $\bm\Psi_{\beta\beta}^{-1}\bm\Psi_{\rm BC}\bm\Psi_{\beta\beta}^{-1}$, and the convergence rate is $\sqrt n_{\rm min}$, $\sqrt n_{\rm min}$, and $\sqrt{n_{\min}^2/m}$, with respect to scenarios $(i)$, $(ii)$, and $(iii).$ In addition, although our method focuses on the multivariable MR model, the theoretical results can be readily extended to the univariable MR model. To the best of our knowledge, this is the first theoretical work to demonstrate how the convergence rate and asymptotic normal distributions vary with the sample sizes of multiple GWAS cohorts and the number of IVs for univariable and multivariable MR.
In this subsection, we discuss how to estimate the bias-correction terms $\bm\Sigma_{W_\beta W_\beta}$ and $\boldsymbol\sigma_{W_\beta w_\alpha}$ in practice. Specifically, we apply the method provided by zhu2015meta to estimate the covariance matrix $\bm\Sigma_{W_\beta\times w_\alpha}$ of the vector $(\boldsymbol w_{\beta_j}^\top,w_{\alpha_j})^\top$ from insignificant GWAS summary statistics. Let $\mathbf G^{\{0\}}=(g^{\{0\}}_{ij})_{n_1\times M},\dots,\mathbf G^{\{p\}}=(g^{\{p\}}_{ij})_{n_s\times M}$ be the sample matrices of $M$ insignificant and independent genetic variants. The insignificance means that the $p$-value of the genetic variants are larger than 0.05 for all exposures and outcome, and independence means that these variants are in LE. The insignificant GWAS statistics are estimated by
for $s=1,\dots,p$. With these insignificant effect sizes, $\bm\Sigma_{W_\beta \times w_\alpha}$ can be estimated by
because $\hat\alpha_j^*$ and $\hat{\beta}^*_{js}$ follow the same distributions of $w_{\alpha_j}$ and $w_{\beta_{js}}$, respectively. Here, $\hat{\bm\Sigma}_{W_\beta W_\beta}$ is the first $(p\times p)$ sub-matrix of $\hat{\bm\Sigma}_{W_\beta \times w_\alpha}$ and $\boldsymbol\sigma_{W_\beta w_\alpha}$ consists of the first $p-1$ elements of the last column of $\hat{\bm\Sigma}_{W_\beta \times w_\alpha}$.
Theorem (ref) shows that $\hat{\bm\Sigma}_{W_\beta \times w_\alpha}$ has a $O(\sqrt M)$ convergence rate after adjusting the scale of $\bm\Sigma_{W_\beta \times w_\alpha}$. As there may be more than 1 million independent variants in the whole genome, $\hat{\bm\Sigma}_{W_\beta \times w_\alpha}$ has high precision. In addition, $n_0,n_1,...,n_p\to\infty$ are required such that $\sqrt{n_0}\hat\alpha_j^*$ and $\sqrt{n_s}\hat{\beta}^*_{js}$ are asymptotically normally distributed. In addition, many popular GWAS methods such as cross-phenotype association analysis (CPASSOC, zhu2015meta) and multi-trait analysis of GWAS (MTAG, turley2018multi) need to estimate the covariance matrix of the estimation errors of GWAS summary statistics. As far as we are concerned, this theorem is the first one to theoretically guarantee that this covariance matrix can be consistently estimated from the GWAS insignificant statistics.
In this subsection, we illustrate how to estimate the covariance matrix of $\hat{\boldsymbol\theta}_{\rm BEE}$, i.e., cov($\hat{\boldsymbol\theta}_{\rm BEE}$)=$\bm\Sigma_{\rm BEE}(\boldsymbol\theta)$, through the famous sandwich formula liang1986longitudinal:
Here, the outer matrix $\textbf{F}_{\rm BEE}$ is the Fisher information matrix, i.e., the expectation of the Hessian matrix of $\boldsymbol{S}_{\rm BEE}(\boldsymbol\theta)$:
The inner matrix $\textbf{V}_{\rm BEE}(\boldsymbol\theta)$ is the covariance matrix of $\boldsymbol{S}_{\rm BEE}(\boldsymbol\theta)$:
where
A consistent estimate of $\bm\Sigma_{\rm BEE}(\boldsymbol\theta)$ is
where
and $\hat{\bm\Sigma}_{W_\beta W_\beta}$ and $\hat{\boldsymbol\sigma}_{W_\beta w_\alpha}$ are estimated through ((ref)).
Theorem (ref) shows that $\hat{\bm\Sigma}_{\rm BEE}(\boldsymbol\theta)$ has a $\min(\sqrt n_{\rm min},\sqrt{n^2_{\rm min}/m},\sqrt{m/\log m})$ convergence rate when $m/n^2_{\rm min}\to 0$. The first two convergence rates are brought by $||\hat{\mathbf F}_{\rm BEE}-\mathbf F_{\rm BEE}||_2$, while the third convergence rate is yielded by $||\hat{\textbf{V}}_{\rm BEE}(\hat{\boldsymbol\theta}_{\rm BEE})-\textbf{V}_{\rm BEE}(\boldsymbol\theta)||_2$, where the non-asymptotic analysis tool of random matrices are used to derive them vershynin2018high. Note that the SE estimation should be of the same importance as the causal effect estimation. Although the inference is made based on an unbiased estimate, it could still be invalid if the SE estimate is not reliable. Our simulations show that the vast majority of current univariable and multivariable MR approaches are unable to provide accurate SE estimates, e.g., MR-median consistently overestimates the SE and others have a tendency to underestimate it. In contrast, the sandwich formula, whose dependability has been extensively investigated empirically, is a reliable technique to obtain the SE estimate for MRBEE. This is yet another advantage of MRBEE over current approaches.
Due to the complexity of GWAS data, we cannot completely rule out the possibility of the existence of UHP and CHP even in the case of modeling multiple exposures. Specifically, if UHP and CHP exist,
where $\gamma_{u_j}$ is a UHP satisfying E($\gamma_{u_j}\boldsymbol\beta_j)=\mathbf0$ and $\gamma_{c_j}$ is a CHP satisfying E($\gamma_{c_j}\boldsymbol\beta_j)\neq\mathbf0$. Conventional pleiotropy detection methods such as MR-Robust, MR-PRESSO, and IMRP do not distinguish between UHP and CHP as long as they resemble outliers. Recently, some novel methods such as CAUSE and MR-CUE have been developed to separate vertical pleiotropy, UHP and CHP by using a mixture model, allowing slightly larger proportions of UHP and CHP. However, both the conventional and novel methods only focus on one exposure, failing to realize that most CHP and UHP may disappear automatically after specifying all the relevant exposures.
In this paper, we assume that we have excluded all CHP by including all the relevant exposures and we adopt IMRP zhu2021iterative to detect UHP. First, we define UHP as
In particular, we assume that $\gamma_j$ has a product structure $\gamma_j=\gamma_j^*b_j$, where $\gamma^*_j$ is a fixed number and $b_j$ is a non-random binary indicator. Let $\mathcal O=\{j:\ b_j\neq0\}$ be the set of UHP. The number of elements in $\mathcal O$ (i.e., $|\mathcal O|$) should be relatively small, otherwise the UHP cannot be regarded as outliers. We specify the following variant-specific hypothesis test:
A natural estimate of $\gamma_j$ is
where $\epsilon_j=w_{\alpha_j}-{\boldsymbol w}_{\beta_j}^\top\boldsymbol\theta+{\boldsymbol w}_{\beta_j}^\top(\hat{\boldsymbol\theta}_{\rm BEE}-\boldsymbol\theta)$. It is easy to see that $\text{E}(\epsilon_j)=0$ and $ \text{var}(\epsilon_j)=\boldsymbol\theta^\top\mathbf\Sigma_{W_\beta w_\alpha}\boldsymbol\theta+\sigma_{\omega_\gamma\omega_\gamma}-2\boldsymbol\theta^\top\boldsymbol\sigma_{W_\beta w_\alpha}.$ As a result, $t_{\gamma_j}=\hat\gamma^2_j/\text{var}(\epsilon_j)$ can be chosen as a feasible testing statistic for the hypothesis in ((ref)), which follows a central $\chi^2_1$-distribution under the null hypothesis. In practice, $\text{var}(\epsilon_j)$ can be estimated by
where $\hat{\boldsymbol\vartheta}_{\rm BEE}=(\hat{\boldsymbol\theta}^\top_{\rm BEE},-1)^\top$, $\mathbf{SE}_j=\text{diag}(\text{se}(\hat\beta_{j1}),\dots,\text{se}(\hat\beta_{jp}),\text{se}(\hat\alpha_j))$, and $\hat{\mathbf R}_{W_\beta \times w_\alpha}$ is the correlation matrix of $\hat{\mathbf\Sigma}_{W_\beta \times w_\alpha}$. Then $\gamma_j$ is considered as an outlier if
where $F_{\chi^2_1}(\cdot)$ is the CDF of $\chi^2_1$-distribution, $\hat t_{\gamma_j}=\hat\gamma_j^2/\widehat{\text{var}}(\epsilon_j)$, and $\kappa$ is a given threshold.
Theorem (ref) indicates that there is a theoretical threshold $\kappa= F_{\chi^2_1}(C_0\log m)$ to consistently identify all UHP. This threshold increases with a rate $O(\log m)$ to reduce the false discovery rate (FDR) and its concrete value can be chosen by a FDR control method benjamini1995controlling. In practice, MRBEE will iteratively apply the hypothesis test ((ref)) to remove the outliers and use the remaining IVs to estimate $\boldsymbol\theta$. The stable estimate is regarded as $\hat{\boldsymbol\theta}_{\rm BEE}$.
In this section, we conduct numerical comparisons between MRBEE and existing MR methods. Full details of simulation settings and additional simulation results are shown in the supplementary material.
We briefly introduce the simulation settings for univariable MR. First, we generate a binomial variable from $\text{Binom}(2,b_j)$ where $b_j\sim\text{Unif}(0.05,0.5)$ and standardize it as $g_{ij}$, the direct effect $\beta_j$ from $\mathcal{N}(0,1/m)$, and $u_i,v_i$ from a normal distribution with correlation coefficient $0.5$. The variances of $u_i$ and $v_i$ are chosen such that the IV-heritabilities are $\sigma_{\beta\beta}/\sigma_{xx}=0.3$ and $\theta^2\times(\sigma_{\beta\beta}/\sigma_{yy})=0.15$, respectively. We specify the causal effect $\theta=0.3/\sqrt{2}$. We compare MRBEE with IVW, DIVW, MR-Egger, MR-Lasso, MR-Median, IMRP, MR-ConMix, and MR-MiX, where most are implemented by using the R package MendelianRandomization yavorska2017mendelianrandomization. Additionally, the IMRP procedure is incorporated into MRBEE in which the threshold $\kappa$ is chosen by R package FDRestimation murray2020false. The so-called overlapping fraction is $n_{01}/n_0$, where the special fraction such that $\text{E}(S_{\rm IVW}(\theta))=0$ is $n_{01}/n_0\approx0.77.$ The number of independent replications is 1000.
First, we study the influences of overlapping fraction $n_{01}/n_0$ and the number of IVs $m$, with the results displayed in Figure (ref). Here, we fix $n_0=n_1=20000$, specify $n_{01}$ according to the overlapping fraction, and assume no UHP or CHP. It is easy to see that in general, only MRBEE is able to yield an unbiased estimate of $\theta$. For a special overlapping fraction (placed in the second column of Figure (ref)), all approaches become unbiased except DIVW. DIVW performs badly because it will further remove IVs based on their significance levels and consequently introduces an extra IV selection bias. In addition, the SE of causal effect estimate for all methods increases as the overlapping fraction decreases but remains unchanged by the increase of $m$. The results are consistent with our theoretical expectation and asymptotic properties of MRBEE.
As for the standard error, we display the boxplot of $\hat{\text{se}}(\hat\theta)-\text{se}(\hat\theta)$ where $\text{se}(\hat\theta)$ is approximated by the empirical SE calculated from the independent replications. It is evident that the SE estimates produced by all approaches have reduced variances as $m$ grows. However, only MRBEE and DIVW can provide consistent SE estimates, confirming the accuracy of MRBEE and DIVW's SE formulas. Additionally, MR-ConMix is extremely likely to underestimate the standard error, while MR-Egger, MR-Lasso, MR-Median, and MR-Mix constantly overestimate it. As for IVW, it underestimates the SE when the fraction is large and overestimates it when the fraction is small.
The coverage frequency refers to the frequency that the confidence interval covers the true causal effect among simulations. Here, this confidence interval is constructed by doubling $\hat{\rm se}(\hat\theta)$, which means that the coverage frequency corresponding to neither an inflated type-I error nor an inflated type-II error should be around 0.95. We observed that only MRBEE enjoys a coverage frequency around 0.95. When $m=250$, MR-Egger, MR-Lasso, and MR-Median suffer from inflated type-II error rates, likely because these methods cannot estimate the SE properly. These approaches also result in inflated-type I error rates caused by weak instrument bias as $m$ increases. Additionally, because MR-Mix overestimates the SE, it consistently exhibits a substantially inflated type-II error rate. Furthermore, IMRP and MR-ConMix consistently have inflated type I error rates because they frequently underestimate the SE.
We next verify if the asymptotic normal distributions in Theorem (ref) and Theorem (ref) are correct. For a general estimate $\hat\theta$, the asymptotic bias and SE are $\sqrt s_n(\hat\theta-\theta)$ and $\sqrt s_n$se($\hat\theta$), respectively, where $\sqrt s_n$ is the convergence rate of $\hat\theta$. If this estimate is strongly asymptotically unbiased, the asymptotic bias $s_n(\hat\theta-\theta)$ should also be 0. Besides, if two estimates have equal asymptotic SEs, they are equally powerful in terms of statistical efficiency. We select MRBEE, IVW, MR-Median, and MR-Lasso to compare, only consider two overlapping fractions: 100% and 0%, set $n_0=n_1=n_{\rm min}$, and fix the causal effect $\theta=0.5$. As for $m$ and $n_{\rm min}$, we focus on the following four cases:
Note that we directly generate the estimation errors $\mathbf W_\beta$ and $\boldsymbol w_\alpha$ according to Theorem (ref) because $n_{\rm min}$ in cases (3) and (4) can be larger than one million. The calculations involving individual-data are extremely time-consuming in these cases.
Figure (ref) demonstrates the simulation results. In case (1), $\hat\theta_{\rm BEE}$ is unbiased while the other three estimates suffer from non-removable biases. As for the asymptotic SE, $\sqrt{n^2_{\rm min}/m}~{\rm se}(\hat\theta_{\rm BEE})$ remains unchanged when $n_{\rm min}$ and $m$ are sufficiently large (e.g., the bars colored in blue), verifying conclusion $(iii)$ in Theorem (ref). However, the coverage frequency of MRBEE is a little larger than 0.95, meaning that the SE of $\hat\theta_{\rm BEE}$ is overestimated in this extreme case. This phenomenon is reasonable because Theorem (ref) points out that the convergence rate of the sandwich formula is $\min(\sqrt n_{\rm min}$,$\sqrt{n_{\rm min}^2/m},\sqrt{m/\log m})$, which slows down as $m$ increases. In case (2), the direct bias of $\hat\theta_{\rm IVW}$ is unchanged as $n_{\rm min}$ tends to infinity, confirming conclusion $(iii)$ in Theorem (ref). As for $\hat\theta_{\rm BEE}$, its asymptotic SE is a little larger than $\hat\theta_{\rm IVW}$, verifying item $(ii)$ in Theorem (ref).
In case (3), the asymptotic bias of $\hat\theta_{\rm IVW}$ is constant as $n_{\rm min}$ goes to infinity, illustrating that $\hat\theta_{\rm IVW}$ is not strongly asymptotically unbiased. As a result, the coverage frequencies of $\hat\theta_{\rm IVW}$ are significantly smaller than 0.95, confirming our claim that any inference made based on $\hat\theta_{\rm IVW}$ is invalid. Besides, the asymptotic SEs of $\hat\theta_{\rm BEE}$ and $\hat\theta_{\rm IVW}$ are essentially the same, indicating that $\hat\theta_{\rm BEE}$ and $\hat\theta_{\rm IVW}$ are equally efficient as long as $m/n_{\rm min}\to0$. In case (4), the asymptotic bias of IVW, MR-Median, and MR-Lasso vanish as $n_{\rm min}$ increases and their coverage frequencies are around 0.95, which is consistent with conclusion $(i)$ in Theorem (ref). The equal asymptotic SEs also indicate that $\hat\theta_{\rm BEE}$ and $\hat\theta_{\rm IVW}$ are equally efficient in this scenario. In addition, IVW, MR-Median, and MR-Lasso suffer from the same degree of bias when there is no pleiotropy, while MR-Median not only suffers from a large asymptotic SE but also is likely to overestimate it. To understand why MR-Median is always less efficient than IVW when there is no pleiotropy, its asymptotic behavior is worthy of future investigation.
For multivariable MR, we consider $p=6$ exposures and set the causal effect vector to be $\boldsymbol\theta=(0.3,0.3,-0.3,-0.3,0,0)^\top$. All of the exposures' IV-heritabilities are 0.3, while the outcome's IV-heritability is 0.15. We set an AR(1) structured genetic correlation matrix with coefficient $\rho=-0.5$ for the genetic effect $\boldsymbol\beta_j$, while considering a more intricate correlation structure for the noise terms $\boldsymbol u_i$ and $v_i$. In order to better mimic real data analysis, we take into account the scenario of completely overlapping GWAS samples (i.e., $n_{sk}=n_s=n_k$ for all $s,k$). Other cases of sample overlaps and details of the simulation settings are present in the supplementary materials.
Figure (ref) presents the comparison between the multivariable versions of IVW, MR-Egger, MR-Lasso, MR-Median, and MRBEE. In general, MRBEE is the only method that can produce unbiased causal effect estimates in all cases. As $m$ increases, the SE of $\hat{\boldsymbol\theta}_{\rm BEE}$ remains the same, while the estimation error of the SE estimate becomes smaller. However, a very large $m$ may conversely reduce the accuracy of the SE estimate in multivariable MR. For example, the SE estimates of all approaches in the cases of $m=1000$ have larger empirical variances than those in the cases of $m=500$. This phenomenon can be explained by Theorem (ref), which indicates that the convergence rate of the sandwich formula is $\min(\sqrt n_{\rm min}$,$\sqrt{n_{\rm min}^2/m},\sqrt{m/\log m})$. Hence, a larger $m$ may result in a worse SE estimate if $n_{\rm min}$ is not increased as $m$.
All the multivariable MR methods except MRBEE suffer from larger weak instrument biases with the increase of $m$. The SE estimates provided by these methods, in particular MR-Median, are less reliable than that of MRBEE. Thus, causal inferences based on the existing multivariable MR methods could be even more unreliable than univariable MR methods. In addition, $\hat{\boldsymbol\theta}_{\rm IVW}$ can have a bias toward any direction in multivariable MR. For example, the bias of $\hat\theta_{5,\rm IVW}$ is positive while the bias of $\hat\theta_{6,\rm IVW}$ is negative. The actual directions are jointly determined by the correlations of confounders and genetic effects.
We also examine the impact of omitting some important exposures. We conduct simulations when 1, 3, and all 6 exposures are included in the multivariable MR model, respectively. Figure (ref) illustrates the results of the simulations. We observed that if associated exposures are omitted, the causal effect estimates can suffer severe biases. The degree of the biases is jointly determined by the genetic covariance matrix and covariance matrix of confounders. In conclusion, even though MRBEE has eliminated the estimation error bias and weak instrument bias, OVB still exists if any relevant exposure is not specified in the multivariable MR model.
For univariable MR, we also investigated the effects of sample sizes, type-I error, winner's curse, and outlier detection. Regarding multivariable MR, we investigated the impact of different sample overlaps. In addition, the precision of estimating $\bm\Sigma_{W_\beta w_\alpha}$ by insignificant GWAS statistics is also studied. Only by increasing the sample sizes of the exposure and outcome cohorts simultaneously, the accuracy of MRBEE can be improved. The traditional MR methods suffer from inflated type-I errors when the overlapping fraction is large. After accounting for the weak instrument bias and estimate error bias, MRBEE is almost free of the winner curse's bias when the overlapping fraction is high. Furthermore, by applying the iterative method in IMRP, MRBEE can efficiently eliminate pleiotropic outliers and produce an accurate causal effect estimate. In addition, the estimation error $\bm\Sigma_{W_\beta w_\alpha}$ decreases with the increase of the number of insignificant variants $M$. Finally, multivariable MRBEE is accurate regardless of sample overlap. We summarized the findings with the simulation details in the supplementary material.
Cardiovascular disease including coronary artery disease (CAD) is one of the leading causes of death for both men and women worldwide. There are many epidemiological studies and MR analyses based on GWAS summary data dedicated to identifying the causal risk factors for CAD. However, the causal effects of the risk factors on CAD are less clear and the existing evidence can be contradictory. For example, elevated low-density lipoprotein cholesterol level (LDL-C) is a well-established causal risk factor for CAD scandinavian1994randomised, whereas wang2022mendelian concluded by multivariable MR analysis that LDL-C is not causally related to CAD in Europeans. Additionally, substantial observational analyses and molecular experiments have suggested that uric acid (UA) and red blood cell counts (RBC) contribute to the development of CAD bujak2015prognostic,yu2020uric. Nevertheless, wang2022mendelian did not observe significant causal effects of the two risk factors on CAD in Europeans. Furthermore, numerous MR analyses have concluded that body mass index (BMI) has a positive causal effect on CAD zhu2020mendelian,wang2022mendelian. However, recent literature indicates that BMI is likely to influence CAD through the mediation with diseases such as diabetes and hypertension gill2021risk. These contradictions may be due to biases in MR methods, including OVB, weak instrument bias, estimation error bias, etc.
We conducted two data analyses to estimate the causal effects of select risk factors on CAD. The first analysis uses the 11 exposures in wang2022mendelian, including BMI, hemoglobin (HB), hemoglobin a1c (Hba1c), hematocrit (HT), high-density lipoprotein cholesterol level (HDL-C), height, LDL-C, RBC, systolic blood pressure (SBP), triglycerides (TG), and UA. In wang2022mendelian, these 11 exposures were divided into two groups and analyzed separately. In contrast, we analyzed them in one multivariable MR model to avoid the OVB. In the second analysis, we replace HB, Hba1c, HT, and RBC with alcohol consumption (alcohol), diabetes, lifetime never smoking status (never.smoking), and sleeplessness. All the GWAS summary statistics used in our analyses were downloaded from the Neale lab (\url{http://www.nealelab.is/uk-biobank/}). Quality controls (QCs) are presented in the supplementary material. The total numbers of instrumental variants for the first and second analyses are 5345 and 5301, respectively.
Figure (ref) displays the causal effect estimates with 95% confidence intervals. MRBEE confirms the causal effects of LDL-C, RBC, and UA on CAD. Here, HB, HT, and RBC have high mutual correlations: $\widehat{\text{cor}}(x_\text{HB},x_{\rm HT})=0.89$, $\widehat{\text{cor}}(x_\text{HB},x_{\rm RBC})=0.63$, and $\widehat{\text{cor}}(x_\text{HT},x_{\rm RBC})=0.72$, and thus the inferences obtained by the existing methods are not reliable. For example, the existing MR methods suggest that RBC is not significant, HB has a significant positive effect, and HT has a significant negative effect, which contradicts the fact that HT and CAD are positively associated sorlie1981hematocrit. MRBEE corrects the estimation error bias and thus leads to a reasonable conclusion -- HB and RBC have positive causal effects on CAD while HT has a positive but insignificant causal effect on CAD. For the second analysis, MRBEE reveals that BMI is likely to affect CAD through the mediation of SBP and diabetes. In addition, MRBEE indicates that never.smoking is protective against CAD, whereas sleeplessness is associated with increasing CAD risk. Furthermore, due to the weak instrument bias and estimation error bias, the existing methods overestimate the effects of HDL-C and height and underestimate the causal effects of diabetes, LDL-C, never.smoking, SBP, and sleeplessness. By using MRBEE, we are able to obtain reliable causal effect estimates and therefore make valid inferences on the causal risk factors of CAD.
In this paper, we first investigated the asymptotic behavior of the multivariable IVW estimate. Since almost all MR methods are based on the IVW method, understanding the asymptotic behavior of the IVW estimate has very far-reaching implications for the theoretical and empirical studies of MR methods. We found that the bias of the multivariable IVW estimate is the product of weak instrument bias and estimation error bias. Also, we revealed that estimation error bias is a linear combination of measurement error bias and confounder bias, in which the sample overlaps trade off the proportion of these two components of estimation error bias. In the literature, although the phenomenon that the IVW estimate suffers from bias has been observed, a quantitative explanation for its existence is still absent. Our work fills the gap, which is a significant theoretical contribution to MR.
Subsequently, in this paper, we describe MRBEE that can yield the unbiased causal effect estimate $\hat{\boldsymbol\theta}_{\rm BEE}$. We point out that $\hat{\boldsymbol\theta}_{\rm BEE}$ is strongly asymptotically unbiased in all scenarios, indicating that $\hat{\boldsymbol\theta}_{\rm BEE}$ is asymptotically valid when making causal inferences. We also discuss how to perform MRBEE in practice, including how to estimate the bias-correction terms, how to estimate the sandwich formula, and how to identify possible UHP when multiple exposures are included. We present corresponding theorems to confirm that the estimates involved in the implementation of MRBEE are consistent in theory. In simulations, we show that MRBEE simultaneously estimates causal effects and the SE unbiasedly, and identifies UHP consistently. In section 5 and also in lorincz-comi2022mrbee, the practical advances of MRBEE are further demonstrated.
It is worth offering guidance on how to properly perform MR analysis from our perspective. First, we suggest applying the multivariable MR approach instead of the univariable MR approach because the causal effect estimates obtained by the univariable MR approach are unreliable due to OVB, regardless of the presence of UHP and CHP in the model. Second, rather than selecting the optimal number of instrumental variants such that the F statistics and conditional F statistics are larger than 10 burgess2011avoiding,sanderson2021testing, we advise including all the independent instrumental variants that are significantly associated with one or more exposures. Our theory illustrates that the asymptotic variance of a causal effect estimate is related to the cumulative variance explained by all specified IVs instead of the average variance explained by each IV. In particular, there is no need to worry about the issue of weak IVs because MRBEE has demonstrated efficiency to eliminate weak instrument bias through our simulations and theory. Third, when performing multivariable MR analysis, it is not necessary to remove variants that are pleiotropic between the exposures. For example, wang2022mendelian observed that LDL-C was insignificantly associated with CAD in Europeans, which is unlikely to be true because this risk causality has been well established in randomized clinical trials scandinavian1994randomised. The potential reason for this false negative is that wang2022mendelian excluded the IVs associated with RBC, HB, HT, and UA in their multivariable MR analysis. We believe that the proper way to perform multivariable MR analysis is to simultaneously include all the relevant exposures, as the multivariable regression can automatically account for the pleiotropic variants shared by the specified exposures. Fourth, among the existing multivariable MR approaches including IVW, MR-Egger, and MR-Lasso, we recommend MRBEE as the primary analysis approach because it has been proven to be the only one that enjoys strongly asymptotic unbiasedness in the presence of many weak IVs.