key: cord-0900446-f17d5fmi
authors: Xu, S.; Fung, W. K.; Liu, Z.
title: pIVW: A novel Mendelian Randomization Method Accounting for Weak Instruments and Horizontal Pleiotropy with Applications to the COVID-19 Outcomes
date: 2021-09-27
journal: nan
DOI: 10.1101/2021.09.25.21264115
sha: 650e4e9b0a219e6eb1657bba9673ea5078632792
doc_id: 900446
cord_uid: f17d5fmi

Mendelian randomization (MR) utilizes genetic variants as instrumental variables (IVs) to estimate the causal effect of an exposure variable on an outcome of interest even in the presence of unmeasured confounders. However, many MR methods including the most popular inverse-variance weighted (IVW) estimator could be biased by the weak IVs that are weakly associated with the exposure. In this article, we develop a novel method called penalized inverse-variance weighted (pIVW) estimator, where we adjust the IVW estimator to account for the weak IVs by a proposed penalization method to prevent the denominator of the pIVW estimator from being close to zero. Moreover, we account for the horizontal pleiotropy|a widespread phenomenon in human genome that could bias the inference for the causal effect|by adjusting the variance estimation of the pIVW estimator. The proposed pIVW estimator can reduce to the debiased IVW (dIVW) estimator|another extension of the the IVW estimator|when the number of IVs and the IV strength increase. More generally, we prove that the pIVW estimator can achieve smaller bias and variance than the dIVW estimator under some regularity conditions. We also illustrate the improved performance of the proposed pIVW estimator over competing MR methods through a comprehensive simulation study. Further, we analyze the causal effects of the obesity-related traits and diseases on the Coronavirus disease 2019 (COVID-19). Notably, we find that hypertensive disease is associated with increased risk of hospitalized COVID-19, while peripheral vascular disease and higher body mass index are associated with increased risks of COVID-19 infection, hospitalized COVID-19 and critically ill COVID-19. The R package for the pIVW method is publicly available at https://github.com/siqixu/mr.pivw.

It is of major interest in health studies to identify causal risk factors associated with various clinical outcomes. For example, identifying the causal risk factors of the Coronavirus disease 2019 is currently one of the most pressing global public health problems (Jordan et al., 2020; Zheng et al., 2020) . The COVID-19 pandemic, caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has posed a serious threat to human health all over the world (Pascarella et al., 2020) . As of July 1, 2021, the COVID-19 pandemic has led to over 180 million cases and more than 3.9 million deaths, which has brought unprecedented medical and economic burdens worldwide . To reduce the risks of COVID-19 incidence and mortality, it is crucial to identify the causal risk factors for the development of public health policies and clinical strategies for prevention and intervention. So far, associations between some exposure variables such as obesity and the COVID-19 outcomes have been reported by several epidemiological studies (Popkin et al., 2020; Stefan et al., 2020) . However, the associations identified from the observational data might be subject to unmeasured confounding of the exposure-outcome relationship.

To address the unmeasured confounding issue in observational studies, Mendelian randomization (MR) utilizes genetic variants as instrumental variables (IVs) to estimate the causal effect of an exposure variable on an outcome of interest even in the presence of unmeasured confounders Ebrahim, 2003, 2004; Sheehan et al., 2008) . Due to the availability of summary-level data from the genome-wide association studies (GWASs), a plenty of methods were developed for GWAS summary-level data using multiple genetic variants as IVs (Zheng et al., 2017; Lawlor, 2016) , which facilitate the wide use of MR analysis in health studies. However, the validity of MR strictly depends on the following three core assumptions defining a valid IV (Lawlor et al., 2008; Didelez and Sheehan, 2007) :

• IV Relevance: the IV is associated with the exposure;

• IV Independence: the IV is independent of any confounders of the exposure-outcome relationship;

• Exclusion Restriction: the IV only affects the outcome via the exposure.

When either of these IV assumptions is violated, the MR analysis may also yield biased estimation of the causal effect. In particular, the first IV assumption (IV Relevance) can be nearly violated when the IVs are only weakly associated with the exposure variable Davies et al., 2015) . In MR studies, the weak IV bias may occur when the genetic variants only explain a small proportion of variance for the exposure variable. On the other hand, the widespread horizontal pleiotropy in human genome can also lead to the violation of the third IV assumption (Exclusion Restriction) (Verbanck et al., 2018; Hemani et al., 2018) , which is a phenomenon that the genetic variants directly affect the outcome not mediated by the exposure variable (see Figure 1 for a graphical illustration).

