{'fraction_non_alphanumeric': 0.050215747297520784, 'fraction_numerical': 0.0207027197690677, 'mean_word_length': 4.266349960411718, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
GENIUS-MAWII: For Robust Mendelian Randomization with Many Weak Invalid Instruments

Mendelian randomization (MR) has become a popular approach to study causal effects by using genetic variants as instrumental variables. We propose a new MR method, GENIUS-MAWII, which simultaneously addresses the two salient phenomena that adversely affect MR analyses: many weak instruments and widespread horizontal pleiotropy. Similar to MR GENIUS (Tchetgen Tchetgen et al., 2021), we achieve identification of the treatment effect by leveraging heteroscedasticity of the exposure. We then derive the class of influence functions of the treatment effect, based on which, we construct a continuous updating estimator and establish its consistency and asymptotic normality under a many weak invalid instruments asymptotic regime by developing novel semiparametric theory. We also provide a measure of weak identification, an overidentification test, and a graphical diagnostic tool. We demonstrate in simulations that GENIUS-MAWII has clear advantages in the presence of directional or correlated horizontal pleiotropy compared to other methods. We apply our method to study the effect of body mass index on systolic blood pressure using UK Biobank.


Challenges in Mendelian randomization
Mendelian randomization (MR) is a method of using genetic variants -typically single nucleotide polymorphisms (SNPs) -as instrumental variables (IVs) to infer the causal effect of a modifiable exposure on an outcome in the presence of unmeasured confounding Ebrahim, 2003, 2004;Lawlor et al., 2008;Davey Smith and Hemani, 2014;Burgess et al., 2015a,b;Zheng et al., 2017). As a powerful tool to disentangle causal relationship from complex environmental confounding, MR has become a popular method for establishing high-quality causal evidence based on observational data (Pingault et al., 2018;Adam, 2019).
For reliable causal inference using MR, genetic variants must be valid IVs that satisfy three key assumptions (Angrist and Krueger, 2001;Baiocchi et al., 2014;Hernan and Robins, 2020): (i) (relevance) they are associated with the exposure; (ii) (independence) they are independent of any unmeasured confounder of the exposure-outcome relationship; (iii) (exclusion restriction) they affect the outcome exclusively through the exposure. The first assumption (relevance) is usually satisfied by selecting SNPs that are significantly associated with the exposure. A well-established challenge one is often faced with in MR, is the possibility that individual SNPs are only weakly associated with the exposure, resulting in weak IV bias (Stock et al., 2002; and extreme sensitivity to minor violations of the other two assumptions (Small and Rosenbaum, 2008;Wang et al., 2018). The second assumption (independence) is plausible within the framework of parent-offspring studies because of the random assortments of genes from parents to offspring. The independence assumption also approximately holds in population data such as the UK Biobank as individuals share much common ancestry (Davey Smith et al., 2020). Among the three core IV assumptions, the exclusion restriction assumption is the most disputable, as emerging evidence has suggested that pleiotropy -a phenomenon in which a genetic variant may affect multiple phenotypic traits (Solovieff et al., 2013;Verbanck et al., 2018) -is widespread. In fact, studies have identified hundreds of genetic variants from genome-wide association studies (GWASs) that are associated with multiple traits (Sivakumaran et al., 2011;Parkes et al., 2013;Gratten and Visscher, 2016;Pickrell et al., 2016;Grassmann et al., 2017;Webb et al., 2017). For example, a variant (rs2075650 in the APOE locus) is found to be significantly associated with several traits and diseases, including Body Mass Index (BMI), Alzheimer's disease, C-reactive protein, high-density lipoprotein cholesterol, low-density lipoprotein cholesterol, plasma triglycerides, waist circumference, hip circumference and waist/hip ratio (Verbanck et al., 2018). Hence, using this variant to study the effect of BMI on systolic blood pressure (SBP) will likely violate the exclusion restriction assumption because the variant may affect SBP via other traits outside of the pathway of BMI. Failure to account for such horizontal pleiotropy (i.e., SNPs having direct effects on the outcome) can lead to spurious findings.

Prior work
In this article, we focus on the two salient challenges in MR: many weak IVs and widespread horizontal pleiotropy. These two challenges rarely act alone but rather interact with each other, because weak IVs can amplify bias from pleiotropy (Small and Rosenbaum, 2008), and later in Sections 2-3 we will see that the proposed method accounts for pleiotropy by exploiting heteroscedasticity and thus may be more susceptible to weak IV bias than usual (Lewbel, 2012). Therefore, it is important to address these two challenges simultaneously for desirable practical performance.
The issue of weak IV has been extensively studied in econometrics (Staiger and Stock, 1997;Chao and Swanson, 2005;Hansen et al., 2008;Newey and Windmeijer, 2009;Stock et al., 2002). Typically, in linear models, an IV is considered weak if the first-stage F statistic is below 10 (Stock et al., 2002). Recent papers by Zhao et al. (2019Zhao et al. ( , 2020; Wang and Kang (2021) and Ye et al. (2021) also develop methods that are robust to weak IVs in two-sample summary-data MR under an assumption of no systematic exclusion restriction violation. A common message from many of these works is that having many weak IVs can greatly circumvent the difficulty from each individual IV being only weakly associated with the exposure and can improve estimation accuracy.
There has also been a rapidly growing development of statistical methods to address widespread horizontal pleiotropy, which mostly fall into the following two strands. The first strand of methods assumes that a certain proportion of candidate IVs are valid. For example, Han (2008); Kang et al. (2016); Bowden et al. (2016) and Windmeijer et al. (2019) propose methods that can recover the causal effect provided less than 50% of IVs are invalid. Hartwig et al. (2017); Guo et al. (2018); Guo (2021) and Windmeijer et al. (2021) develop methods based on the plurality rule, assuming that the number of valid IVs is larger than any number of invalid IVs sharing the same ratio estimator limit. Other proposals in this first strand include Qi and Chatterjee (2019) and Verbanck et al. (2018).
Clearly, none of these methods apply to the situation when pleiotropy is pervasive.
The second strand of work allows for all the IVs to be pleiotropic but effectively restricts the effects of IVs on the exposure and outcome. Within the second strand, it is also helpful to distinguish between two types of horizontal pleiotropy: uncorrelated pleiotropy, also known as the instrument strength independent of direct effect (InSIDE) assumption (Bowden et al., 2015), which says that the direct effects of the IVs on the outcome are uncorrelated with their effects on the exposure, and correlated pleiotropy, which says that the direct effects of the IVs on the outcome are correlated with their effects on the exposure.
Of the two types, uncorrelated pleiotropy is easier to deal with, based on which multiple methods have been developed, including Kolesár et al. (2015) and Bowden et al. (2015) for direct effects with nonzero mean (directional horizontal pleiotropy), and Zhao et al. (2019Zhao et al. ( , 2020 and Ye et al. (2021) for direct effects with zero mean (balanced horizontal pleiotropy).
Correlated pleiotropy is more challenging. Morrison et al. (2020) and Wang et al. (2021) allow a small proportion of genetic variants to exhibit correlated pleiotropy arising from one or several well-understood pleiotropic pathways. Tchetgen Tchetgen et al. (2021) and Sun et al. (2022) tackle this challenge from a different perspective; without assuming a certain structure underlying the correlated pleiotropy, the identification extends a novel strategy proposed in Lewbel (2012Lewbel ( , 2018 that exploits heteroscedasticity of the exposure variable. The details are reviewed in Section 2. Other proposals in this second strand include Burgess and Thompson (2015), Spiller et al. (2019) and Liu et al. (2022).

Our contributions
In this work, we propose a new MR method, GENIUS-MAWII, that simultaneously addresses many weak IVs and widespread horizontal pleiotropy. We deal with widespread horizontal pleiotropy by leveraging heteroscedasticity of the exposure, and we account for many weak IVs by establishing the consistency and asymptotic normality of the continuous updating estimator (CUE) obtained from using the derived influence functions as moment conditions under many weak moment asymptotics. We also provide GENIUS-MAWII with a measure of weak identification, an overidentification test, and a graphical diagnostic tool.
We demonstrate in simulations and a real example using UK Biobank the clear advantages of GENIUS-MAWII in the presence of directional or correlated horizontal pleiotropy compared to other methods. Furthermore, our work makes important advances in the theory of generalized method of moments (GMM) involving unknown nuisance parameters under many weak moment conditions, which to our knowledge has not been studied in the literature. This is a challenging task due to two main reasons. First, with the number of moment conditions growing to infinity, the number of nuisance parameters also grows to infinity. Second, the many weak moment asymptotics, which is well suited for MR studies with a large number of SNPs, is fundamentally different from the classical asymptotics (with a fixed number of "strong" moment conditions). Importantly, under the classical asymptotics, it is well known that utilizing the influence function which belongs to the ortho-complement of the nuisance tangent space and estimating the nuisance parameters at a fast enough rate ensure us that the impact of estimating the nuisance parameters is negligible; this is the key insight that drives many other successful applications of using the influence function to handle nuisance parameters (Newey, 1994;Ackerberg et al., 2014;Ning and Liu, 2017;Robins et al., 2017;Chernozhukov et al., 2018;Bravo et al., 2020). To our surprise, we find that this appealing property does not hold under many weak moment asymptotics in general, but still holds for GENIUS-MAWII because its moment conditions are linear in the parameter of interest (see Theorem 2). In addition, our proof handles infinite-dimensional nuisance parameters.
The rest of the article proceeds as follows. In Section 2, we introduce the invalid IV model and review the GENIUS identification strategy. In Section 3, we derive the class of influence functions which are shown to be multiply robust, and the efficient influence function. In Section 4, we consider estimation and inference of the treatment effect. In Section 5, we provide a measure of weak identification, an overidentification test, and a graphical diagnostic tool. The article is concluded with simulations in Section 6, a real data application in Section 7, and more discussion in Section 8. All technical proofs are in the supplementary materials. The R code for the proposed methods is publicly available at https://github.com/tye27/mr.genius.

Review of the GENIUS identification strategy
Suppose that we observe an independent and identically distributed sample (O 1 , . . . , O n ) with O = (Z, X, A, Y ), where Z = (Z 1 , . . . , Z m ) T is a column vector including m SNPs, each taking on values from the set {0, 1, 2} which represents the number of minor alleles, X is a vector of observed covariates which can be empty when there are no observed covariates, A and Y are continuous exposure and outcome variables. We emphasize that the m SNPs Z 1 , . . . , Z m are not required to be independent, i.e., we allow the m SNPs to be in linkage disequilibrium. We are interested in the causal effect of A on Y , denoted by β 0 , in the presence of unmeasured confounders U .
When there are no observed covariates, we consider the following structural equations: where α, γ, ξ y , ξ a are unspecified functions, and Z ⊥ U . In particular, α(Z) encodes the direct effect of Z on Y , and α(Z) = 0 indicates that the exclusion restriction assumption is violated. Lewbel (2012) also considers models (1)-(2). Kolesár et al. (2015) and Bowden et al. (2015) consider the special case with α(Z) = m j=1 α j Z j and γ(Z) = m j=1 γ j Z j , and assume that α j , γ j , j = 1, . . . , m are random effects satisfying α j ⊥ γ j (commonly referred to as the InSIDE assumption or uncorrelated pleiotropy), which is likely violated when there are SNPs affecting the exposure and outcome through common pathways (Morrison et al., 2020). In contrast, we make no such restrictions. Furthermore, as reviewed in Section 1, many existing MR methods, including Kang et al. (2016); Bowden et al. (2016); Hartwig et al. (2017); Guo et al. (2018);Windmeijer et al. (2019) and Guo (2021), rely on the assumption that pleiotropy only sparsely involves a small proportion of SNPs, whereas we allow every SNP to be pleiotropic.
Assume (1)-(2) and Z ⊥ U , it is shown in Tchetgen Tchetgen et al. (2021) that β 0 is the unique solution to We elaborate the key of identification in (3). With Z being potentially invalid IVs that have a direct effect on the outcome, the usual IV-based identification formula no longer In fact, when having a direct effect on the outcome, the invalid IVs Z are nothing more than observed confounders that are independent of U . If the effect of Z on the outcome is not modified by U , as is the case under (1), then (Z − E(Z))c(U ) for any function c(·) satisfying E(c(U )) = 0 can be conceptualized as "valid IVs" satisfying E{ are uncorrelated with any function of U and any function of Z, and do not have a direct effect on the outcome. The conceptualized valid IVs (Z − E(Z))c(U ) are infeasible as U is unobserved, but under (2) a noisy version of which can be constructed as an additive interaction between conditionally centered exposure and centered IVs (Z − E(Z))R A and is used as the feasible "valid IVs".
There are three comments about the above intuition. First, the idea of using geneenvironment interactions as valid IVs also appears in Spiller et al. (2019) (2), its magnitude is usually small compared to the main effects of Z and U , then the bias of GENIUS is also relatively small. Finally, the key condition encoded by (1)-(2), i.e., the effects of Z on the exposure and outcome not being modified by U is stronger than needed and can be relaxed to some extent (see Section 3.1 of the supplementary materials).
We can also circumvent this restriction by collecting information about the part of U that interacts with SNPs and adjust for them as part of the observed covariates. This will be discussed further in Section 3.
When m = 1 (i.e., one SNP), β 0 identified via (3) can be rewritten as a Wald ratio where the numerator is the effect of (Z − E(Z))R A on Y , the denominator is the effect of (Z − E(Z))R A on A, and β 0 is simply the ratio. When m > 1 (i.e., multiple SNPs), β 0 is over-identified. Moreover, identification using (3) requires that E{(Z − E(Z))R A A} = 0, which is analogous to the relevance assumption in the IV literature, except here we con- , which means that identification using (3) requires A being heteroscedastic, i.e., Var(A | Z) depends on at least some Z. We remark that the condition E{(Z − E(Z))R A A} = 0 can be empirically checked since Cov(Z, R 2 A ) can be estimated by the sample covariance between Z and the squared residuals from fitting a linear regression of A on Z. One can also apply tests for heteroscedasticity such as the tests in Koenker (1981) and White (1980). Heteroscedasticity can be due to gene-environment interactions (Paré et al., 2010); see Wang et al. (2019) and Sulc et al. (2020) for some recent discoveries.

Semiparametric theory
The identification result in Section 2 can be easily extended when there is an observed covariate vector X. Consider the following structural equations: where α, γ, ξ y , ξ a are unspecified functions and Z ⊥ U | X. Then, as shown in Section 3.2 of the supplementary materials, β 0 is the unique solution to is the conditionally centered exposure. Hence, identification by (6) requires that Var(A | Z, X) depends on some Z. Note that the GENIUS identification strategy can be extended to binary exposure and/or binary outcome that follow semiparametric log-linear models; see Section 6 of the supplementary materials for details.
Comparing structural equations (4)-(5) and Z ⊥ U | X with their unconditional counterparts, we see that to satisfy these assumptions, X should include covariates that (i) are correlated with Z; (ii) modify the effect of Z on the outcome (which should be rare as Z should primarily influence the exposure); (iii) are confounders of the exposure-outcome relationship and modify the effect of Z on the exposure. Another interesting type of covariates is those that do not affect the outcome but modify the effect of Z on the exposure.
Adjusting for these covariates can weaken heteroscedasticity and thus weaken identification, and may even make the exposure effect unidentifiable if conditioning on all such covariates.
However, identification of exposure effect can still be achieved if there is residual latent heterogeneity in the effect of Z on A within all levels of X. A diagram of how to choose X is in Section 1.1 of the supplementary materials.
We derive the class of influence functions and the efficient influence function (Bickel et al., 1993) under the sole observed data restriction , let h(Z, X) be any scalar-valued function, the class of influence where andZ is a column vector of all the dummy variables for the joint levels defined by Z.
The proof is given in the supplementary materials. Theorem 1 includes the results without observed covariates as a special case by setting X to be empty.
Identification using the influence function (7) is in fact multiply robust. As shown in the supplementary materials, the influence function (7) evaluated at β = β 0 has expectation zero when either one of the following three sets of the models is correctly specified: Therefore, in classical settings, multiply robust estimation and inference about β 0 is straightforward via the classical GMM results (Hansen, 1982).

Estimation and inference with many invalid IVs
Since we have m SNPs, it is natural to use all m SNPs to improve estimation accuracy.
Without loss of generality, we consider a m-dimensional vector of moment conditions where ∆ is defined in Theorem 1, η denotes the vector of nuisance parameters, and η 0 is its true value. With β 0 being over-identified by (8), we can solve for β 0 using GMM.
We introduce some additional notations. Let Note that G i (η) does not depend on β because g i (β, η) is linear in β. This largely simplifies the problem.
As always, asymptotic theory is useful if it provides a good approximation to finitesample performance in applications. In MR with a large number of SNPs while each individual SNP is only weakly related to the exposure, many weak moment asymptotics is well-suited and provides an improved approximation to finite sample behavior of invalid IV robust inference than the classical asymptotics with a fixed number of "strong" moment conditions as n goes to infinity (see simulations in Section 6). Now we are ready to give the formal characterization of the many weak moment asymptotics.
Assumption 1 provides an improved approximation when the many moment conditions are weak. When µ n = √ n and m is finite, it agrees with the classical asymptotics (a finite number of "strong" moment conditions). More discussion on µ n is in Section 1.3 of the supplementary materials.
The many weak moment asymptotics is fundamentally different from the classical asymptotics. Analogous to the weak IV bias arising from linear models, it has also been recognized that many weak moment conditions can make the usual GMM inference inaccurate (Stock et al., 2002). For example, Newey and Windmeijer (2009) find that the two-step GMM is biased and has non-normal asymptotic distribution, while estimators in the generalized empirical likelihood (GEL) family (Smith, 1997;Parente and Smith, 2014) are consistent and asymptotically normal but have larger asymptotic variance than usual.
Furthermore, as outlined in Section 1.3, the many weak moment asymptotics poses several technical difficulties on dealing with nuisance parameters, which to our knowledge has not been addressed in the literature. For the rest of this section, we develop novel semiparametric theory to handle unknown nuisance functions under many weak moment conditions.
These theoretical developments enable fast and stable estimation and inference about β 0 .
For estimation purposes, we will assume linear SNP (Zhang and Sun, 2021), exposure, and outcome models in Assumption 2.
Assumption 2. (nuisance parameters) Suppose that Z, X are bounded, X ∈ R dx with where ω 0 and θ 0 are unspecified functions, and the first component of X is 1 representing the intercept term.
Under Assumption 2, η = (π T 1 , . . . , π T m , µ T , λ T , ω(x; µ, λ), θ(x; µ)) T collects all the nuisance parameters, η 0 is the true value of η. Write the estimator of η 0 asη, which includes the least squares estimatorsπ 1 , . . . ,π m ,μ,λ from fitting the linear models in Assumption 2, and the kernel estimators with plug-in estimated parametersω(x;μ,λ),θ(x;μ). We choose the more flexible kernel estimators for ω 0 (x; µ, λ) and θ 0 (x; µ) to avoid modeling the second moment terms. Alternatively, one can assume that E(R A R Y | X) and E(R 2 A | X) follow parametric models, for example, linear models that include a full set of quadratic terms of X or saturated models when X consists of only categorical variables. Then all the nuisance parameter estimatorsη can be obtained from the least squares estimation. Either way, the estimated nuisance parametersη converge to their true values under the assumed conditions. A special case is when not adjusting for covariates, i.e., X only includes the intercept term, ω(x; µ, λ) and θ(x; µ) become two one-dimensional parameters, and their estimators degenerate to simple averages.
We focus on the continuous updating estimator (CUE) -a member of the GEL familyin this article, because its objective function has an explicit form and its empirical performance is similar to the other estimators in the GEL family (Newey and Smith, 2004). We propose the following GENIUS estimator that leverages MAny Weak Invalid IVs (GENIUS-MAWII), which is obtained using the influence function g IF (O; β, η) defined in (8) with a plug-in nuisance parameter estimatorη (defined above), where B is a compact set of parameter values. Note that the form of the CUE is similar to the familiar two-step GMM estimator, except that the objective function is simultaneously minimized over β in the optimal weighting matrixΩ(β,η). This is key in eliminating the many weak moment bias of two-step GMM estimator (Newey and Windmeijer, 2009).
The following theorem establishes the asymptotic properties ofβ defined in (9).
Taylor expansion of the first-order condition ∂Q(β,η)/∂β| β=β = 0 gives whereβ is some value between β 0 andβ. We prove that nµ −1 n ∂Q(β,η)/∂β| β=β 0 is asymptotically equivalent to the sum of a usual GMM term nµ −1 n G T Ω −1ĝ (β 0 , η 0 ) and a U-statistic term (nµ n ) −1 i =j U T i Ω −1 g j that is no longer negligible due to the many weak moment asymptotics, and both terms are mean zero. The asymptotic normality of nµ −1 n ∂Q(β,η)/∂β| β=β 0 then follows from the U-statistic term being uncorrelated with the usual GMM term, the asymptotic variance of the usual GMM term being nµ −2 n G T Ω −1 G, the asymptotic variance of the U-statistic term being µ −2 n E[U T i Ω −1 U i ], and the central limit theorem. The result in (10) follows from showing nµ −2 n ∂ 2Q (β,η)/∂β 2 | β=β = nµ −2 n G T Ω −1 G + o p (1). Interestingly, the outline above implies thatβ defined in (9) is asymptotically equivalent to arg min β∈BQ (β, η 0 ), which means that estimation of η 0 does not affect the asymptotic distribution of the CUE, even with multiple complications arising from the many weak moment asymptotics, the number of nuisance parameters growing to infinity, and nonparametric kernel estimators that themselves involve estimated parameters. This result relies on the estimated nuisance parameters converging to the truth. But still, this is an unusual property that does not hold in general, but holds in the current setting due to three factors: (i) the number of moment conditions m grows to infinity at a rate slower than n 1/3 ; (ii) the use of influence functions as the moment conditions and the fast convergence rate of the estimated nuisance parameters which imply that √ n ĝ(β 0 ,η)−ĝ(β 0 , η 0 ) = o p (1); and (iii) the where · is the 2 vector norm. Crucially, we need (iii) to make sure that the impact of estimating η 0 is negligible for the U-statistic term; in contrast, (iii) is not needed under the classical asymptotics because the U-statistic term is a higher order term that is negligible.
In Theorem 2, the number of SNPs m is required to grow slower than the sample size n, which is more restrictive than the limited information maximum likelihood (LIML) estimator (Chao and Swanson, 2005) where m can grow at the same rate as n or even faster. The reason behind this difference is that LIML assumes homoscedasticity, while CUE makes no such assumption. Consequently, for consistency of the CUE, m 2 /n → 0 seems necessary given the need to consistently estimate the heteroscedastic weight matrix Ω which has m 2 elements.
According to Theorem 2, under the many weak moment asymptotics, as long as m 3 /n is small, the CUEβ is consistent and asymptotically normal, and the convergence rate is µ n . The asymptotic variance ofβ consists of the limit of two terms: The first term is the classical GMM variance, while the second term is the variance contribution due to the variability of the moment derivative G i which does not vanish under many weak moment asymptotics. Specifically, when m/µ 2 n → 0 or G i is a constant, the second term is negligible compared to the first term, so that (10) agrees with classical GMM theory; otherwise, the additional variance is not negligible and results in larger variance of the CUE. Interestingly, in this many weak moment asymptotic regime where identification is weak, the impact of estimation of the weight matrix Ω is small compared to estimation of G and does not appear in the variance formula.
In practice, we can estimate the asymptotic variance ofβ usingV /n, wherê Here,Ĥ is an estimator of G T Ω −1 G, the middle termD TΩ−1D is an estimator of the asymptotic variance of √ n∂Q(β,η)/∂β| β=β 0 . Notice thatĜ TΩ−1Ĝ cannot be used in place ofĤ becauseĜ TΩ−1Ĝ is biased under many weak moment asymptotics (Newey and Windmeijer, 2009). Based on the variance estimator, we can test the hypothesis H 0 : β = β * using the Wald statistic T = √ n(β − β * )/V 1/2 or we can construct a 1 − α confidence interval based on normal approximation. Other identification robust statistics, e.g., Lagrange multiplier statistic, conditional likelihood ratio statistics can also be applied here, and they are in fact asymptotically equivalent to the Wald statistic under many weak moment asymptotics; see Newey and Windmeijer (2009) for details. It is worth noting that under strong identification, Theorem 2 and the variance formula in (12) are asymptotically equivalent to classical GMM counterparts, and thus remain applicable. Hence, the results in this section are not only well-suited for MR analysis but also applicable to a wider range of regimes compared to classical GMM results.
The exposure and outcome models in Assumption 2 can be extended to include SNP-SNP and SNP-covariate interactions, by simply fitting linear models with those interaction terms included. Then, consistency in Theorem 2 continues to apply when the numbers of regressors in the exposure and outcome models go to infinity slower than n 1/2 , and the asymptotic normality result in Theorem 2 continues to apply when the numbers of regressors in the exposure and outcome models go to infinity slower than n 1/3 .
Finally, as discussed after Theorem 1, estimators based on influence functions enjoy a multiple robustness property in classical settings with a finite number of "strong" moment conditions, andβ achieves the semiparametric efficiency bound (Ackerberg et al., 2014).
However, multiple robustness and semiparametric efficiency under many weak moment asymptotics are more delicate and will be interesting to pursue in future work.

Measure of weak identification and diagnosis
We have developed a new MR method, GENIUS-MAWII, which simultaneously addresses the two salient challenges in MR: many weak IVs and widespread horizontal pleiotropy.
We account for many weak IVs by establishing consistency and asymptotic normality of the CUE obtained from using the derived influence functions as moment conditions under many weak moment asymptotics, and we deal with widespread horizontal pleiotropy by leveraging heteroscedasticity of the exposure based on the GENIUS identification strategy.
However, GENIUS-MAWII will break down if identification is too weak for many weak moment asymptotics to provide good approximation, resulting in estimation bias; or if the untestable assumptions (4)-(5) and Z ⊥ U | X do not hold, resulting in identification bias. Therefore, to enhance the reliability of GENIUS-MAWII in MR analysis, besides using domain knowledge, it is useful to have tools to gauge whether identification is strong enough for the promised asymptotic results to kick in and whether there is any evidence to falsify the assumptions. For these two purposes, we present a measure of weak identification, an overidentification test, and a graphical diagnostic tool in this section.

Measure of weak identification
Detection of weak identification is an important part in IV analysis, because weak identification can result in unreliable estimation and inference (Stock et al., 2002). Until now, several formal procedures are readily available for weak IV detection based on linear IV models.
For example, Stock and Yogo (2001) propose to use the first-stage F-statistic to assert whether IVs are weak; Hahn and Hausman (2002) develop a specification test for strong IVs. In this section, we provide a measure of weak identification for GENIUS-MAWII, which measures the extent of heteroscedasticity and can serve as a helpful diagnosis for reliable inference.
Write the moment equations in (8) as which can be viewed as using as the derived outcome and R 2 A − E(R 2 A | X) as the exposure. This motivates using the heteroscedasticity-robust F-statistic in the regression of R 2 as a measure of weak identification. This F-statistic can also be thought of as a Koenker test for heteroscedasticity (Koenker, 1981), specifically for whether Var(R A | Z, X) depends on Z. In practice, we replace the unknown quantities by their estimators and denote the heteroscedasticity-robust F-statistic in the regression ofR 2 represents an estimator of the conditional expectation E(B | C). As according to Stock and Yogo (2001 , Table 4) and Stock et al. (2002, Chap. 6.2), for a large m, the F-statistic needs to be larger than 2 for LIML inference to be reliable. Because CUE is a GMM-like generalization of LIML, we recommend check to make sure F GENIUS is larger than 2. This strategy is empirically evaluated in Section 6.2.

Overidentification test and graphical diagnosis
In the GMM literature, it is common to perform overidentification tests to test whether (Hansen, 1982). A popular statistic is simply a scaled minimized CUE objective function 2nQ(β,η), which is often called the J-statistic.
Theorem 3. Under the same conditions in Theorem 2 and m 3 /n → 0, when the null Theorem 3 shows that we can reject H 0 if 2nQ(β,η) ≥ χ 2 1−α (m − 1), which is the same as the overidentification test under the classical setting with a fixed number of "strong" moment equations. In Section 1.2 of the supplementary materials, we analytically show that the overidentification test has power to detect assumption violations in typical MR applications.
In addition to the overidentification test, another diagnosis approach is as follows. Note that our method relies on untestable assumptions (4)-(5) and Z ⊥ U | X. These assumptions imply the conditional moment restriction We have used part of its implications to identify β 0 in (8). But this conditional moment restriction has many other implications that we can use for falsification. This motivates a graphical diagnostic tool which plots the "residual" is a pre-specified function that is non-linear in z to avoid using duplicated information as that used for identification in (8). If modeling assumptions hold and all estimators have negligible bias, thent i should be centered around zero across different values of f (Z i , X i ); evidence thatt i is not centered around zero indicates violation of the assumptions. Note that unlike the typical residual plots for diagnostic in regression models, variance difference oft i across different values of f (Z i , X i ) does not violate our assumptions. As a final remark, the graphical diagnostic tool is also applicable to the situation where there is no observed covariates by setting X to be empty. In Section 7, we set f

A simulation when assumptions for GENIUS-MAWII hold
We conduct a simulation study to evaluate the finite-sample performance of GENIUS-MAWII when its assumptions hold, i.e., under models (4)-(5) and Z ⊥ U | X. Its performance is compared to the GENIUS estimators obtained from using ordinary CUE and two-step GMM. All three GENIUS estimators are based on the influence functions in (8), with nuisance parameters estimated using the models in Assumption 2. GENIUS-GMM is the two-step GMM from the gmm package in R. GENIUS-MAWII and GENIUS-CUE have the same point estimators defined by (9) and are computed in the same way using the optimize and uniroot functions in R. Specifically, we first minimize the objective function in (9) using optimize function with the specified boundary (-10, 10), and if the returned value is very close to the specified boundary, we recompute by applying uniroot to the derivative of the objective function in (9). The main purpose of this extra step is to stabilize the numerical optimization. The difference between GENIUS-MAWII and GENIUS-CUE is in the variance estimators: GENIUS-CUE uses the classical textbook variance estimator, while GENIUS-MAWII uses (12).
GENIUS-MAWII is also compared to two IV estimators and five other MR estima-tors in the literature: the two-stage least squares (2SLS), limited information maximum likelihood (LIML), the inverse variance-weighted (IVW) estimator (Burgess et al., 2013), robust adjusted profile score estimator (MR-raps) (Zhao et al., 2020), MR-Egger regression (Bowden et al., 2015), MR-median (Bowden et al., 2016), and MR-mode (Hartwig et al., 2017). These five MR methods are developed as two-sample MR methods, but we apply them to our one-sample setting regardless to see their performance. 2SLS is implemented using the AER package, LIML using the ivmodel package, IVW using the mr.divw package, MR-raps using the mr.raps package, and MR-Egger, MR-median, and MR-mode using the MendelianRandomization package.
We generate m independent SNPs Z 1 , . . . , Z m with P (Z j = 0) = 0.25, P (Z j = 1) = 0.5, P (Z j = 2) = 0.25, j = 1, . . . , m. We generate the exposure and outcome from where m is the total number of SNPs, ϕ γj , ϕ δj 's are constants that are generated once from a standard normal distribution. Here, h 2 can be interpreted as the proportion of variance in A that is attributed to E(A | Z), and κ controls the level of heteroscedasticity.
As illustrated in Figure 1, we consider three types of SNPs: S 1 , S 2 , S 3 with proportion p 1 , p 2 , p 3 , where p 1 + p 2 + p 3 = 1. Specifically, S 1 consists of valid IVs with α j = 0 (i.e., no direct effect on the outcome); S 2 consists of invalid IVs with uncorrelated pleiotropic effects and α j ∼ i.i.d. N ( τ 2 0 , τ 2 0 ) (i.e., InSIDE is satisfied), where τ 2 0 = h 2 /(1.5m); S 3 consists of invalid IVs that affect A and Y through a common factor C, which leads to correlated pleiotropic effects and α j = γ j /2. We consider four settings: Figure 1: Illustration of three types of SNPs: S 1 consists of valid IVs, S 2 consists of invalid IVs with uncorrelated pleiotropic effects (i.e., InSIDE is satisfied), and S 3 consists of invalid IVs with correlated pleiotropic effects.

(40% invalid IVs
We set h 2 = 0.2 and κ = 1. We consider two values of the sample size: n = 10, 000 and n = 100, 000. The results with 1,000 Monte Carlo repetitions are in Table 1, which summarizes (i) the Monte Carlo mean and Monte Carlo standard deviation (SD) of each estimator, (ii) average of standard errors (SEs), and (iii) coverage probability (CP) of 95% confidence intervals from normal approximation.
From Table 1, the performance of GENIUS-MAWII is similar across Setting 1-4. When n = 10, 000, the average F-statistic F GENIUS = 1.83, GENIUS-MAWII shows nominal coverage probability, but has some attenuation bias and its SD is slightly more than √ 10 = 3.16 times larger than the SD with n = 100, 000. This is because when the sample size is small, there are some outliers in some simulation runs due to instability of numerical optimization. Numerical optimization becomes more stable as the sample size becomes larger, essentially resolving this issue. When n = 100, 000, the average F-statistic F GENIUS = 9.48 and GENIUS-MAWII shows negligible bias and nominal coverage, which agrees with our theoretical assessment that GENIUS-MAWII performs well when the identification is not too weak. Across all scenarios, the SEs calculated using (12) are close to Monte Carlo SDs of GENIUS-MAWII. Notice that with m = 100 and n = 10, 000, m 2 /n = 1 and m 3 /n = 100 are not small, but the GENIUS-MAWII estimator still performs quite well, indicating that our method is able to work well in typical MR studies with around 100 SNPs and 10,000-500,000 sample size.
Under all scenarios, the SEs underestimate the Monte Carlo SDs of the GENIUS-CUE estimator, which is also reflected by the fact that the CPs are below the nominal level 95%.
This is expected because according to our Theorem 2, a higher order variance term is no longer negligible under many weak moment asymptotics. Comparing the GENIUS-CUE estimator and the GENIUS-MAWII estimator, we see that many weak moment asymptotics indeed provides a better finite sample approximation.
Across all simulation scenarios, the GENIUS-GMM estimator has a larger bias than GENIUS-MAWII and GENIUS-CUE estimators, especially when n = 10, 000. Moreover, it is not difficult to derive that the ordinary least squares (OLS) estimates obtained from regressing Y on intercept, A, and Z 1 , . . . , Z m is approximately β 0 + 1/(1 + h 2 ) = 1.23.
Hence, we see that the GENIUS-GMM estimator is in fact biased towards the OLS. This is a GMM version of the well-known phenomenon that 2SLS is biased towards the OLS The MR-Egger estimator is valid under Setting 3. However, in Setting 3, MR-Egger is biased even when n = 100, 000, and its SD and SE are larger compared to GENIUS-MAWII.
Under Setting 2 and 4, the MR-Egger estimator is also biased.
In Section 2 of the supplementary materials, we conduct similar simulation studies where there is an observed covariate or when the SNPs are dependent (i.e., SNPs are in linkage disequilibrium). Similar to the results presented in Table 1, we generally find that the GENIUS-MAWII estimator has desirable performance with negligible bias and nominal coverage probability.

A simulation under assumption violation and weak identification
We conduct more simulations for GENIUS-MAWII under assumption violation and weak identification, and demonstrate the use of F-statistic and diagnostics tools to help identify situations where GENIUS-MAWII can be reliably applied. The setting is identical to that in Section 6 except that we generate the outcome from (4) is violated, and we generate the exposure from and ϕ Y ηj , ϕ Aηj 's are constants that are generated once from a standard normal distribution. In other words, 20% SNPs can have interactions with the unmeasured confounder, and the magnitude of which is similar to their interactions with A .
We consider four situations: no model assumption is violated, only the exposure or the outcome model assumption is violated, and both model assumptions are violated. In each situation, we consider κ = 0, 0.1, 0.5, 1 for increasing level of heteroscedasticity, and n = 10, 000, 50, 000 and 100, 000. Note κ = 0 means there is no model violation. The results are in Table 2.
From Table 2, when F GENIUS > 2, GENIUS-MAWII has negligible bias and nominal coverage probability when there is no model violation. When there is model misspecification and F GENIUS > 2, GENIUS-MAWII has nontrivial power (more than 50% power) to detect model misspecification. Therefore, although assumption violations and/or weak heteroscedasticity can severely bias the GENIUS-MAWII estimator, the combined use of F GENIUS and overidentification test can effectively identify those situations and provide guidance about when GENIUS-MAWII can be reliably applied.
7 Application to the UK Biobank data UK Biobank is a large-scale ongoing prospective cohort study with around 500,000 participants aged 40-69 at recruitment from 2006 to 2010. Participants provided biological samples, completed questionnaires, underwent assessments, and had nurse led interviews.
Follow up is chiefly through cohort-wide linkages to National Health Service data, including electronic, coded death certificate, hospital, and primary care data (Sudlow et al., 2015).
Prevalent disease was coded using ICD-9 and ICD-10, and cause of death was coded using ICD-10. Genotyping was performed using two arrays, the Affymetrix UK BiLEVE (UK Biobank Lung Exome Variant Evaluation) Axiom array (about 50,000 participants) and Affymetrix UK Biobank Axiom array (about 450,000 participants). The SNPs included for analysis were directly genotyped or imputed using the Haplotype Reference Consortium panel. To reduce confounding bias due to population stratification, we restrict our analysis to people of genetically verified white British descent, as in previous studies (Tyrrell et al., 2016). For quality control, we exclude participants with (1) excess relatedness (more than 10 putative third-degree relatives), or (2) mismatched information on sex between genotyp-ing and self-report, or (3) sex-chromosomes not XX or XY, or (4) poor-quality genotyping based on heterozygosity and missing rates > 2%.
We are interested in estimating the causal effect of body mass index (BMI) on systolic blood pressure (SBP). We also exclude participants who are taking blood pressure medication based on self report. In total, the sample size for the final analysis is 292,757. We use 93 SNPs that are associated with BMI at genome-wide significance level (Locke et al., 2015).
We apply our method to the UK Biobank data and the results are summarized in Table 3. The implementation details are the same as in Section 6. From Table 3, the Fstatistic for the standard IV is F IV = 59.7. The 2SLS, LIML, IVW, and MR-raps all have point estimates larger than that from GENIUS-MAWII and MR-Egger, which is likely due to failing to account for horizontal pleiotropic effects with nonzero mean. Moreover, the Sargan test (Sargan, 1958)  Finally, we run the tools developed in Section 5 to assess the strength of identification and plausibility of the assumptions. The heteroscedasticity robust F-statistic for GENIUS is F GENIUS = 11.1, large enough for application of GENIUS-MAWII. The overidentification test statistic is 2nQ(β,η) = 108.3, smaller than the critical value χ 2 0.95 (92) = 115.4. In addition, the diagnostic plot in Figure 2 shows that the blue line, which is the estimated conditional mean using the smoothing splines, is close to a straight horizontal line through zero. Therefore, both diagnosis approaches find no evidence of assumption violation in this application.

Discussion
In this paper, we have developed GENIUS-MAWII, a new method for Mendelian randomization (MR) which simultaneously addresses the two salient phenomena that adversely affect MR analyses: many weak IVs and widespread horizontal pleiotropy.
We show via theory and simulations that GENIUS-MAWII can incorporate a large number of SNPs, allows for every SNP to be pleiotropic, and is able to account for directional or correlated horizontal pleiotropy. These features make GENIUS-MAWII stand out with clear advantages over existing methods in the presence of directional or correlated horizontal pleiotropy. In an application to the UK biobank data to study the effect of BMI on SBP, GENIUS-MAWII produces a plausible effect size estimate (β = 0.140, 95% CI: [0.005, 0.275]), whereas 2SLS, LIML, IVW, and MR-raps produce larger effect size estimates (β ranges from 0.277 to 0.482) that are likely due to failing to account for horizontal pleiotropic effects with nonzero mean. In addition, MR-Egger produces an effect size estimate (β = 0.175) that is of similar magnitude to that from GENIUS-MAWII, but is much less precise and fails to yield statistical significance.
GENIUS-MAWII leverages heteroscedasticity to identify the causal effect, which can occur due to gene-environment interactions and is plausible for many situations (Paré et al., 2010;Wang et al., 2019;Sulc et al., 2020). However, if the degree of heteroscedasticity is not very strong or the sample size is not large, estimation and inference may become challenging due to weak identification, and certain deviations away from the assumptions can generate large biases. Therefore, we recommend to perform the overidentification test and the graphical diagnosis to check for any evidence of assumption violation, and check to make sure the GENIUS F-statistic is larger than 2. These tests are very useful in determining whether GENIUS-MAWII can be applied reliably.
Finally, we have developed novel semiparametric theory for handling unknown nuisance parameters under many weak moment conditions, which to our knowledge has not been studied in the literature. Our theory addresses three main technical challenges: (i) the number of weak moment conditions grows to infinity with the sample size; (ii) the number of nuisance parameters grows to infinity with the sample size; and (iii) there exist infinitedimensional nuisance parameters. Our theoretical developments enable fast and stable estimation and inference about the causal effect of interest.

Supplementary Materials
The supplementary materials contain all technical proofs, identification results for binary exposure and/or binary outcome with the exponential link, and additional analytical and simulation results. Table 1: Simulation results based on 1,000 Monte Carlo repetitions with β 0 = 0.4 and m = 100; F GENIUS represents the average F-statistics for GENIUS-MAWII discussed in Section 5.1, SD is the Monte Carlo standard deviation, SE is the average standard error, CP is the coverage probability of 95% asymptotic confidence interval. For each setting, the methods of which the assumptions hold are in boldface.

Identification under assumption violation
Consider the following exposure and outcome models: For any SNP Z j , the identified parameter is When there exist Z-U interactions in the specification of Y and A, i.e., when ξ a (Z) and/or ξ a (Z) are not constants, β j would typically depend on j in MR applications, so there exists no β that satisfies all the moment conditions, under which the overidentification test has power to detect assumption violations.
One situation when the overidentification test has no power to detect assumption violations is when the β j is equally biased for all j. This can happen under a peculiar situation when ξ a (Z) = ξ a0 + ξ az m j=1 Z j , ξ y (Z) = ξ y0 + ξ yz m j=1 Z j , and σ(Z) = σ 0 + σ z m j=1 Z j , and Z j 's are mutually independent and identically distributed, i.e., when the role of each Z j is homogeneous.

Interpretation of µ n
The key quantity that determines the asymptotic variance and convergence rate ofβ − β 0 in (11) is nG T Ω −1 G. Moreover, from Assumption 1, the many weak moment asymptotics requires µ n → ∞, which corresponds to nG T Ω −1 G going to infinity. We give a simple example to provide an interpretation of nG T Ω −1 G. Suppose that . . , Z m are mutually independent. We again emphasize that the mutual independence of Z 1 , . . . , Z m is not needed for our method but is assumed here for ease of illustration. The influence function in (8) suppressing X evaluated at β = β 0 is By definition, expressions for G, Ω and nG T Ω −1 G can be written as where the approximation for Ω in (S2) holds when σ(Z) − E{σ(Z)} is small in magnitude compared to E{σ(Z)} (shown below). Hence, for identification to be strong in the sense that nG T Ω −1 G is of order n, we need m j=1 {Cov(Z j , Var(A | Z))} 2 /Var(Z j ) to be of a constant order. A closer look at the formula indicates that m j=1 {Cov(Z j , Var(A | Z))} 2 /Var(Z j ) = m j=1 a 2 j Var(Z j ), which measures the total variance of Var(A | Z) that can be explained by the set of m SNPs, where a j = Cov(Z j , Var(A | Z))/Var(Z j ) is the population coefficient from the least squares regression of Var(A | Z) on Z j . Evidently, the use of many SNPs that are predictive of Var(A | Z) can strengthen identification and improve estimation accuracy. Now we show the approximation for Ω in (S2) when σ(Z) − E{σ(Z)} is small in magnitude compared to E{σ(Z)}. Letξ a (U ) = ξ a (U ) − E{ξ a (U )} andξ y (U ) = ξ y (U ) − E{ξ y (U )}, and we have It is not difficult to show that where the fifth equality is from A ⊥ (Z, U, Y ), the sixth equality uses the condition that σ(Z) − E{σ(Z)} is small compared to E{σ(Z)}.

When there are observed covariates
We evaluate the finite sample performance of the proposed GENIUS-MAWII estimator when there exists an observed covariate and we estimate E(R A R Y | X = x) and E(R 2 A | X = x) using the nonparametric kernel regression or least squares estimation (LSE) with quadratic terms. The setting is identical to that in Section 6.1 except that we generate the outcome from Y = β 0 A+ m j=1 α j Z j +η Y U (1+X)+ Y , the exposure from A = m j=1 γ j Z j + η A U (1 + X) + 1 + m j=1 δ j Z j A , where X ∼ N (0, 0.6(1 − h 2 )). The kernel regression is implemented using R package np with default options, where the bandwidth is selected via cross-validation. Table S1 presents the results of GENIUS-MAWII (kernel), GENIUS-MAWII (LSE), 2SLS, and LIML under Setting 4; the other two GENIUS estimators are not included because we have shown in Table 1 that GENIUS-MAWII has better performance, the other four MR methods are not included because they cannot adjust for observed covariates. The conclusions from Table S1 are similar to those in Table 1. Specifically, when n = 100, 000, the two GENIUS-MAWII estimators are similar, showing negligible bias and nominal coverage, and the SEs calculated using (12) are close to the simulation SDs; when n = 10, 000, F GENIUS is too small and the two GENIUS-MAWII estimators have some attenuation bias. In contrast, the 2SLS and LIML are biased because of failing to address the horizontal pleiotropy.   Table S2, based on which the conclusions are similar to that in Table 1.   (3) We present a weaker assumption under which β 0 is the unique solution to (3). Here, α, γ, ξ y , ξ a are unspecified functions satisfying α(U, 0) = γ(U, 0) = 0. (6) Under (4)-(5) and Z ⊥ U | X, then

Proof of
where the last equality is from Z ⊥ U | X. Therefore, To show β 0 is the unique solution to (6), note that

Proof of Theorem 1
Let P θ denote a parametric submodel with P 0 = P , where P is the true distribution of the observed data O = (Y, A, Z, X). The score corresponding to P θ with density dP θ is for any scalar-valued function h(Z, X). Following Newey (1994), we obtain the influence functions for estimation of β 0 by deriving the pathwise derivatives ∂β(θ)/∂θ| θ=0 based on (S4). For each h, differentiating under the integral yields ]. Then, we use the following identities repeatedly: Consider the terms A 1 , A 2 , A 3 separately: Since this derivation holds for each h, the set of influence functions of β 0 is given by {g IF (O; h)}. In addition, we note that L h = E[var(A | Z, X)R h Z ]; a necessary condition for non-singularity of L h is the dependence of var(A | Z, X) on Z, an assumption made in the main manuscript.
(b) First, notice that because every Z j is discrete and takes values 0,1,2, any function h(Z, X) can be expressed as C(X) TZ , whereZ contains all the dummy variables of the strata defined by Z. Define d(Z, X) =Z − E(Z | X), all the influence functions of β 0 can be characterized by From the property of the efficient influence function (Tsiatis, 2007), for any C(X), where U (C opt ) corresponds to the efficient influence function. Thus, Since the above equation holds for any C(X), it holds when which implies that (S8) equals zero almost surely. Therefore,

Proof of multiple robustness
We will show that the influence function in (7) evaluated at β = β 0 has expectation zero when either one of the following three sets of the models is correctly specified: In the following, we use a bar notation to denote a general specification of the model that is not necessarily correct.
In the first scenario, In the second scenario,Ē(Y |Z, In the third scenario, where the last line is from the derivations in the first scenario.

Proof of results in Section 4
For any vector x, we denote its 1 , 2 and ∞ norms by Let η = (η T p , η T np ) T ∈ R q collect all the nuisance parameters, with η p = (π T 1 , . . . , π T m , µ T , λ T ) T collecting the finite-dimensional parameters and η np = (ω(x), θ(x)) T collecting the infinitedimensional functional parameters. Denote the true value of η as As mentioned in the main article, if one would like to include interactions terms in the exposure and outcome models, we can simply include those interaction terms in V and then the following proof still goes through as long as dim(V ) = O(m).
For two sequences of real numbers a n and b n , we write a n = O(b n ) if |a n | ≤ Cb n for all n and some C > 0, a n = o(b n ) if a n /b n → 0 as n → ∞, a n = Θ(b n ) if Cb n ≤ |a n | ≤ C b n for all n and some C, C > 0. We use p − → to denote convergence in probability, d − → to denote convergence in distribution. For random variables X and Y , we denote is bounded in probability. We will use w.p.a.1. as abbreviation for with probability approaching 1.
Interestingly, this property also holds for G i (η) because g i (β, η) is linear in β, which can be seen from We will show later that this property is crucial for the estimation of η to have negligible impact on the distribution ofβ.

Regularity conditions
Assumption 3 (Kernel). (i) K(u) is bounded, K(u) is zero outside a bounded set, K(u)du = 1 and uK(u)du = 0; (ii) E(Y 8 | X)f 0 (X) and E(A 8 | X)f 0 (X) are bounded, and the density f 0 (x) is bounded away from zero in the support of X; (iii) The bandwidth of kernel estimator σ satisfies σ 4 √ nm → 0 and σ dx √ n/( √ m log n) → ∞ as n → ∞.

Lemmas
We will first prove some lemmas, which will be used in the proof of Theorem 2 in Section 4.3. The organization of the proof is illustrated in Figure 4.
The first lemma is important for global identification of β 0 .
We can similarly show that As B is compact, we can definê Then, by triangle inequality, it follows that for all β ∈ B, This concludes the proof.

Lemma S2. Under Assumptions 1 and 5, there is
Proof. (i) From Assumptions 1 and 5, we have (ii) From (S11), we have that is established in the proof of Lemma S1. This concludes the proof.

Proof. (i) Let
Under Assumption 2, we can write and an intermediate term ). We will show that when m 2 /n → 0, For A 1 , from the Strong Schwartz Matrix Inequality 1 , Next, we analyzeπ j − π j0 for j = 1, . . . , m. Asπ j = arg min π j Z j −Xπ j 2 , wherẽ Using Hölder's inequality and Cauchy-Schwartz inequality, we have On the other hand, with σ min (X) being the minimum singular value ofX, we have X (π j − π j0 ) 2 ≥ σ 2 min (X) π j − π j0 2 ≥ Cn π j − π j0 2 , where the last inequality is because ξ min (n −1X TX ) ≥ C from Assumption 4. In addition, from applying Lemma 8 in Chernozhukov et al. (2015) and the boundedness of Z and X, we have Then from Markov inequality, we know that max j=1,...,m X T (Z j −Xπ j0 ) ∞ = O p ( √ n log m).
Combining the above derivations, we have for j = 1, . . . , m, and thus,

This concludes the proof of
For A 2 , again using the Strong Schwartz Matrix Inequality, We use Tropp (2015, Theorem 1) and Markov inequality to construct a bound for A 21 . From Theorem 1 of Tropp (2015), we need to calculate the matrix variance parameter v = E(A 21 A T 21 ) and the large deviation where the third line is from V i being bounded, the last line is from Assumption 5. Also note that and the last inequality is from Bühlmann and Van De Geer (2011, Lemma 14.12). This implies that From the above analysis and Theorem 1 of Tropp (2015), we know the matrix variance parameter v is driving the upper bound and Then, from Markov inequality, we know that A 21 2 = O p (m log m/n). Using the same argument as the proof of π j −π j0 in (S13), we can show that μ−µ 0 = O p ( m log m/n).
Thus, A 2 = O p (m log m/n). The last term A 3 is bounded using the same argument.
Finally, as the remainder term Rem a consists of higher order terms, we can use the above arguments to show that Rem a is negligible. This concludes the proof of part (a).
For B 2 , we show that The other term in B 2 can be shown in the same way. Letγ . Using triangle inequality and Strong Schwartz Matrix Inequality, we write We can use the same argument as in the proof of (S14) to show that Since the remainder term Rem c contains the higher order terms, we conclude that (S15) is (ii) Because g i (β, η) is linear in β, G i (η) inherits its nice property in terms of the nuisance parameters. In other words, we still have This key property ensures that the claim in (ii) is true and can be proved in the same way as (i). The details are omitted.

The result follows from (i)-(ii) and the compactness of B.
A direct implication of the following result (i) is that when Next, under an extra assumption that where Q is a vector that includes all quadratic terms of X, i.e., they follow linear models that include a full set of quadratic terms, we can then estimate these nuisance parameters using least squares (denoted asω andθ) instead of nonparametric kernel. We prove the following result: Lemma S5': Suppose that Assumptions 2-5 hold. Also suppose that where Q is a vector that includes all quadratic terms of X. When m 2 /n → 0, Proof. (i) It suffices to prove part (b) in the proof of Lemma S5 (i). The rest are the same as the proof of Lemma S5. Writê Note that Q has a finite dimension, and thuŝ Here, B 1 = O p (n −1/2 ) by finite-dimensional linear model theory. For the second term, using the Strong Schwartz Matrix Inequality, Similar to the proof of Lemma S5, we use Tropp (2015, Theorem 1). We calculate the Also note that This implies that From the above analysis and Theorem 1 of Tropp (2015), we know the matrix variance parameter ν is driving the upper bound and From Markov inequality, we know that 1 We showed in the proof of Lemma S5, μ − µ 0 = O p ( m log m/n). Thus, B 2 = O p (m log m/n). The last term B 3 is bounded using the same argument. As the remainder term consists of higher order terms and can be shown to be negligible. Therefore, we havê Then, using a similar argument as that for A 2 in Lemma S5, we have Thus, we see that part (b) is negligible compared to part (a), so ĝ(β 0 ,η) −ĝ(β 0 , η 0 ) is asymptotically equivalent with part (a).
The proof of (ii) and (iii) are the same as that of Lemma S5.
Lemma S6. Under Assumptions 2-5, Then By triangle inequality, A key technique we use here is to vectorize the matrix, then a proof very similar to the proof of Lemma S5 will establish the result. Specifically, since I 1 ≤ I 1 F = vec(I 1 ) , where vec(A) = (a T 1 , . . . , a T m ) T ∈ R mn for a matrix A = (a 1 , . . . , a m ) ∈ R n×m . Let ⊗ be the Kronecker product. For vec(I 1 ), vec(I 1 ) For the above first term, from the Strong Schwartz Matrix Inequality, where the last expression is obtained using the same argument as in the derivation of A 21 in the proof of Lemma S5, and π −π 0 2 = O p (m log m/n) derived after (S13). The second term in (S16) has the same norm as the first term in (S16). For the third term in (S16), ≤ Cn −1 m and the large deviation parameter where the last inequality is from Assumption 4 and Bühlmann and Van De Geer (2011).
Hence, we have from Tropp (2015, Theorem 1) that This combined with Markov inequality, we have that m j=1 n i=1 S ij 2 = O p (m 2 log m/n).
This implies the third term in (S16) is The other proofs are very similar to the derivations in the proof of Lemma S5 and are omitted. Hence, we have shown that I 1 = O p (m 3/2 log m/n).
We similarly vectorize I 3 and analyze vec(I 3 ), following the steps in the proof of Lemma S5. First, we show that Following the notations in (S14). We obtain the linearization of 1 and the remainder term from the linearlization is o p ( m/n). Then, we have 1 n n i=1 L(O i ;γ− γ 0 ; µ) = O p (m/n). Thus, we have that (S17) is of order o p ( m/n). We can also show The other terms in I 3 can be shown similarly. Therefore, I 3 ≤ I 3 F = vec(I 3 ) = o p ( m/n) + o p ( m 3 log m/n).
Next, notice that I 2 can be bounded by where Rem I is a higher order term, the above first term can be bounded in the same way as I 1 , the above second term can be bounded in the same way as I 3 .
Combining the above arguments, we have that (ii)-(iii) The proof follows the same steps as in part (i) and is omitted.
The result follows from (i)-(iii) and the compactness of B.
(v) The proof follows the same steps as in part (iv).
Similar to Lemma S5', Lemma S6' is a parallel result to Lemma S6 assuming that where Q is a vector that includes all quadratic terms of X. When m 2 /n → 0, Lemma S6 (i)-(v) hold.
Proof. (i) The proof is similar to the proof of (i) in Lemma S6. It suffices to prove the term Hence, I 3 is negligible compared to the other terms, and Lemma S6 (i) still holds.
(ii)-(v) The proof is the same as those in Lemma S6.

Other Exposure and Outcome Types
Many MR applications consider binary outcomes (Holmes et al., 2014(Holmes et al., , 2017. Binary exposure is not very common in MR studies (Burgess and Labrecque, 2018), but are still of interest (Nead et al., 2015;Gage et al., 2017;Larsson et al., 2017;Vaucher et al., 2018). In this section, we extend the methods to consider binary exposure and/or binary outcome.
We focus on identification and leave formal treatment of inference under many weak invalid IVs to future work.
In this section, when the exposure (or outcome) variable is continuous, we consider the linear model (i.e., the identity link function); when the exposure (or outcome) variable is binary, we consider the log-linear model (i.e., the exponential link function). Therefore, for continuous outcome with f y (x) = x, β 0 a = E[Y | A = a, U, Z, X] − E[Y | A = 0, U, Z, X] encodes the treatment effect on the outcome mean upon increasing the exposure by a unit; for binary outcome with f y (x) = exp(x), β 0 a = log{P (Y = 1 | A = a, U, Z, X)} − log{P (Y = 1 | A = 0, U, Z, X)} encodes the log risk ratio. Other types of link function (e.g., logistic or probit) are not considered because of the noncollapsibility (Baiocchi et al., 2014) and the effect of U and A, Z are not easily separable (Clarke and Windmeijer, 2012).
Consider the following structural equations: where β 0 , α 0 , γ 0 are unknown true parameters, ξ y , ξ a are unspecified functions, f y , f a are pre-specified link functions, and Z ⊥ U |X.
Next, we state our identification results. Let R A = A − E(A | Z, X) when f a (x) = x; when f a (x) = exp(x). Let r Y = Y − βA when f y (x) = x; r Y = Y exp(−βA − α T Z) when f y (x) = exp(x); r Y 0 denote r Y when β = β 0 and α = α 0 . We use η to denote the nuisance parameters which may be different for each scenario and let η 0 be the true values.
Proposition 1 provides identification formulas for β 0 for binary exposure and/or binary outcome. The proof will show that E(R A r Y 0 |Z, X) = E(R A r Y 0 |X) holds almost surely in all cases and E(r Y 0 |Z, X) = E(r Y 0 |X) holds almost surely when f y (x) = exp(x). In all the cases, R A is the residual in A after netting out the effect of Z; r Y 0 is the residual in Y after netting out the effect of A when f y (x) = x, and is the residual in Y after netting out the effect of A and Z when f y (x) = exp(x). Note that when f a (x) = x, the exposure model (S25) can be relaxed, and identification in Proposition 1 remains true as long as E(A | U, Z, X) can be expressed as γ 0 (Z, X) + ξ a (U, X), where γ 0 , ξ a are unspecified functions. Similarly, when f y (x) = x, the outcome model (S24) can be relaxed, and identification in Proposition 1 remains true as long as E(Y |U, Z, X) can be expressed as β 0 A + α 0 (Z, X) + ξ y (U, X), where α 0 , ξ y are unspecified functions.
For binary exposure and continuous outcome (i.e., f a (x) = exp(x), f y (x) = x), γ 0 can be identified a prior from a separate set of estimation equations E[(Z − E(Z|X))R A ] = 0.
In Proposition 1(b) where the outcome is binary (i.e., f y (x) = exp(x)), α 0 cannot be identified a prior and needs to be identified simultaneously with β 0 , so that the estimation equations for β 0 and α 0 are stacked.
Proof. (a) The case where f a (x) = f y (x) = x is already established in the main article.
The rest of the proof follows the same step as the proof of continuous exposure and binary outcome.