The inverse-variance weighted (IVW) estimator (Burgess et al., 2013) is the most popular MR method being widely used in health studies. It has a simple and explicit expression, which combines the estimated causal effects from multiple IVs into a weighted average with the idea borrowing from the fixed-effect meta-analysis literature. Despite its widespread popularity, studies Ye et al., 2019) pointed out that the IVW estimator can be seriously biased by the weak IVs. Besides the IVW estimator, many common MR methods were found to be subject to the weak IV bias, such as MR-Egger (Bowden et al., estimator, which was shown to be robust to weak IVs. However, in contrast to the IVW estimator, MR-RAPS has no closed form of solution and might have multiple roots. Recently, the debiased IVW (dIVW) estimator (Ye et al., 2019) was proposed to account for the weak IV by a simple modification to the IVW estimator. Comparing to the IVW estimator, the dIVW estimator was proved to be consistent under weaker conditions that allow the presence of many weak IVs. Nevertheless, as a ratio estimator, the dIVW estimator is likely to yield an inflated estimate when its denominator is close to zero, which might occur in the presence of many weak IVs. Moreover, when the denominator is close to zero, a ratios estimator may have heavy tailed distributions and may not possess finite moments, such as the ratio of two normal random variables (Piegorsch and Casella, 1985; Marsaglia et al., 2006; Press, 1969) .

In the economic literature, many IV methods as ratio estimators also encounter similar issues especially under the the weak IV setting (Nelson and Startz, 1990; Andrews et al., 2019) .

In this article, we develop a novel method called penalized inverse-variance weighted (pIVW) estimator, where the IVW estimator is further adjusted by a proposed penalized log-likelihood function. Through the penalization, we prevent the denominator in the ratio estimator from being close to zero and thus provide improved causal estimation. Moreover, we account for the horizontal pleiotropy by adjusting the variance estimation of the pIVW estimator. The proposed pIVW estimator have some nice features. First, our theoretical and numerical results show that the proposed pIVW estimator can achieve smaller bias and variance than the dIVW estimator under some regularity conditions. Second, it is consistent and asymptotically normal even in the presence of many weak IVs, and requires no further assumptions than the dIVW estimator. Third, it has a unique and explicit form of solution, whereas some other robust MR methods like MR-RAPS might have multiple roots and diverging solutions. We also illustrate the improved performance of the proposed pIVW estimator over the competing MR methods by a comprehensive simulation study.

Furthermore, we focus our interest on the causal effects of the obesity-related exposures on the COVID-19 outcomes. Some recent epidemiological studies have found that obesity and some obesity-related diseases (e.g., hypertension) were associated with increased risks of COVID-19 infection and severity (Popkin et al., 2020; Nakeshbandi et al., 2020; Klang et al., 2020; Zhang et al., 2020; de Almeida-Pititto et al., 2020) , however, their associations might arise from the unobserved confounders rather than the causality. Although some MR studies have also found associations between body mass index (BMI) on COVID-19 outcomes (Leong et al., 2021; Ponsford et al., 2020; Aung et al., 2020; Freuer et al., 2021) , but the associations between the obesity-related diseases and COVID-19 outcomes have seldom been studied. Moreover, current MR analyses mostly relied on the traditional MR methods such as the IVW estimator that may suffer from the weak IV bias. In particular, the weak IVs are likely to exist when the exposure variables are some complex traits and diseases. Therefore, we aim to provide more robust estimation by utilizing the proposed pIVW estimator to account for the weak IVs. Specifically, we apply the pIVW estimator to analyzing the causal effects of five obesity-related traits and diseases (i.e., peripheral vascular disease, dyslipidemia, hypertensive disease, type 2 diabetes and BMI) on three COVID-19 outcomes (i.e., COVID-19 infection, hospitalized COVID-19 and critically ill COVID-19).

Among them, we find that hypertensive disease is significantly associated with increased risk of hospitalized COVID-19, while peripheral vascular disease and higher BMI are significantly associated with increased risks of COVID-19 infection, hospitalized COVID-19 and critically ill COVID-19.

Suppose that there are p independent genetic variants {G j } p j=1 . When there exists no horizontal pleiotropy, the relationships among the genetic variants G j s, the exposure X, the outcome Y and the unmeasured confounder U depicted in Figure 1 can be formulated by the linear structural models as follows (Bowden et al., 2015) :

where γ j is the genetic effect of G j on X, β is the causal effect of our interest, and X and Y are mutually independent random errors. Let Γ j denotes the effect of G j on Y , then we 5 . CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 have Γ j = βγ j by substituting Equation (1) for X in Equation (2). When Y is a binary outcome and Equation (2) is replaced by a logistic model, Γ j = βγ j can still be justified by a probit approximation of the logistic model Vansteelandt et al., 2011) .

Letγ j andΓ j be the estimates of γ j and Γ j with the variances σ 2 γ j and σ 2 Γ j , respectively.

In the two-sample MR design (Lawlor, 2016) , γ j , σγ j p j=1 and Γ j , σΓ j p j=1 can be obtained from two independent GWASs. Since the GWASs generally involve large population, it is common to assume thatγ j andΓ j are independently distributed asγ j ∼ N γ j , σ 2 Figure 1 : The relationships among the jth genetic variant G j , the exposure X, the outcome Y and the unmeasured confounder U . The effect of G j on X is γ j , the direct effect (pleiotropic effect) of G j on Y is α j , and the causal effect of X on Y is β.

The popular inverse-variance weighted (IVW) estimator combines the estimated causal effectsβ j =Γ j /γ j from each genetic variant with the weights w j = σ −2 Γ jγ 2 j as follows:

We have β = µ 1 /µ 2 since Γ j = βγ j under models (1)-(2). Then, as shown by Zhao et al. (2020) , the IVW estimator can be approximated byβ

When there is no measurement error for γ j (i.e., σ 2 γ j = 0), we haveβ IV W ≈ µ 1 /µ 2 = β.

. CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 However, studies (Ye et al., 2019; Zhao et al., 2020) indicated that the IVW estimator can be seriously biased toward zero for ignoring the measurement errors of γ j especially in the presence of many weak IVs that have small σ −2 γ j γ 2 j . To handle the bias due to weak IVs, the debiased IVW (dIVW) estimator replaces the denominator in the IVW estimator with an unbiased estimatorμ 2 of µ 2 aŝ

The dIVW estimator was shown to be consistent and asymptotically normal under weaker conditions than the IVW estimator that allow the presence of weak IVs. However, the dIVW estimator is more likely to yield extreme estimates in the presence of weak IVs, because its denominatorμ 2 has larger coefficient of variation and thus higher probability of being close to zero than that of the IVW estimator. Therefore, in the following Section 3, we extend the IVW estimator to account for the weak IVs by proposing a penalized log-likelihood function for µ 1 and µ 2 , which can prevent the estimator of µ 2 from being close to zero and then provide improved causal estimation.

Assume that the estimatorsμ 1 = p j=1 σ −2

We propose a penalized log-likelihood function to adjust the estimates of µ 1 and µ 2 as

where f (μ 1 ,μ 2 ) denotes the bivariate normal density function of (μ 1 ,μ 2 ) and λ > 0 is the penalty parameter. We adopt the penalty term log |µ 2 | inspiring by Wang et al. (2020) , which was also used to handle the issue of near-zero denominator when deriving the confidence interval for a ratio estimator. Through the penalty, the penalized log-likelihood l p (µ 1 , µ 2 ) will become smaller when µ 2 approaches to zero. Therefore, an estimate of µ 2 being close to zero is less preferable by l p (µ 1 , µ 2 ). Specifically, by maximizing l p (µ 1 , µ 2 ), we can obtain the adjusted estimators of µ 1 and µ 2 as

are the estimates of v 12 and v 2 , respectively. Then, we propose the penalized IVW (pIVW) estimator as a ratio ofμ 1 andμ 2 , that iŝ

Note that when the penalty parameter λ = 0, we haveμ 2 =μ 2 and thereforeβ pIV W reduces toβ dIV W . When λ > 0,β pIV W will also be close toβ dIV W whenv 2 /μ 2 2 approaches to zero. To see this, we define the IV strength as

Further, we make the following assumptions:

Assumption 1. The number of IVs, p, diverges to infinity.

away from zero and infinity for all j. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10. 1101 We notice that Assumptions 1 and 2 are also required for the consistency of the dIVW estimator (Ye et al., 2019) . They are reasonable in the two-sample MR design since the GWASs often involve large population and a large amount of genetic variants. And the independence among γ j ,Γ j p j=1 can be guaranteed by the linkage-disequilibrium clumping.

Then, under Assumptions 1 and 2,v 2 /μ 2 2 = O p 1 κp + 1 κ 2 p and henceβ pIV W reduces tô

More importantly, in the following Theorem 3.1 (a)-(b), we show thatβ pIV W can achieve smaller bias and variance thanβ dIV W under proper choice of the penalty parameter λ. Further, in the Theorem 3.1 (c), we show thatβ pIV W is consistent and asymptotically normal under some regularity conditions.

Theorem 3.1. Suppose that models (1)-(2) and Assumptions 1-2 hold. Then, as κ √ p → ∞,

we have the following conclusions.

(a) The approximate bias ofβ dIV W is of order

Then,β pIV W is consistent and asymptotically normally distributed as

The consistency and asymptotic normality ofβ pIV W still hold, when we replace V with the

Under the setting of many weak IVs, we may have κ → 0 as p → ∞ because more weak

IVs are likely to be included into the analysis as the number IVs p increases, which may reduce the IV strength κ. Nevertheless, the above theorem can still hold in this case as long as κ √ p → ∞, which means that it allows the presence of many weak IVs. On the other hand, although many ratio estimators may not have finite moments when their denominators 9 . CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10. 1101 are close to zero (Piegorsch and Casella, 1985; Marsaglia et al., 2006; Press, 1969) , for the denominatorμ 2 ofβ dIV W , we have P (|μ 2 | > ) → 1 for any > 0 as κ √ p → ∞. And we have same result for the denominatorμ 2 ofβ pIV W sinceμ 2 is in larger magnitude than µ 2 due to Equation (6). Therefore, we can still approximate the biases and variances of

that the approximate bias ofβ pIV W converges to zero in a faster rate than that ofβ dIV W when λ = 1, while the approximate variance ofβ pIV W is also smaller than that ofβ dIV W , which implies that we may choose λ = 1 for the pIVW estimator in practice to achieve improved performances over the dIVW estimator. Further, in the Theorem 3.1 (c), we show thatβ pIV W is consistent and asymptotically normal, which requires no further assumptions comparing toβ dIV W .

In the following sections, we show that these nice properties of the proposed pIVW estimator still hold when we conduct selection to exclude some weak IVs from the analysis and account for the horizontal pleiotropy.

To handle the weak IV bias, it is common in MR studies to exclude the weak IVs with small σ −2 γ jγ 2 j from the analysis. However, the IV selection based on the exposure dataset γ j , σγ j p j=1 might lead to selection bias or winner's curse that would also bias the causal estimation . To address this problem, in the three-sample MR design (Zhao et al., 2019) , IV selection is performed on a third dataset γ * j , σ * γ j p j=1 which is independent of the exposure dataset and the outcome dataset. Specifically, an IV is included into the analysis when |γ * j | > δσ * γ j with a pre-set threshold δ > 0. Ye et al. (2019) showed that IV selection with an appropriate threshold δ could reduce the bias of the IVW estimator and improve the efficiency of the dIVW estimator. They also recommended a threshold δ = √ 2 log p to guarantee small probability of selecting any null IVs (i.e., γ j = 0).

When IV selection is performed at a threshold δ, we define the IV strength as

. CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ;  where q δ,j = P |γ * j | > δσ * γ j and p δ = p j=1 q δ,j . Let S δ = j : |γ * j | > δσ * γ j be the set of selected IVs. Then, we can estimate κ δ byκ δ =p −1 δ j∈S δ γ 2 j − σ 2 γ j σ −2 γ j , wherep δ is the number of selected IVs within S δ . To study the theoretical properties of the proposed pIVW estimator under IV selection, we have the following Assumption 3 for the summary-level data in the selection dataset.

are mutually independent andγ * j ∼ N γ j , σ * 2 γ j with known variance σ * 2 γ j for every j. The ratio of variances σ 2 γ j /σ * 2 γ j is bounded away from zero and infinity for all j.

Given a selection threshold δ, we evaluate the dIVW estimatorβ δ,dIV W =μ 1,δ /μ 2,δ and the proposed pIVW estimatorβ δ,pIV W =μ 1,δ /μ 2,δ using the selected IVs within the set S δ . In the following Theorems 3.2, we show that,β δ,pIV W can still have smaller bias when λ = 1 and has smaller variance when λ > 0 thanβ δ,dIV W , andβ δ,pIV W is consistent and asymptotically normal under some regularity conditions.

Theorem 3.2. Suppose that models (1)-(2) and Assumptions 1-3 hold, and κ δ √ p δ / max(1, δ 2 ) → ∞. Then, we have the following conclusions.

(a) The approximate bias ofβ δ,dIV W is of order

The approximate variance ofβ δ,pIV W is smaller than that ofβ δ,dIV W at the order O 1 κ δ p δ + 1+δ 4 κ 2 δ p δ when λ > 0.

(c) Further assume that max j γ 2 j σ −2 γ j q δ,j /(κ δ p δ + p δ ) → 0. Then,β δ,pIV W is consistent and asymptotically normally distributed as

The consistency and asymptotic normality ofβ δ,pIV W still hold, when we replace V δ with the 11 . CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint 

Under the IV selection at a threshold δ, Theorem 3.2 provides similar conclusions to Theorems 3.1. It indicates that the pIVW estimator can have smaller bias and variance than the dIVW estimator, and it is consistent and asymptotically normal regardless of whether or not the IV selection is conducted. In comparison, Theorem 3.2 also involves the IV selection threshold δ which may vary with the number of IV p, for instance, δ = √ 2 log p. In the special case when δ = 0 (i.e., no IV selection), Theorems 3.2 is equivalent to Theorems 3.1.

When there exists horizontal pleiotropy (i.e., direct effect of G j on Y not mediated by X), the linear structural model (2) can be modified as

where α j denotes the direct genetic effect of G j on the outcome Y (i.e., pleiotropic effect).

In this case, we have Γ j = βγ j + α j . We follow a common practice in many MR methods to assume that the horizontal pleiotropy is balanced (i.e., the pleiotropic effect has mean zero) (Ye et al., 2019; Zhao et al., 2020; Bowden et al., 2017; Cheng et al., 2020) and treat α j as random effect following α j ∼ N (0, τ 2 ). Then,Γ j ∼ N βγ j , σ 2 Γ j + τ 2 has larger variance in the presence of horizontal pleiotropy due to the non-zero τ 2 . To account for the horizontal pleiotropy, we adjust the variance estimation ofβ δ,pIV W by 

12 . CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint 

To establish the theoretical results for the bias, the variance and the asymptotic property

We further assume that τ 2 < c 1 σ 2 Γ j with a constant c 1 for all j, and that max j σ −2

for a constant c 2 in Theorem 3.2 (c). Then, Theorem 3.2 still hold with V δ andV δ replaced by V * δ andV * δ , respectively.

We randomly generate the summary-level data for 1000 IVs fromγ j ∼ N γ j , σ 2

independently. For the true γ j , we consider three scenarios as follows.

• Scenario A (some strong IVs and many null IVs): 10 IVs have γ j ∼ N (0, 0.03 2 ) and the rest of IVs have γ j = 0. As such, the IV strength of 10 non-null IVs is around 15.54 and that of all 1000 IVs is around 0.16.

• Scenario B (many weak IVs and many null IVs): 100 IVs have γ j ∼ N (0, 0.02 2 ) and the rest of IVs have γ j = 0. As such, the IV strength of 100 non-null IVs is around 3.11 and that of all 1000 IVs is around 0.31.

• Scenario C (all weak IVs):all IVs have γ j ∼ N (0, 0.02 2 ). The IV strength is around 4.07.

Then, we let Γ j = βγ j + α j , where α j ∼ N (0, τ 2 ). We set β = 0.5, and τ = 0 and 0.01 representing the absence and the presence of horizontal pleiotropy, respectively. For the variances σ 2 γ j and σ 2 Γ j , we calculate them by σ 2

where n X and n Y denote the sample sizes of the GWASs for the exposure and the outcome, respectively. We set n X = n Y = 100, 000. For var(G j ), we let G j ∼ Bin(2, M AF j ) and randomly generate the minor allele frequencies by M AF j ∼ U (0.1, 0.5). For var(X) and var(Y ),

. CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint (1) and (6) with the variances of U , X and Y being 2, respectively. Furthermore, we generate an independent dataset withγ * j ∼ N γ j , σ 2 γ j /2 for the IV selection at threshold δ = √ 2 log p = 3.72. The simulation is based on 10,000

replicates.

In the simulation study, we first investigate the impact of the penalty parameter λ on the The performances among various methods are compared in terms of the relative bias (bias divided by the true β) and the empirical standard error of the estimated causal effect, as well as the coverage probability of 95% confidence interval.

In Table 1 , we investigate the performances of the pIVW estimator with various penalty parameter λ when no horizontal pleiotropy exists and no IV selection is conducted. Note that, when λ = 0, the pIVW estimator reduces to the dIVW estimator. The pIVW estimator has smallest bias at λ = 1 under Scenarios A and B. Specifically, the bias decreases as λ increases from 0 to 1 and increases as λ increases from 1 to 2.5. On the other hand, the empirical standard error decreases as the increase of λ. The coverage probability of confidence interval is generally around 95%. In Scenario C where there is larger IV strength than Scenarios A and B, the pIVW estimator has nearly the same performance when λ varies from 0 to 2.5. We observe similar impact of λ on the performances of the pIVW estimator when there exists horizontal pleiotropy and when the IVs are selected from an independent dataset (Supplementary Tables S1.1-S1.3).

In Table 2 , we compare the proposed pIVW estimator with the competing MR methods when no horizontal pleiotropy exists and no IV selection is conducted. Here, we use the penalty parameter λ = 1 for the pIVW estimator. In Scenarios A and B, the IVW,

14 . CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; Table 1 : The performances of the pIVW estimator with various penalty parameter λ. The relative bias (bias divided by the true causal effect), the empirical standard error and the coverage probability of the 95% confidence interval. The true causal effect β = 0.5. No horizontal pleiotropy exists (τ = 0). No IV selection is conducted. The simulation is based on 10,000 replicates In comparison, the MR-RAPS and dIVW estimators have smaller biases, but we find that they are likely to produce extreme estimates which result in their relatively large empirical standard errors in Scenarios A (see Supplementary Figure S1 .1 for the box plot of the estimates from 10,000 simulation replicates). In contrast, the proposed pIVW estimator has smallest bias among six methods and also has smaller empirical standard errors than the dIVW estimators. In Scenario C where the IV strength is larger, all the methods have improved performances. However, the IVW, MR-Egger and MR-Median estimators still have substantial biases, whereas the proposed pIVW estimator is nearly unbiased. The MR-RAPS and dIVW estimators perform similarly to the pIVW estimator in this case. In addition,

we study the impact of the number of IVs (i.e., p) by increasing p from 1000 to 5000 in all three scenarios and using 50 and 500 non-null IVs in Scenarios A and B, respectively (Supplementary Tables S1.4). We find that the increase in p has little impact on the biases 15 . CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101/2021.09.25.21264115 doi: medRxiv preprint Table 3 : Comparison among six methods. The relative bias (bias divided by the true causal effect), the empirical standard error and the coverage probability of the 95% confidence interval. The true causal effect β = 0.5. No horizontal pleiotropy exists (τ = 0). The IV selection threshold δ = √ 2 log p. The simulation is based on 10,000 replicates (Supplementary Tables S1.5-S1.6), we find similar results to those in Tables 2-3 except that all the methods have increased standard errors in these cases.

On the other hand, although the IVW estimator suffers from substantial bias in the presence of weak IVs, Ye et al. (2019) indicated that the IVW estimator was still consistent when the the true causal effect β was zero. Therefore, we conduct further simulation under Table S1 .7). In this case, we find that all six methods are nearly unbiased even in the the presence of weak IVs under Scenarios A-C. However, the MR-RAPS and dIVW estimators still have some extremely large estimates and relatively large empirical standard errors in Scenario A (Supplementary Figure S1.2) . Overall, the proposed pIVW estimator is applicable to more general situations than the competing methods, which has smallest bias among all six methods and has smaller standard error than the dIVW estimator 17 . CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 under the situations considered.

In this section, we focus our interest on the causal effects of five obesity-related diseases and traits (i.e., peripheral vascular disease, dyslipidemia, hypertensive disease, type 2 diabetes and BMI) on COVID-19. Specifically, we consider three categories of COVID-19 outcomes including (1) (Banda et al., 2015) with 53,991 individuals of European ancestry and UK BioBank (Abbott et al., 2018) with 108,039 individuals of European ancestry, respectively. More detailed data description is provided in the Supplementary Table S2 .1. To exclude correlated IVs, we perform the linkage-disequilibrium clumping to remove the correlated genetic variants within 10Mb pairs and with the linkage disequilibrium r 2 < 0.001. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 . No IV selection is conducted. In each cell, the first row displays the estimated causal effect and the estimated standard error (in bracket), and the second row displays the P value which is labelled with " * " when less than 0.05 significance level. The pIVW estimator with penalty parameter λ = 1 19 . CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ;  significance level, which are concordant with some recent epidemiological studies indicating that hypertensive disease and higher BMI were significantly associated with increased risks of the COVID-19 outcomes (Popkin et al., 2020; Nakeshbandi et al., 2020; Klang et al., 2020; Zhang et al., 2020; de Almeida-Pititto et al., 2020) . Some previous MR studies have also found significant causal effects of BMI on the COVID-19 outcomes (Leong et al., 2021; Aung et al., 2020; Freuer et al., 2021; Ponsford et al., 2020) , but there is still no MR analysis on the hypertensive disease to the best of our knowledge. Additionally, the pIVW estimator suggests that peripheral vascular disease is significantly associated with higher risks of three COVID-19 outcomes under the IV selection at threshold δ = √ 2 log p (see Supplementary Figure S2 .1), which may improve the efficiency of the pIVW estimator. To our knowledge, there also lack epidemiological studies and MR studies on the associations between peripheral vascular disease and the COVID-19 outcomes, despite a high incidence of peripheral vascular disease as complication in COVID-19 patients (García-Ortega et al., 2021; Hanff et al., 2020; Piazza and Morrow, 2020) . For type 2 diabetes and dyslipidemia, the pIVW estimator gives no evidence of their associations with the COVID-19 outcomes regardless of the IV selection, while their associations with COVID-19 remain unclear due to discordant findings among previous epidemiological studies and MR studies (Ponsford et al., 2020; Leong et al., 2021; Zheng et al., 2020; Hariyanto and Kurniawan, 2020; Yang et al., 2021; Aung et al., 2020) .

Like the pIVW estimator, the competing MR methods also have significant findings for the causal effects of peripheral vascular disease, hypertension disease and BMI on the is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ;  dIVW estimator is very close to zero in the presence many weak IVs. For dyslipidemia (κ = 0.82) and type 2 diabetes (κ = 0.31), all methods have insignificant findings and have similar estimates even in the presence many weak IVs. It is likely because the competing MR methods could be nearly unbiased when the true causal effect is zero as shown in the simulation study. is smaller discrepancy among the methods. The pIVW estimator also reduces to the dIVW estimator in this case and therefore they provide nearly the same results.

In MR analysis, the estimation of causal effect can be biased by the presence of weak IVs.

The IVW estimator is one of the most popular MR methods but also suffers from substantial bias due to the weak IVs. In this paper, we develop a penalized IVW (pIVW) estimator, where the IVW estimator is adjusted by a proposed penalization method which prevents the denominator of the ratio from being zero in order to decrease the bias and the variance of the estimator in the presence of many weak IVs. Moreover, we allow for the horizontal pleiotropy by adjusting the variance estimation of the pIVW estimator.

The proposed pIVW estimator have some nice features. First, our theoretical and numerical results show that the proposed pIVW estimator can achieve smaller bias and variance compared to the dIVW estimator-another recently proposed extension of the IVW estimator-in the presence of many weak IVs. Second, it is consistent and asymptotically normal even in the presence of many weak IVs, and requires no further assumptions than the dIVW estimator. Third, it has a unique and explicit form of solution, whereas many competing MR methods that are robust to the weak IVs may not have a closed form of solution and may have multiple diverging roots. We also illustrate the superior performance 21 . CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101/2021.09.25.21264115 doi: medRxiv preprint of the pIVW estimator over other popular MR methods through comprehensive simulation study. Our theoretical and numerical results also provide guidance for the choice of the penalty parameter λ in the pIVW estimator, which plays a role in the trade-off between the bias and the variance. We show that the pIVW estimator has smallest bias when λ = 1 and generally has decreased variance as the increase of λ. Thus, to guarantee the pIVW estimator to achieve smallest bias and valid inference, we recommend to choose the penalty parameter to be 1 in practice. Note that, the dIVW estimator is also a special case of the pIVW estimator when λ = 0. As the number of IVs and the IV strength increase, the choice of λ tends to have less influence on the performance of the pIVW estimator, and the pIVW estimator can also reduce to the dIVW estimator in this case. Therefore, our pIVW estimator is applicable to more general situations and is recommended for practical use.

Further, we analyze the causal effects of five obesity-related diseases and traits (i.e., peripheral vascular disease, dyslipidemia, hypertensive disease, type 2 diabetes and BMI) on three COVID-19 outcomes (i.e., COVID-19 infection, hospitalized COVID-19 and critically ill COVID-19). We find significant positive causal effects of hypertensive disease on hospitalized COVID-19 and BMI on three COVID-19 outcomes by the pIVW estimator at 0.05 significance level, which are accordant with the findings from some recent epidemiological studies. Additionally, we find significant positive causal effects of peripheral vascular disease on three COVID-19 outcomes under IV selection, while there still lack epidemiological and MR studies on the association between peripheral vascular disease and COVID-19 despite its high incidence in COVID-19 patients. And we find no evidence supporting the associations of type 2 diabetes and dyslipidemia with the COVID-19 outcomes, while their roles on the COVID-19 outcomes also remain unclear due to discordant findings among previous epidemiological and MR studies. Overall, our findings suggest that the preventive strategies aiming at obesity and some related diseases, such as control of body weight and prevention of hypertensive disease and peripheral vascular disease, may help to reduce the risks of COVID-19 infection and severity. On the other hand, comparing to other MR estimators, the proposed pIVW estimator may have smaller bias in the presence of many weak IVs.

For instance, the peripheral vascular disease has small IV strengths implying the existence of many weak or null IVs. In this case, the IVW estimator has an estimate very close to zero, but which might be underestimated as shown in the simulation studies. In contrast,

. CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 the dIVW estimator yields a very large estimate with an extreme standard error, but which might also be inflated because its denominator could be close to zero when many weak IVs exist. The causal effect estimate from the proposed pIVW estimator may be more reliable and have smaller bias, which lies between those of the IVW and dIVW estimators and has much smaller standard error than that of the dIVW estimator.

In future study, it would also be possible to extend the proposed penalization method in the pIVW estimator to more flexible model frameworks for handling the weak IVs and other problems in MR simultaneously, such as the extension to MR-Egger for the unbalanced pleiotropy.

. CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101/2021.09.25.21264115 doi: medRxiv preprint . CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101/2021.09.25.21264115 doi: medRxiv preprint . CC-BY-NC 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 

Round 2 gwas results of thousands of phenotypes in the uk biobank

Genome-wide association study identifies 112 new loci for body mass index in the japanese population

Weak instruments in instrumental variables regression: Theory and practice

Causal inference for genetic obesity, cardiometabolic profile and covid-19 susceptibility: A mendelian randomization study

Characterizing race/ethnicity and genetic ancestry for 100,000 subjects in the genetic epidemiology research on adult health and aging (gera) cohort

A framework for the investigation of pleiotropy in two-sample summary data mendelian randomization

Mendelian randomization with invalid instruments: effect estimation and bias detection through egger regression

Consistent estimation in mendelian randomization with some invalid instruments using a weighted median estimator

Mendelian randomization analysis with multiple genetic variants using summarized data

Bias in causal estimates from mendelian randomization studies with weak instruments

Avoiding bias from weak instruments in mendelian randomization studies

Mrldp: a two-sample mendelian randomization for gwas summary statistics accounting for linkage disequilibrium and horizontal pleiotropy

The covid-19 host genetics initiative, a global initiative to elucidate the role of host genetic factors in susceptibility and severity of the sars-cov-2 virus pandemic

Mapping the human genetic architecture of covid-19 by worldwide meta-analysis

The many weak instruments problem and mendelian randomization

Severity and mortality of covid 19 in patients with diabetes, hypertension and cardiovascular disease: A meta-analysis

Mendelian randomization as an instrumental variable approach to causal inference

An interactive web-based dashboard to track covid-19 in real time. The Lancet infectious diseases

Impact of body composition on covid-19 susceptibility and severity: A two-sample multivariable mendelian randomization study

Incidence, risk factors, and thrombotic load of pulmonary embolism in patients hospitalized for covid-19 infection

Thrombosis in covid-19

Dyslipidemia is associated with severe coronavirus disease 2019 (covid-19) infection

Evaluating the potential role of pleiotropy in mendelian randomization studies

Covid-19: risk factors for severe disease and death

Severe obesity as an independent risk factor for covid-19 mortality in hospitalized patients younger than 50

Commentary: Two-sample mendelian randomization: opportunities and challenges

Mendelian randomization: using genes as instruments for making causal inferences in epidemiology

Cardiometabolic risk factors for covid-19 susceptibility and severity: A mendelian randomization analysis

Ratios of normal variables

The impact of obesity on covid-19 complications: a retrospective cohort study

The distribution of the instrumental variables estimator and its t-ratio when the instrument is a poor one

Covid-19 diagnosis and management: a comprehensive review

Diagnosis, management, and pathophysiology of arterial and venous thrombosis in covid-19

The existence of the first negative moment

Cardiometabolic traits, sepsis, and severe covid-19: A mendelian randomization investigation

Individuals with obesity and covid-19: A global perspective on the epidemiology and biological relationships

The t-ratio distribution

Mendelian randomisation and causal inference in observational epidemiology

mendelian randomization': can genetic epidemiology contribute to understanding environmental determinants of disease?

Mendelian randomization: prospects, potentials, and limitations

Obesity and impaired metabolic health in patients with covid-19

On instrumental variables estimation of causal odds ratios

Detection of widespread horizontal pleiotropy in causal relationships inferred from mendelian randomization between complex traits and diseases

Penalized fieller's confidence interval for the ratio of bivariate normal means

Lack of significant association between dyslipidemia and covid-19 mortality

Debiased inverse-variance weighted estimator in two-sample summary-data mendelian randomization

Association of hypertension with the severity and fatality of sars-cov-2 infection: A meta-analysis

Powerful three-sample genomewide design and robust statistical inference in summary-data mendelian randomization

Statistical inference in two-sample summary-data mendelian randomization using robust adjusted profile score

Recent developments in mendelian randomization studies

Risk factors of critical & mortal covid-19 cases: A systematic literature review and meta-analysis

Causal associations between risk factors and common diseases inferred from gwas summary data