On the instrumental variable estimation with many weak and invalid instruments

On the instrumental variable estimation with many weak and invalid instruments

\affil[a]{Department of Statistics, The Chinese University of Hong Kong} \affil[b]{Department of Statistics, University of Oxford} \affil[c]{Department of Economics, The Chinese University of Hong Kong}

{\bf{Keywords}}: {Invalid Instruments, Model Identification, Non-convex Penalty, Treatment Effect, Weak Instruments.}

Introduction

Recently, estimation of causal effects with high-dimensional observational data has drawn much attention in many research fields such as economics, epidemiology and genomics. The instrumental variable (IV) method is widely used when the treatment variable of interest is endogenous. As shown in Figure (ref), the ideal IV needs to be correlated with the endogenous treatment variable (C1), it should not have a direct effect on the outcome (C2) and should not be related to unobserved confounders that affect both outcome and treatment (C3).

Our research is motivated by the difficulty of finding IVs that satisfy all the above conditions. In applications, invalid IVs (violation of C2 or C3) davey2003mendelian,kang2016instrumental,windmeijer2019use and weak IVs { (concerning the weak correlation in C1)} Bound1995,staiger1997instrumental are prevalent. A strand of literature studies the “many weak IVs" problem Stock2002,chao2005consistent. With the increasing availability of large datasets, IV models are often high-dimensional Belloni2012sparse,lin2015,fan2018nonparametric, and have potentially weak IVs Andrews2019, and invalid IVs guo2018confidence,windmeijer2019confidence. Among those problems, we mainly focus on the invalid IV problem, while allowing for potential high-dimensionality and weak signals.

Related Works

Most related works fall into two main categories: robust estimation with invalid IVs and estimation which can select valid IVs without any prior knowledge of validity. The first strand of literature allows all IVs to be invalid. For example, kolesar2015identification restricted the direct effects of IVs on treatment and response being random effects and independent. In practice, this assumption might be difficult to justify. lewbel2012using,tchetgen2021genius,guo2022causal utilized conditional heteroskedasticity or heterogeneous curvatures to achieve robustness with potentially all IVs invalid. However, their performances are not satisfactory once the identification condition is not evident.

The second strand focused on unknown invalid IVs, while imposing certain identification conditions on the number of valid IVs. kang2016instrumental proposed a Lasso type estimator (sisVIVE). windmeijer2019use pointed out the inconsistent variable selection of sisVIVE under a relative IV strength condition and proposed an adaptive Lasso estimator, which has asymptotic oracle properties under the assumption that more than half of the IVs are valid, also called the majority rule. guo2018confidence,windmeijer2019confidence further developed two-step (the first one for relevance, the second one for validity) selection approaches, Two-Stage Hard Thresholding (TSHT) and Confidence Intervals IV (CIIV), respectively, under the plurality rule conditional on the set of relevant IVs. The plurality rule states that the valid IVs form the largest group. However, the approaches mentioned above are not robust to many weak IVs due to the restriction of the majority/plurality rule amongst the strong IVs instead of all IVs. Our method closely follows this strand of literature. Instead of a two-step selection, we require the plurality rule for valid IVs for a one-step selection procedure, thus considerably relaxing the requirement of valid IVs in theory and most practical scenarios.

The study of many (weak) IVs originated from empirical motivations but often assumed known validity. For example, staiger1997instrumental,hansen2008estimation,newey2009generalized,hansen2014instrumental considered different estimators for situations with many (weak) valid IVs but fixed the number of known covariates. kolesar2015identification,kolesar2018minimum allowed the number of covariates to grow with the sample size. { seng2022structural introduced an alternative model averaging method to handle weak IVs in low- and high-dimensional settings.} We consider the weak IV issues that are prevalent in empirical studies.

Main Results and Contributions

We propose a {\bf{W}}eak and {\bf{I}}nvalid IV robust {\bf{T}}reatment effect (WIT) estimator. The sparsest rule is sufficient for identification and is operational in numerical optimization. The proposed procedure has a selection stage (regarding IV validity) and a post-selection estimation stage. The selection stage is a penalized IV-regression via minimax concave penalty (MCP, zhang2010nearly), a proper surrogate penalty aligned with the identification condition to achieve model selection consistency of valid IVs under much weaker technical conditions than existing methods guo2018confidence, windmeijer2019confidence. In the estimation stage, we utilize the limited information maximum likelihood (LIML) estimator to handle the weak IVs staiger1997instrumental. An efficient computational algorithm for the optimal solution is provided. {{The computer codes for implementing the WIT estimator are available at \url{https://github.com/QoifoQ/WIT}.}}

The key contributions of this paper are summarized as follows.

The article is organized as follows. In Section (ref), we describe the model with some invalid IVs and analyze identification conditions in a general way. In Section (ref), we present the methodology and the novel WIT estimator. We establish the theorems to identify the valid IVs, estimation consistency, and asymptotic normality. Section (ref) shows the finite sample performance of our proposed estimator using comprehensive numerical experiments. Section (ref) applies our methods to { study the effect of BMI on diastolic blood pressure using Mendelian Randomization}. Section (ref) concludes. All the technical details and proofs are provided in the appendix.

Model and Identification Strategy

Potential Outcome Model with Some Invalid IVs

{ For $i =1,2,\ldots,n$, we have the random sample $(Y_i,D_i,\boldsymbol{Z}_{i.})$, where $Y_i\in \mathbb{R}^1$ is the outcome variable, $D_i\in\mathbb{R}^1$ is the (endogenous) treatment variable and $\boldsymbol{Z}_i\in\mathbb{R}^p$ are the potential IVs. Following the same model setting as in small2007sensitivity,kang2016instrumental,guo2018confidence,windmeijer2019confidence,windmeijer2019use, we consider a linear functional form between treatments $D_i$ and instruments $\boldsymbol{Z}_{i.}$ as the first-stage specification; meanwhile, a linear exposure of $Y_i$ and $D_i$ and $\boldsymbol{Z}_i$ is assumed as follows:

where $\epsilon_{i},\eta_{i}$ are random errors.

Following kang2016instrumental, we define the valid instruments as follows,

The validity of the $j$-th IV is quantified by $\alpha^*_j$, which captures the direct effect of the potential IV $\boldsymbol{Z}_j$ on the outcome (C1) as well as its influence on unmeasured confounders (C2). More details can be found in kang2016instrumental.} Further, we define the valid IV set $\mathcal{V}^* = \{j: \alpha^*_j = 0\}$ and invalid IV set $\mathcal{V}^{c*} = \{j: \alpha^*_j \neq 0\}$. Let ${p_{\mathcal{V}^*}} = |\mathcal{V}^*|$, ${p_{\mathcal{V}^{c*}}} = |\mathcal{V}^{c*}|$ and $p = {p_{\mathcal{V}^{*}}}+{p_{\mathcal{V}^{c*}}}$. Notably, ${p_{\mathcal{V}^*}} \geq 1$ refers to the existence of an excluded IV, thus satisfying the order condition wooldridge2010econometric. Let the $n \times p$ matrix of observations on the instruments be denoted by $\boldsymbol{Z}$, and the $n$-vectors of outcomes and treatments by $\boldsymbol{Y}$ and $\boldsymbol{D}$, respectively. We consider the cases of many and weak IVs in (ref) and make the following model assumptions:\\ {\bf{Assumption}} 1 (Many valid and invalid IVs): $p<n$, ${p_{\mathcal{V}^{c*}}}/n \rightarrow \upsilon_{p_{\mathcal{V}^{c*}}}+o(n^{-1/2})$ and ${p_{\mathcal{V}^*}}/n\rightarrow \upsilon_{p_{\mathcal{V}^*}}+o(n^{-1/2})$ for some non-negative constants $\upsilon_{p_{\mathcal{V}^{c*}}}$ and $\upsilon_{p_{\mathcal{V}^*}}$ such that $0\leq\upsilon_{p_{\mathcal{V}^*}}+\upsilon_{p_{\mathcal{V}^{c*}}}<1$.\\ {\bf{Assumption}} 2: Assume $\boldsymbol{Z}$ is standardized. It then has full column rank and $\|\boldsymbol{Z}_j\|^2_2\leq n$ for $j = 1,2,\ldots,p$.\\ {\bf{Assumption}} 3: Let $\boldsymbol{u}_{i}=\left(\epsilon_{i}, \eta_{i} \right)^{\top}$. $\boldsymbol{u}_i\mid \boldsymbol{Z}_i$ are i.i.d. and follow a multivariate normal distribution with mean zero and positive definite covariance matrix $ \boldsymbol{\Sigma} = \left( {

} \right) $. The elements of $\boldsymbol{\Sigma}$ are finite and $\sigma_{\epsilon,\eta}\neq 0$.\\ {\bf{Assumption}} 4 (Strength of valid IVs): The concentration parameter $\mu_n$ grows at the same rate as $n$, i.e., $\mu_n \coloneqq \boldsymbol{\gamma}^{*\top}_{\boldsymbol{Z}_{\mathcal{V}^*}}\boldsymbol{Z}_{\mathcal{V}^*}^{\top}M_{\boldsymbol{Z}_{\mathcal{V}^{c*}}}\boldsymbol{Z}_{\mathcal{V}^*}\boldsymbol{\gamma}^{*}_{\boldsymbol{Z}_{\mathcal{V}^*}}/\sigma_\eta^2 {\rightarrow} \mu_0 n$, for some $\mu_0>0$.

Assumption 1 is identical to the assumption of many instruments in kolesar2015identification,kolesar2018minimum. It relaxes the conventional many IVs assumptions bekker1994alternative,chao2005consistent that only allow the dimension of valid IVs ${p_{\mathcal{V}^*}}$ to grow with $n$. Also, it has not been considered in the literature on selecting valid IVs kang2016instrumental,guo2018confidence,windmeijer2019confidence,windmeijer2019use. Assumption 2 is standard for data preprocessing and scaling $\boldsymbol{Z}_j$. Assumption 3 follows guo2018confidence,windmeijer2019confidence to impose the homoskedasticity assumption and endogeneity of treatment $D_i$. {

} {Assumption 4 implies a strong identification condition in terms of the concentration parameter bekker1994alternative,newey2009generalized. In the fixed $p$ case, it indicates the presence of a constant coefficient $\gamma_j$, $\exists j$, and the rest of IVs could be weak staiger1997instrumental. Specifically, we model weakly correlated IVs as $\gamma = Cn^{-\tau}$, $0<\tau\le 1/2$, which is the “local to zero” case staiger1997instrumental. Essentially this is a mixture of constant $\gamma$-type and asymptotically diminishing $\gamma$-type instruments for fixed $p$. For $p\to \infty$ in the same order as $n$, we allow all the IVs to be weak in the “local to zero” case with specified rates. This IV strength assumption can be further weakened along the lines of hansen2008estimation to have weak identification asymptotics. In this paper, we focus on the individually weak (diminishing to zero as in staiger1997instrumental) signals model in high-dimensionality. Notice our model allows for much weaker individually weak IVs regardless of their validity (as long as the concentration parameter satisfies Assumption 4), unlike that of guo2018confidence. Nevertheless, the constant $\mu_0$ can be a small number to accommodate empirically relevant finite samples with many individually weak IVs.}

Identifiability of Model (ref)

The following moment conditions can be derived from Model (ref):

where $\boldsymbol{\Gamma}^* = E(\boldsymbol{Z}^\top\boldsymbol{Z})^{-1}E(\boldsymbol{Z}^\top\boldsymbol{Y})$ and $\boldsymbol{\gamma}^* = E(\boldsymbol{Z}^\top\boldsymbol{Z})^{-1}E(\boldsymbol{Z}^\top\boldsymbol{D})$, both are identified by the reduced form models. Without the exact knowledge about which IVs are valid, kang2016instrumental considered the identification of $(\boldsymbol{\alpha}^*,\beta^*)$ via the unique mapping of

Notice that the moment conditions ((ref)) consist of $p$ equations, but $(\boldsymbol{\alpha}^*,\beta^*) \in \mathbb{R}^{p+1}$ need to be estimated and is hence under-identified without further restrictions. kang2016instrumental proposed a sufficient condition, called majority rule (first proposed by han2008), such that ${p_{\mathcal{V}^*}} \geq \lceil{p}/{2} \rceil$, to identify the model parameters without any prior knowledge of the validity of individual IVs. However, the majority rule could be restrictive in practice. guo2018confidence further relaxed it to the plurality rule as follows:

which was stated as an “if and only if” condition of identification of $(\boldsymbol{\alpha}^*,\beta^*)$.

We re-examine the identifiability problem from the model DGP perspective. Given first-stage information: $\{\boldsymbol{D},\boldsymbol{Z},\boldsymbol{\gamma}^*\}$, without loss of generality, we denote the { true} DGP with $\{\beta^*, \boldsymbol{\alpha}^*, \boldsymbol{\epsilon}\}$ in ((ref)) as DGP $\mathcal{P}_0$ that generates $\boldsymbol{Y}$. Given this $\mathcal{P}_0$, for $ j \in \mathcal{V}^{c*}$, we { obtain an alternative representation}: $ \boldsymbol{Z}_j{\alpha}^*_j = \frac{\alpha^*_j}{\gamma_j^*}(\boldsymbol{D}-\sum_{l\neq j}\boldsymbol{Z}_l\gamma^*_l-\boldsymbol{\eta})$. Denote $\mathcal{I}_c=\{j\in \mathcal{V}^{c*}:c={\alpha^*_j}/{\gamma^*_j}\}$, where $c\neq 0$. For compatibility, we denote $\mathcal{I}_0 = \mathcal{V}^*$. Thus, we can reformulate $\boldsymbol{Y} = \boldsymbol{D}\beta^*+\boldsymbol{Z}\boldsymbol{\alpha}^*+\boldsymbol{\epsilon}$ in ((ref)) to:

where $\{\tilde{\beta}^{c},\tilde{\boldsymbol{\alpha}}^{c},\tilde{\boldsymbol{\epsilon}}^{c}\} = \{\beta^*+c,\boldsymbol{\alpha}^*-c\boldsymbol{\gamma}^*,\boldsymbol{\epsilon}-c\boldsymbol{\eta}\}$, for some $ j \in \mathcal{V}^{c*}$. Evidently, for different $c\neq 0$, it forms different DGPs $\mathcal{P}_c = \{\tilde{\beta}^{c},\tilde{\boldsymbol{\alpha}}^{c},\tilde{\boldsymbol{\epsilon}}^{c}\}$ that can generate the same $\boldsymbol{Y}$ (given $\boldsymbol{\epsilon}$), which also satisfies the moment condition ((ref)) as $\mathcal{P}_0$ since $E(\boldsymbol{Z}^\top\tilde{\boldsymbol{\epsilon}}^c) = \boldsymbol{0}$. Building on the argument of guo2018confidence, Theorem 1, the additional number of potential DGPs satisfying the moment condition ((ref)) is the number of distinguished $c \ne 0$ for $ j \in \mathcal{V}^{c*}$. We formally state this result in the following theorem.

Theorem (ref) shows there is a collection of model DGPs { (different parametrizations)}

corresponding to the same observation $\boldsymbol{Y}$ conditional on first-stage information. Given some $\mathcal{P}_0$, there are additional $1 \le G\le {p_{\mathcal{V}^{c*}}}$ equivalent DGPs. All members in $\mathcal{Q}$ are related through the transformation procedure ((ref)) and $1 <|\mathcal{Q}| = G+1\leq p$. Notably, the non-overlapping sparse structure among all possible DGPs leads to the sparsest model regarding $\boldsymbol{\alpha}^*$ being equivalent to plurality rule $|\mathcal{V}^*|> {\operatorname{max}}_{c \neq 0} |\mathcal{I}_{c}|$ in the whole set of IVs.

The Sparsest ($\boldsymbol{\alpha}$) Rule

The sparsest rule is conceptually equivalent to the plurality rule on the whole IV set, considering the non-overlapping sparse solutions given by Theorem (ref). To relax the majority rule, guo2018confidence proposed to use the plurality rule based on the relevant IV set:

where $\mathcal{S}^*$ is the set of strong IVs estimated by $\hat{\mathcal{S}}$ via first-step hard thresholding. Thus, TSHT and CIIV explicitly leverage the $\hat{\mathcal{S}}$-based plurality rule to estimate $\mathcal{V}^*_{\hat{\mathcal{S}}}$ and $\beta^*$.

In contrast to earlier studies on invalid IVs, our approach leverages the information from weak IVs. Firstly, weak IVs can be employed to estimate $\beta^*$. In situations where strong IVs are not available, weak IV robust estimators such as LIML prove useful Andrews2019. Secondly, weak IVs can aid in the identification of the valid IVs set, as demonstrated in Theorem (ref). When weak IVs are present, the plurality rule applied after the initial selection of strong instruments may yield unstable estimates of $\mathcal{V}^*$, as exemplified in the following scenario.

The mixture of weak and invalid IVs is ubiquitous in practice, especially in many IVs cases. For the sake of using all instruments' information for estimating $\beta^*$ and identification of $\mathcal{V}^*$, we allow some individual IV strengths to be local to zero chao2005consistent, say $\gamma^*_j\rightarrow 0$, or a small constant that cannot pass the first-stage threshold guo2018confidence unless with a very large sample size. However, we can see that in ((ref)), plurality rule-based methods that rely on first-stage selection are problematic, since $\mathcal{I}_0 = \{j :{\alpha^*_j}/{\gamma_j^*} = 0\}$ is ill-defined asymptotically due to the problem of “${0}/{0}$” if $\gamma_j^*$ is local to zero.

For using weak IVs information and improving finite sample performance, we turn to the sparsest rule that is also operational in computation algorithms. From the multiple DGPs $\mathcal{Q}$, recall $\mathcal{P}_c = \{\tilde{\beta}^{c},\tilde{\boldsymbol{\alpha}}^{c},\tilde{\boldsymbol{\epsilon}}^{c}\} = \{\beta^*+c,\boldsymbol{\alpha}^*-c\boldsymbol{\gamma}^*,\boldsymbol{\epsilon}-c\boldsymbol{\eta}\}$, where $\tilde{\boldsymbol{\alpha}}^c_{\mathcal{I}_c} = \boldsymbol{0}$. For other elements in $\tilde{\boldsymbol{\alpha}}^c$ (corresponding to a different DGP in $\mathcal{Q}$) and $j \in \{j:\alpha^*_j/\gamma^*_j = \tilde{c}\neq c\}$, we obtain

The above $|\tilde{\alpha}^c_j|$ needs to be distinguished from $0$ on the ground of the non-overlapping structure stated in Theorem (ref). To facilitate the discovery of all solutions in $\mathcal{Q}$, we assume: \\ {\bf{Assumption}} 5 (Separation Condition): $|\boldsymbol{\alpha}_{\mathcal{V}^{c*}}^*|_{\text{min}}>\kappa(n)$ and $|\tilde{\alpha}_j^c|>\kappa^{c}(n)$ for $j \in \{j:\alpha^*_j/\gamma^*_j = \tilde{c}\neq c\}$, where $\kappa(n)$ and $\kappa^c(n)$ are a generally vanishing rate specified by a particular estimator.

The conditions described above are comparable to the “beta-min" condition van2009conditions,loh2015regularized. { The rates $\kappa(n)$ and $\kappa^c(n)$ will be further specified in Theorem (ref). We aim to facilitate the understanding of these technical conditions by providing practical examples. In Section 3.3, we introduce two commonly encountered DGPs in both low and high dimensions. These examples demonstrate the validity of Assumption 5 without delving into the intricate specifics of the separation (beta-min) condition.

}

Hence, the identification condition known as the sparsest rule is formally defined as\\ {\bf{Assumption}} 6: {(The Sparsest Rule)}: $\boldsymbol{\alpha}^* = {\operatorname{arg min}}_{\mathcal{P}=\{\beta,\boldsymbol{\boldsymbol{\alpha}},\boldsymbol{\epsilon}\}\in \mathcal{Q}} \,\,\|\boldsymbol{\alpha}\|_0.$ \addtocounter{example}{-1}

In the next subsection, we reexamine the penalized approaches by kang2016instrumental and windmeijer2019use, and discuss a class of existing estimators concerning penalization, identification, and computation. We also explore the general penalization approach that aligns model identification with its objective function.

Penalization Approaches with Embedded Surrogate Sparsest Rule

A Lasso penalization approach was first used in the unknown IV validity context by kang2016instrumental. We extend this to a general formulation and discuss the properties of different classes of penalties.

Consider a general penalized estimator based on moment conditions (ref),

where $p_\lambda^{\text{pen}}(\boldsymbol{\alpha}) = \sum_{j=1}^{p}p_\lambda^{\text{pen}}(\alpha_j)$ and $p_\lambda^{\text{pen}}(\cdot)$ is a general penalty function with tuning parameter $\lambda > 0$ and ${p_\lambda^{\text{pen}}}^{\prime}(\cdot)$ is its derivative that satisfy: $\underset{x \rightarrow 0^+}{\operatorname{lim}}\, {p_\lambda^{\text{pen}}}^{\prime}(x) = \lambda$, $p_\lambda^{\text{pen}}(0) = 0$, $p_\lambda^{\text{pen}}(x) =p_\lambda^{\text{pen}}(-x) $, $(x-y)(p_\lambda^{\text{pen}}(x)-p_\lambda^{\text{pen}}(y))\geq0$, and ${p_\lambda^{\text{pen}}}^{\prime}(\cdot)$ is continuous on $(0,\infty)$.

In the RHS of ((ref)), (I) and (II) correspond to two requirements for the collection of valid DGPs in $\mathcal{Q}$ defined in ((ref)). $(I)$ is a scaled finite sample version of $E\Big((\boldsymbol{Z}^\top \boldsymbol{\epsilon})^\top (\boldsymbol{Z}^\top\boldsymbol{Z})^{-1}(\boldsymbol{Z}^\top \boldsymbol{\epsilon})\Big)$, which is a $(\boldsymbol{Z}^\top\boldsymbol{Z})^{-1}$ weighted quadratic term of condition $E(\boldsymbol{Z}^\top\boldsymbol{\epsilon}) = \boldsymbol{0}$, and $(II)$ is imposed to ensure sparsity structure in $\widehat{\boldsymbol{\alpha}}$.

Further, regarding $(I)$, one can reformulate ((ref)) with respect to $\widehat{\boldsymbol{\alpha}}^{\text{pen}}$ as

where $\tilde{\boldsymbol{Z}} = M_{\widehat{\boldsymbol{D}}}\boldsymbol{Z}$ and $\widehat{\boldsymbol{D}} = P_{\boldsymbol{Z}}\boldsymbol{D} = \boldsymbol{Z}\widehat{\boldsymbol{\gamma}}$, where $\widehat{\boldsymbol{\gamma}}$ is the least squares estimator of $\boldsymbol{\gamma}$, see kang2016instrumental. The design matrix $\tilde{\boldsymbol{Z}}$ is rank-deficient with rank $p-1$ since $\tilde{\boldsymbol{Z}}\widehat{\boldsymbol{\gamma}} = 0$. However, we show that it does not affect the $\boldsymbol{\alpha}$ support recovery using a proper penalty function. {{On the other hand, $\tilde{\boldsymbol{Z}}$ is a function of $\boldsymbol{\eta}, \boldsymbol{\gamma}^*$ and $\boldsymbol{Z}$, hence is correlated with $\boldsymbol{\epsilon}$. This inherited endogeneity initially stems from $\widehat{\boldsymbol{D}}$, in which $E(\widehat{\boldsymbol{D}}^\top \boldsymbol{\epsilon}) = \sigma_{\epsilon,\eta}^2 p/n$ does not vanish in the many IVs model (Assumption 1). The following lemma implies that the level of endogeneity of each $\tilde{\boldsymbol{Z}_j}$ is limited. }}

Concerning $(II)$ in ((ref)), Theorem (ref) shows that model (ref) can be identified by different strategies with non-overlapping results. On the ground of the sparsest rule assumption, the role of the penalty on $\boldsymbol{\alpha}$, i.e.\ $p_\lambda^{\text{pen}}(\boldsymbol{\alpha})$, should not only impose a sparsity structure but also serve as an objective function corresponding to the identification condition we choose. For example, the penalty $\lambda \|\boldsymbol{\alpha}\|_0$ matches the sparsest rule.

To see the roles of a proper penalty function clearly, we rewrite ($\ref{problem1}$) into an equivalent constrained objective function with the optimal penalty $\|\boldsymbol{\alpha}\|_0$ regarding the sparsest rule:

where $\delta$ is the tolerance level which we specify in Section (ref). The constraint above narrows the feasible solutions into $\mathcal{Q}$ because it aligns with the Sargan test statistic, sargan1958estimation, $\|P_{\boldsymbol{Z}}(\boldsymbol{Y}-\boldsymbol{D}\beta-\boldsymbol{Z}\boldsymbol{\alpha})\|_2^2/\|(\boldsymbol{Y}-\boldsymbol{D}\beta-\boldsymbol{Z}\boldsymbol{\alpha})/\sqrt{n}\|_2^2 = O_p(1)$ under the null hypothesis $E(\boldsymbol{Z}^\top \boldsymbol{\epsilon}) = \boldsymbol{0}$ as required in $\mathcal{Q}$; otherwise, the constraint becomes $O_p(n)$ that cannot be bounded by $\delta$. Thus, a properly chosen $\delta$ in (ref) leads to an equivalent optimization problem

that matches the identification condition in Assumption 6. Therefore, the primary optimization object in ((ref)) should also serve as an identification condition: the sparsest rule.

{

}

Due to computational NP-hardness for $\|\boldsymbol{\alpha}\|_0$ in ((ref)), a surrogate penalty function is needed. kang2016instrumental proposed to replace the optimal $l_0$-norm with Lasso in ((ref)), denoting their estimator sisVIVE as $(\widehat{\boldsymbol{\alpha}}^{\text{sis}},\hat{\beta}^{\text{sis}})$. And $\mathcal{V}^*$ is estimated as $\hat{\mathcal{V}}^{\text{sis}} = \{j:\widehat{\alpha}_j^{\text{sis}} = 0\}$. However, the surrogate $\ell_1$ penalty brings the following issues: { (a) Failure in consistent variable selection under some deterministic conditions, namely the sign-aware invalid IV strength (SAIS) condition windmeijer2019use:

The sisVIVE cannot achieve $\mathcal{P}_0$ under the SAIS condition, which is likely to hold when the invalid IVs are relatively stronger than valid ones. (b) kang2016instrumental, Theorem 2, proposed a non-asymptotic error bound $|\hat{\beta}^{\text{sis}}-\beta^*|$ for sisVIVE. The dependence of the restricted isometry property (RIP) constant of $P_{\widehat{\boldsymbol{D}}}\boldsymbol{Z}$ and the error bound is not clear due to the random nature of $\widehat{\boldsymbol{D}}$. (c) The objective function deviates from the original sparsest rule: $g_1(\mathcal{P}) = \|\boldsymbol{\alpha}\|_0$ and $g_2(\mathcal{P}) = \|\boldsymbol{\alpha}\|_1$ correspond to incompatible identification conditions unless satisfying an additional strong requirement

which further impedes sisVIVE in estimating $\beta^* \in \mathcal{P}_0$. }

Problems (a) and (b) relate to the Lasso problem within invalid IVs, while (c) exposes a deeper issue beyond Lasso: a proper surrogate penalty in ((ref)) must align with the identification condition.

windmeijer2019use suggested the Adaptive Lasso zou2006adaptive, using as an initial estimator the median estimator of han2008 to tackle problem (a)'s SAIS issue. This solution also addresses (c) concurrently. However, it necessitates a more demanding majority rule and requires all IVs to be strong in the fixed $p$ case. Like TSHT and CIIV, it is sensitive to weak IVs.

The following proposition outlines the appropriate surrogate sparsest penalty.

The surrogate sparsest penalty requirement aligns with the folded-concave penalized method fan2001variable,zhang2010nearly. We adopt MCP penalty for our proposed approach in the following. In standard sparse linear regression, MCP necessitates $\lambda = \lambda(n) = O(\sqrt{\log p/n})$ and ${p_\lambda^{\text{MCP}}}^\prime (t)= 0$ for $t>C\lambda$, with $C$ a constant, fulfilling Proposition (ref). We further elaborate that this property is applicable in invalid IVs scenarios in the following section, and demonstrate its ability to bypass sisVIVE's shortcomings shown in (a) and (b).

WIT Estimator

Estimation Procedure

We implement the penalized regression framework ((ref)) and specifically employ a concave penalty in ((ref)), the MCP, which provides nearly unbiased estimation. The MCP is favored numerically for its superior convexity in penalized loss and its consistent variable selection property without imposing incoherence conditions on the design matrix loh2017support,feng2019sorted. This makes MCP more suited to two-stage estimation problems than Lasso.

The selection stage is formally defined as:

where $\widetilde{\boldsymbol{Z}} = M_{\widehat{\boldsymbol{D}}}\boldsymbol{Z}$, and $\widehat{\boldsymbol{D}} = P_{\boldsymbol{Z}}\boldsymbol{D} = \boldsymbol{Z}\widehat{\boldsymbol{\gamma}}$ are the same as in (ref). $p_\lambda^{\text{MCP}}(\boldsymbol{\alpha}) = \sum_{j = 1}^p p_\lambda^{\text{MCP}}({\alpha}_j) = \sum_{j = 1}^p\int_0^{|\alpha_j|} \Big(\lambda-{t}/{\rho}\Big)_{+}dt$ is the MCP penalty and $\rho>1$ is the tuning parameter, which also controls convexity level $1/\rho$, and its corresponding derivative is ${p_\lambda^{\text{MCP}}}^{\prime}(t) = (\lambda-{|t|}/{\rho})_+$. Unlike Lasso, MCP imposes no penalty when $|\alpha_j|>\lambda\rho$, resulting in nearly unbiased coefficient estimation and exact support recovery without the SAIS condition (see (ref) for further discussion). Consequently, a consistent estimation of the valid IVs set, i.e., $\hat{\mathcal{V}} = \{ j:\widehat{{\alpha}}^{\text{MCP}}_j = 0 \}$ and $\Pr(\hat{\mathcal{V}} = \mathcal{V}^*) \overset{p}{\rightarrow}1$, is expected to hold under more relaxed conditions. We next illustrate how the WIT blends the benefits of penalized TSLS and LIML estimators at various stages.

The LIML estimator is consistent not only in classic many (weak) IVs bekker1994alternative,hansen2008estimation, but also in many IVs and many included covariates kolesar2015identification,kolesar2018minimum. However, simultaneous estimation of $\hat{\kappa}_{\text{liml}}$ and $\hat{\mathcal{V}^c}$ in (ref) below is difficult. In the selection stage (ref), we use the moment-based objective function. If we do not consider the penalty term (II), the moment-based part (I) of (ref) coincides with TSLS. Furthermore, the bias in TSLS has a limited effect on consistent variable selections (see Theorem (ref)). In the estimation step, however, due to LIML's superior finite sample performance and the issue of TSLS in the presence of many (or weak) IVs even when $\mathcal{V}^*$ is known sawa1969exact,chao2005consistent, we embed the LIML estimator to estimate $\beta^*$ on the basis of estimated valid IVs set via ((ref)). The performance of oracle-TSLS shown in simulations verifies this choice.

Consequently, we proposed the Weak and some Invalid instruments robust Treatment effect (WIT) estimator in the estimation stage,

Note, ((ref)) belongs to the general $k$-class estimators nagar1959bias, whose properties vary upon the choice of $\hat{\kappa}$, i.e.\ $\hat{\kappa} = 0$ refers to OLS and $\hat{\kappa} = 1$ reduces to the TSLS estimator. ((ref)) has a closed-form solution: $\hat{\kappa}_{\text{liml}} = \lambda_{\text{min}}\Big(\{[\boldsymbol{Y},\boldsymbol{D}]^\top M_{\boldsymbol{Z}}[\boldsymbol{Y},\boldsymbol{D}]\}^{-1}\{[\boldsymbol{Y},\boldsymbol{D}]^\top M_{\boldsymbol{Z}_{\hat{\mathcal{V}}^c}}[\boldsymbol{Y},\boldsymbol{D}]\}\Big)$, where $\lambda_{\text{min}}(\cdot)$ denotes the smallest eigenvalue. Focusing on $\hat{\beta}^{\text{WIT}}$ as the primary interest, we reformulate ((ref)) and ((ref)) based on the residuals of $\boldsymbol{Y},\boldsymbol{D},\boldsymbol{Z}_{\hat{\mathcal{V}}}$ on $\boldsymbol{Z}_{\hat{\mathcal{V}}^{c}}$. Denote $\boldsymbol{Y}_{\perp} = M_{\boldsymbol{Z}_{\hat{\mathcal{V}}^c}}\boldsymbol{Y}$, $\boldsymbol{D}_{\perp} = M_{\boldsymbol{Z}_{\hat{\mathcal{V}}^c}}\boldsymbol{D}$ and $\boldsymbol{Z}_{\hat{\mathcal{V}}\perp} = M_{\boldsymbol{Z}_{\hat{\mathcal{V}}^c}}\boldsymbol{Z}_{\hat{\mathcal{V}}}$ and notice $M_{\boldsymbol{Z}_{\hat{\mathcal{V}}^c}}M_{\boldsymbol{Z}_{\hat{\mathcal{V}}\perp}} = M_{\boldsymbol{Z}}$, thus it is equivalent to derive asymptotic results on the following model (ref).

Asymptotic Behavior of WIT Estimator

Throughout this section, we aim to recover the one specific element in $\mathcal{Q}$, denoted as $({\beta^*},\boldsymbol{\alpha}^*,\boldsymbol{\epsilon})$. Though a slight abuse of notation, we use $\widehat{\boldsymbol{\alpha}}$ to denote a local solution of ((ref)) with MCP.

All local solutions of ((ref)) we consider are characterized by the Karush–Kuhn–Tucker (KKT) or first-order condition, i.e.\

Explicitly, to the end of finding valid IVs via comparing with true signal $\boldsymbol{\alpha}^*$, we rewrite ((ref)) as

where the inequalities in the first line stem from the convexity of the MCP penalty and in the last line originate in the sub-derivative of the MCP penalty at the origin.

As discussed in Section (ref), $\tilde{\boldsymbol{Z}}$ is a function of $(\boldsymbol{Z},\boldsymbol{\gamma}^*,\boldsymbol{\eta})$ and thus endogenous with $\boldsymbol{\epsilon}$. The fact that $\tilde{\boldsymbol{Z}}$ inherits the randomness of $\boldsymbol{\eta}$ distinguishes itself from the general assumptions put on the design matrix, and obscures the feasibility of conditions required to achieve exact support recovery in the literature of penalized least squares estimator (PLSE) feng2019sorted,loh2017support,zhang2012general. The sisVIVE method imposed the RIP condition directly on $\tilde{\boldsymbol{Z}}$ to establish an error bound. However, the restricted eigenvalue (RE) condition bickel2009simultaneous is the weakest condition van2009conditions available to guarantee rate minimax performance in prediction and coefficient estimation, as well as to establish variable selection consistency for Lasso penalty. feng2019sorted further adopted the RE condition for non-convex penalty analysis. We then state the conditions on the design matrix $\tilde{\boldsymbol{Z}}$ of (ref) in the following. Define restricted cone $\mathscr{C}(\mathcal{V}^{*} ; \xi)=\left\{\boldsymbol{u}:\left\|\boldsymbol{u}_{\mathcal{V}^{*}}\right\|_{1} \leq \xi\left\|\boldsymbol{u}_{\mathcal{V}^{c*}}\right\|_{1}\right\}$ for some $\xi >0$ that estimation error $\widehat{\boldsymbol{\alpha}}-\boldsymbol{\alpha}^*$ belongs to. The restricted eigenvalue $K_{\mathscr{C}}$ for $\tilde{\boldsymbol{Z}}$ is defined as $K_{\mathscr{C}}$ = $K_{\mathscr{C}}(\mathcal{V}^*,\xi)\coloneqq \underset{\boldsymbol{u}}{\operatorname{inf}}\{\|\tilde{\boldsymbol{Z}}\boldsymbol{u}\|_2/(\|\boldsymbol{u}\|_2 n^{1/2}): \boldsymbol{u}\in \mathscr{C}(\mathcal{V}^* ; \xi) \}$ and the RE condition refers to the condition that $K_{\mathscr{C}}$ for $\tilde{\boldsymbol{Z}}$ should be bounded away from zero.

Lemma (ref) elaborates that the RE condition on $\tilde{\boldsymbol{Z}}$ holds without any additional assumptions on $\tilde{\boldsymbol{Z}}$, unlike the extra RIP condition for sisVIVE. Moreover, this restricted cone is invariant to scaling, thus indicating the accommodation of many weak IVs. These two features suggest the theoretical advantages of penalized methods ((ref)) over existing methods.

Next, we discuss the selection of valid IVs by comparing the local solution of ((ref)) with the oracle (moment-based) counterpart. Define $\widehat{\boldsymbol{\alpha}}^{\text{or}}_{\mathcal{V}^*} = \boldsymbol{0}$ and

where $\hat{\beta}_{\text{or}}^{\text{TSLS}} = [\boldsymbol{D}^\top(P_{\boldsymbol{Z}}-P_{\boldsymbol{Z}_{\mathcal{V}^{c*}}})\boldsymbol{D}]^{-1}[\boldsymbol{D}^\top(P_{\boldsymbol{Z}}-P_{\boldsymbol{Z}_{\mathcal{V}^{c*}}})\boldsymbol{Y}]$ and the two versions of $\widehat{\boldsymbol{\alpha}}^{\text{or}}_{\mathcal{V}^{c*}}$ are equivalent. Notice this $\hat{\beta}^{\text{TSLS}}_{\text{or}}$ is not for the final treatment effect estimation, but to illustrate the selection stage consistency only. To this end, we show the supremum norm of $\boldsymbol{R}^{\text{or}} = \widetilde{\boldsymbol{Z}}^{\top}(\boldsymbol{Y}-\widetilde{\boldsymbol{Z}}{\widehat{\boldsymbol{\alpha}}}^{\text{or}})/n$ is bounded by an inflated order of $O(\sqrt{\log p_{\mathcal{V}^*} /n})$. Denote $\Tilde{\Tilde{\boldsymbol{D}}} = ( P_{\boldsymbol{Z}}-P_{\boldsymbol{Z}_{\mathcal{V}^{c*}}})\boldsymbol{D}$ and $\Tilde{\Tilde{\boldsymbol{Q}}}_n = \boldsymbol{Z}_{\mathcal{V}^*}^\top ( P_{\boldsymbol{Z}}-P_{\boldsymbol{Z}_{\mathcal{V}^{c*}}})\boldsymbol{Z}_{\mathcal{V}^*}/n$, we derive the following lemma.

Based on Lemmas (ref) and (ref), we consider the set $\mathscr{B}(\lambda,\rho) = \{\widehat{\boldsymbol{\alpha}} \text{ in} \,\,(\ref{form 2}): \lambda \geq \zeta, \rho > K^{-2}_{\mathscr{C}}(\mathcal{V}^*,\xi)\lor 1 \}$, in which $\zeta$ is defined in ((ref)) and $\xi$ is guaranteed by Lemma (ref), as a collection of all local solutions $\widehat{\boldsymbol{\alpha}}$ computed in ((ref)) through a broad class of MCP under a certain penalty level $\lambda$ and convexity $1/\rho$. Given that the computed local solutions in practice are through a discrete path with some starting point (see Section (ref)), we further consider the computable solution set $\mathscr{B}_0(\lambda,\rho)$, introduced by feng2019sorted, i.e.,

The connection between $\mathscr{B}_0(\lambda,\rho)$ and $\mathscr{B}(\lambda,\rho)$ is that $\exists \,\, \widehat{\boldsymbol{\alpha}}^{(l)}\in \mathscr{B}(\lambda,\rho)$ with penalty level $\lambda^{(l)}$ increasing with the index $l = 1,2,\ldots$, such that $\widehat{\boldsymbol{\alpha}}^{(0)}-\boldsymbol{\alpha}^*\in \mathscr{C}(\mathcal{V}^*,\xi)$, $\widehat{\boldsymbol{\alpha}} = \widehat{\boldsymbol{\alpha}}^{(l)}$ for large enough $l$ and $\|\widehat{\boldsymbol{\alpha}}^{(l)}-\widehat{\boldsymbol{\alpha}}^{(l-1)}\|_1<a_0 \lambda^{(l)}$, where $a_0$ is specified in Lemma (ref) of (ref). Thus, $\mathscr{B}_0(\lambda,\rho)$ is a collection of approximations of $\boldsymbol{\alpha}$ in all DGPs.

Denote

Also, let $\Tilde{\Tilde{\boldsymbol{Q}}}_n^c$ and $\operatorname{Bias}(\hat{\beta}_{\text{or}}^{c~\text{TSLS}})$ be defined as $\mathcal{P}_c$ version of $\Tilde{\Tilde{\boldsymbol{Q}}}_n$ and $\operatorname{Bias}(\hat{\beta}_{\text{or}}^{\text{TSLS}})$. Then, we provide the asymptotic result of selection consistency of WIT.

In Theorem (ref), $T_1$ is similar to a standard rate $\sqrt{\log p /n}$ in penalized linear regression, while $T_2$ and $T_3$ are the additional terms that only play a role in the many IVs context and vanish fast in the finite strong IVs case (see Corollary (ref)). This result is new to the literature.

Proposition (ref) shows that $T_2$ is limited in the general case where dominant scaled $\gamma^*_j$ does not exist. For example, if we assume $\boldsymbol{Q}_n = \boldsymbol{I}$ and $\boldsymbol{\gamma}_{\mathcal{V}^*}^* = C\boldsymbol{1}_{p_{\mathcal{V}^*}}$, where $C$ is a constant or diminishing to zero, then $\|\boldsymbol{\gamma}_{\mathcal{V}^*}^*\|_\infty/\|\boldsymbol{\gamma}_{\mathcal{V}^*}^*\|_1 = o(\|\boldsymbol{\gamma}_{\mathcal{V}^*}^*\|_1/p_{\mathcal{V}^*}) = o(Cp_{\mathcal{V}^*}/p_{\mathcal{V}^*}) = o(C)$ holds and it follows that $T_2\rightarrow0$.

The following theorem describes the asymptotic behavior of the WIT estimator for many valid and invalid IVs cases by combining Theorem 1 and invariant likelihood arguments in kolesar2018minimum. We further denote two statistics

as the estimates of the covariance matrix of reduced-form error $\boldsymbol{\Omega} = \operatorname{Cov}(\boldsymbol{\epsilon}+\beta^*\boldsymbol{\eta},\boldsymbol{\eta})$ and a variant of concentration parameter, respectively. Also, let $m_{\text{max}} = \lambda_{\text{max}}( \boldsymbol{S}^{-1}\boldsymbol{T})$, $\hat{\mu}_n = \operatorname{max}(m_{\text{max}}-{p_{\mathcal{V}^*}}/n,0)$ and $\widehat{\boldsymbol{\Omega}} = \frac{n-p}{n-{p_{\mathcal{V}^{c*}}}/n} \boldsymbol{S}+\frac{n}{n-{p_{\mathcal{V}^{c*}}}/n}(\boldsymbol{T}-\frac{\hat{\mu}_n}{\hat{\boldsymbol{a}}^\top \boldsymbol{S}^{-1} \hat{\boldsymbol{a}}} \hat{\boldsymbol{a}} \hat{\boldsymbol{a}}^\top )$, where $\hat{\boldsymbol{a}} = (\hat{\beta}^{\text{WIT}},1)$ and {{ $|\boldsymbol{\Sigma}|$ is the determinant of $\boldsymbol{\Sigma}$}}.

Notably, when the number of invalid IVs ${p_{\mathcal{V}^*}}$ is a constant, the variance estimator above is reduced to the one that bekker1994alternative derived for the typical many IVs case. hansen2008estimation showed that it is still valid under many weak IVs asymptotics.

Special Cases in Low & High Dimensions

{ The values of $\kappa(n)$ and $\kappa^c(n)$ outlined in Theorem (ref) provide the general formulae to confirm Assumption 5. In this subsection, we discuss some representative cases in low and high dimensions and verify that Assumption 5 holds.

Finite $p$ Case

First, we show that the WIT estimator is a more powerful tool than the existing methods requiring the majority rule kang2016instrumental,windmeijer2019use also under the finite number of IVs with a mixture of strong and weak IVs settings. WIT achieves the same asymptotic results as windmeijer2019use under more relaxed conditions. Specifically, under finite IVs, Assumption 1 can be reduced to Assumption 1$^\prime$ as follows.

{\bf{Assumption}} 1$^\prime$ (Finite Number of IVs): ${p_{\mathcal{V}^{c*}}}\geq 1$ and ${p_{\mathcal{V}^*}}\geq1$ are fixed constants, and ${p_{\mathcal{V}^{c*}}}+{p_{\mathcal{V}^*}} = p<n$.

{In the finite IVs case, the $T_2$ and $T_3$ terms in Theorem (ref) for selecting valid IVs go to zero fast and the required separation rate reduces to $\kappa(n)\asymp n^{-1/2}$. We present the asymptotic properties for the WIT estimator here. Consider the following mixed IV strengths case. Let $\gamma_j^* = Cn^{-\tau_j}$ for $j = 1,2,\ldots,p$, $\tau_{\mathcal{V}^*} = {\operatorname{argmax}}_{\tau_j}\{\tau_j:j\in\mathcal{V}^*\}$, $\tau_{\mathcal{I}_c} = {\operatorname{argmax}}_{\tau_j}\{\tau_j:j\in\mathcal{I}_c\}$, and $\tau_{\mathcal{I}_{\tilde{c}}} = {\operatorname{argmax}}_{\tau_j}\{\tau_j:j\in\mathcal{I}_{\tilde{c}}\}$, where $c\neq \tilde{c}$.

} Corollary (ref) addresses a scenario where the violation of validity, $\alpha^*_j$, is held constant and thus is of the same order when $\tau_j=0$, or exceeds its strength parameter, $\gamma^*_j$, for $j \in \mathcal{V}^*_c$. This scenario frequently arises when conducting robust estimations, particularly in circumstances that inadvertently incorporate weakly relevant but strongly invalid IVs, or when there is no prior information on the candidate IVs set. Dealing with these situations effectively is a crucial component of robust statistical analysis.

Conversely, as noted by an anonymous reviewer, experienced researchers often select IVs that are argued to be both relevant and valid. A violation of IV validity could happen among relatively strong IVs. This corresponds to a scenario where the ratio of $\alpha^*_j/\gamma^*_j$ is small. In the following proposition, we demonstrate that the WIT estimator can effectively manage cases where $\alpha^*_j/\gamma^*_j = o(1)$.

Proposition (ref) gives the desired asymptotic results for the WIT estimator for DGPs in the $\alpha^*_j/\gamma^*_j=o(1)$ scenario, which aligns with the aforementioned case that some researchers may pick relatively strong IVs in the design phase, and true $\alpha^*_j$ to satisfy the separation condition in Assumption 5. It only requires that $|\gamma_l^*|>\underset{j\in \mathcal{V}^{c*}}{\operatorname{inf}}O(n^{-1/2} \gamma_j^* /\alpha_j^*)$, which means the valid IVs' strength signals are stronger than the rate of $n^{-1/2}$ over $\alpha^*_j/\gamma_j^*$. It can be easily implied from this condition that a smaller value in ratio of $\alpha^*_j$ to $\gamma^*_j$ helps to reduce the needed strength signal $\gamma_l^*$ of the valid IVs. Therefore, $\gamma^*_l = o(\gamma^*_j)$ for $l \in \mathcal{V}^*$ and $j \in \mathcal{V}^{c*}$ is a sufficient condition to verify it. For the second part in Assumption 5: $|\tilde{\alpha}_j^c|>\kappa^{c}(n)$ for $j \in \{j:\alpha^*_j/\gamma^*_j = \tilde{c}\neq c\}$ will hold automatically under $\alpha^*_j/\gamma^*_j = o(1)$. Table (ref) subsequently illustrates a simplified representative example featuring only 4 IVs, two of which are valid, and demonstrates the rates of $\kappa(n)$ and $\kappa^{c}(n)$ in true and transformed DGPs.

As a result, Corollary (ref) and Proposition (ref) provide easily interpretable examples of common DGPs across two distinct scenarios. These examples affirm the validity of Assumption 5, facilitating a more straightforward understanding for practitioners.

High Dimensional (Many IVs) Case

Next, we consider a common scenario Andrews2019,chao2005consistent with a large number (can grow with $n$) of individually weak instruments, where $\gamma^*_j = O(\sqrt{\log p/n^{1-\delta}})$ for $j = 1,2,\ldots,p$, and some small enough $\delta>0$ such that $\log n/n^{-\delta}<\infty$ to avoid an exploded variance in the outcome variable $Y_i$. Regarding invalid IVs, we assume the same rate of $|\boldsymbol{\alpha}_{\mathcal{V}^{c*}}^*|_{\text{min}} = O(\sqrt{\log p/n^{1-\delta}})$.

Consequently, the separation condition $\kappa(n)$ is reduced to:

as per Proposition (ref). This satisfies the first part of Assumption 5, that is, $|\boldsymbol{\alpha}_{\mathcal{V}^{c*}}^*|_{\text{min}}>\kappa(n)$. The second part of Assumption 5, $|\tilde{\alpha}_j^c|>\kappa^{c}(n)$ for $j \in \{j:\alpha^*_j/\gamma^*_j = \tilde{c}\neq c\}$, can be confirmed according to (ref),

where $|\tilde{c}-c|$ has the same order of $\alpha^*_j/\gamma^*_j \asymp O(1)$ for $j \in \mathcal{V}_c^* $. The threshold of $\kappa^c(n)$ in (ref) is now specified as

Therefore, the second part of Assumption 5 holds as well.

We summarize the above discussions as the following proposition.

}

Computational Implementation of WIT Estimator

Through Proposition (ref), we know that MCP belongs to the surrogate sparsest penalty under Assumptions 1-6, and it ensures the global solution in ((ref)) matches the sparsest rule. However, each element in $\mathcal{Q}$ is the local solution of ((ref)) and we can only obtain one local solution from one initial value practically. A multiple starting points strategy is needed to achieve the global solution. Enumerating the whole $\boldsymbol{\alpha}^*\in \mathbb{R}^p$ is impossible. Therefore, it is important to develop an efficient algorithm for the MCP penalty.

In light of practical use, we adopt the iterative local adaptive majorize-minimization (I-LAMM) algorithm fan2018lamm, which satisfies ((ref)) as shown in feng2019sorted Section 2.1, with different initial values to achieve the local solution of $\widehat{\boldsymbol{\alpha}}$ in ((ref)). See more technical details and derivation in (ref).

Motivated by the individual IV estimator windmeijer2019confidence such that $\hat{\beta}_j = \hat{\Gamma}_j/\hat{\gamma}_j\overset{p}{\rightarrow}\beta^*+\alpha^*_j/\gamma^*_j$, where $\hat{\boldsymbol{\Gamma}} = (\boldsymbol{Z}^\top \boldsymbol{Z})^{-1}\boldsymbol{Z}^\top \boldsymbol{Y}$ and $\hat{\boldsymbol{\gamma}} = (\boldsymbol{Z}^\top \boldsymbol{Z})^{-1}\boldsymbol{Z}^\top \boldsymbol{D}$, we can construct

with any initial value $\check{\beta}$. Thus $\| \check{\boldsymbol{\alpha}}(\check{\beta})-\boldsymbol{\alpha}^*\|_1 \overset{p}{\rightarrow}|\beta^*-\check{\beta}| \cdot\|\hat{\boldsymbol{\gamma}}\|_1$. It is valid to replace $(\beta^*,\boldsymbol{\alpha}^*)$ by $(\tilde{\beta}^c,\tilde{\boldsymbol{\alpha}}^c)$, we know that $\| \check{\boldsymbol{\alpha}}(\check{\beta})-\tilde{\boldsymbol{\alpha}}^c\|_1$ is asymptotically controlled by $|\tilde{\beta}^c-\check{\beta}| \cdot\|\hat{\boldsymbol{\gamma}}\|_1$. Thus, varying $\check{\beta}\in\mathbb{R}^1$ is equivalent to varying $\check{\boldsymbol{\alpha}}(\check{\beta})\in {\mathbb{R}}^{p}$ of a close DGP $\mathcal{Q}$.

Here we provide a simple heuristic procedure to efficiently choose a proper $\check{\beta}$. Let $\hat{\beta}_{[j]}$ be the sorted $\hat{\beta}_j$. A fuzzy MCP regression is conducted to clustering $\hat{\beta}_j$:

where $\bar{\boldsymbol{\beta}}\in \mathbb{R}^p$, $\bar{\lambda}$ is a prespecified penalty level and $\grave{\boldsymbol{\beta}}$ is initialized by $\grave{{\beta}}_j = \hat{\beta}_{[j]}$. Therefore, $\bar{\boldsymbol{\beta}}$ consists of less than $p$ distinct values. Thus we can choose the initial $\widehat{\boldsymbol{\alpha}}$ in (ref) such that

where $\check{\beta}$'s are chosen in $\bar{\boldsymbol{\beta}}$ with the priority given to the largest (remaining) cluster and $\widehat{{\alpha}}_j^{(0)}(\check{\beta}) = 0$ for $j$ to be the unsorted index of $\hat{\beta}_j$ such that $\check{\beta} = \bar{{\beta}}_j$.

Following (ref), which provides a theoretical guideline for the tuning parameter, we look for a data-driven tuning procedure that has good performance in practice. From numerical studies, results are found to not be sensitive to the choice of $\rho$, which is fixed to $2$ for most applications. However, $\lambda$ is required to be tuned. Cross-validation is implemented in sisVIVE, but is known to be time-consuming and to select too many IVs as invalid. windmeijer2019use,windmeijer2019confidence proposed to use the Sargan test under low dimensions to choose the tuning parameter consistently and obtain a good finite sample performance. However, the Sargan test is designed for fixed $p$ and cannot handle many IVs. Therefore, we propose the modified Cragg-Donald (MCD, kolesar2018minimum) test-based tuning procedure, which extends the Sargan test to allow high-dimensional covariates and IVs.

Specifically, consider a local solution $\widehat{\boldsymbol{\alpha}}$ in (ref), Denote ${{p_{\widehat{\mathcal{V}}}}} = |\{j:\widehat{{\alpha}}_j = 0\}|$ and ${{p_{\widehat{\mathcal{V}}^{c}}}} = |\{j:\widehat{{\alpha}}_j \neq 0\}|$. Let $m_{\text{min}}$ be the minimum eigenvalue of $\boldsymbol{S}^{-1}\boldsymbol{T}$, where $\boldsymbol{S}$ and $\boldsymbol{T}$ are defined as a function of $\widehat{\boldsymbol{\alpha}}$ in (ref). Then the MCD test is given by $nm_{\text{min}}$. According to kolesar2018minimum Proposition 4, the MCD test with asymptotic size $\varrho_n$ would reject the null of ${\boldsymbol{\alpha}}_{\hat{\mathcal{V}}}= \boldsymbol{0}$, when

where $\Phi(\cdot)$ denotes the CDF of the standard normal distribution and $\chi^2_{{{p_{\widehat{\mathcal{V}}^{c}}}} -1}(\varrho_n)$ is $1-\varrho_n$ quantile of $\chi^2_{{{p_{\widehat{\mathcal{V}}^{c}}}} -1}$ distribution. This property holds regardless of whether ${p_{\mathcal{V}^{c*}}}$ is fixed or grows with $n$. For the sake of the model selection consistency, the size of the MCD test needs to be $o(1)$. Following Belloni2012sparse,windmeijer2019use,windmeijer2019confidence, we adopt a scaled rate $\varrho_n = 0.5/\log n$ that works well in simulation studies. Thus, with $\widehat{\boldsymbol{\alpha}}^{(0)}$ in (ref) and a sequence of $\boldsymbol{\lambda}^{\text{seq}} =\boldsymbol{C} \sqrt{\log p/n}$ where $\boldsymbol{C} = 0.1\times(1,\ldots,20)^\top$, we propose to use the MCD test to select the proper $\lambda\in \boldsymbol{\lambda}^{\text{seq}}$ that is not rejected by (ref) with size $\varrho_n$ and largest ${{p_{\widehat{\mathcal{V}}}}} $.

To sum up, we provide Algorithm (ref) to demonstrate the detailed implementation.

Numerical Simulations

In this section, we conduct numerical studies to evaluate the finite sample performance of the proposed WIT estimator. In the design of the simulation experiments, we consider scenarios corresponding to different empirically relevant problems.

We consider the same model in Section (ref), $$ \boldsymbol{Y} = \boldsymbol{D}\beta^*+\boldsymbol{Z}\boldsymbol{\alpha}^*+\boldsymbol{\epsilon}, \quad \boldsymbol{D} = \boldsymbol{Z}\boldsymbol{\gamma}^*+\boldsymbol{\eta}. $$ Throughout all settings, we fix true treatment effect $\beta^* = 1$. $\boldsymbol{Z}$ is the $n\times p$ potential IV matrix and $\boldsymbol{Z_{i.}}\overset{i.i.d.}{\sim}\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}^{\boldsymbol{Z}})$, where $\boldsymbol{\Sigma}^{\boldsymbol{Z}}_{jj} = 0.3$ and $\boldsymbol{\Sigma}^{\boldsymbol{Z}}_{jk} = 0.3|j-k|^{0.8}$ for $i = 1,\ldots,n$ and $k,j = 1,\ldots,p$. Denote $\boldsymbol{\epsilon}= (\epsilon_1,\epsilon_2,\ldots,\epsilon_n)^\top$ and $\boldsymbol{\eta} = (\eta_1,\eta_2,\ldots,\eta_n)^\top$ and generate $(\epsilon_i,\eta_i)^\top \overset{\operatorname{i.i.d.}}{\sim}\mathcal{N}\left(\boldsymbol{0},\left(

\right)\right)$. We let $\sigma_\epsilon = 0.5$ and $\operatorname{corr}(\epsilon_i,\eta_i) = 0.6$ in all settings but vary $\sigma_\eta^2$ to get different concentration parameters concerning strong or weak IVs cases.

We compare the WIT estimator with other popular estimators in the literature. Specifically, sisVIVE is computed by R package sisVIVE; Post-ALasso windmeijer2019use, TSHT and CIIV are implemented using codes on Github guo2018confidence,windmeijer2019confidence. TSLS, LIML, oracle-TSLS and oracle-LIML (the truly valid set $\mathcal{V}^*$ is known a priori) are also included. Regarding our proposed WIT estimator, the MCD tuning strategy is implemented to determine $\lambda$, and we fix $\rho=2$. In the I-LAMM algorithm, we take $\delta_c = 10^{-3}$ and $\delta_t = 10^{-5}$ as the tolerance levels. We report results based on 500 simulations.

We measure the performance of all estimators in terms of median absolute deviation (MAD), standard deviation (STD), and coverage probability (CP) based on $95\%$ confidence intervals. Moreover, we provide measurements on the estimation of $\boldsymbol{\alpha}^*$ and IV model selection. Specifically, we measure the performance of invalid IVs selection by false positive rate (FPR) and false negative rate (FNR). To be concrete, denote the number of incorrect selections of valid and invalid IVs as FP and FN, respectively, and the number of correct selections of valid and invalid as TP and TN, respectively. Thus, FPR $=\text{FP} /(\text{FP}$ $+\text{TN})$ and $\text{FNR}=\text{FN} /(\text{FN}+\text{TP})$.

Case 1: Low dimension

We first consider the low dimension scenario:

In the above two cases, we maintain $\sigma_\eta = 1$ and vary the sample size $n = 200$ to $500$. Case 1(I) refers to all strong IVs case, but SAIS (ref) condition holds. Case 1(II) refers to Example (ref) with mixed strong and weak IVs.

Table (ref) presents the detailed results. In Case 1(I), high FPR or FNP indicates that sisVIVE and Post-Alasso mistarget $\boldsymbol{\alpha}^*$ because of lack of majority rule and SAIS holds. Their performances do not improve much with sample size $n$. Due to finite strong IVs, WIT performs similarly to TSHT, CIIV, and oracle-LIML. In low dimension settings, the oracle-TSLS is very close to oracle-LIML. Case 1(II) shows how weak IVs break the strong IVs-based plurality rule. As shown in Example (ref), TSHT and CIIV worsen in terms of all measures though the sample size increase from $n=200$ to $500$. Post-Alasso also fails due to the failure of the majority rule. As the analysis in Example (ref), such a mixture of weak IVs does not impede penalized methods. The WIT estimator outperforms the other estimators when $n=200$ and approaches oracle-LIML when $n=500$. Interestingly, the comparably low FPR and MAD imply sisVIVE correctly target true $\boldsymbol{\alpha}^*$ since the additional requirement of matching objective function (ref) happens to hold in this example. However, its FNR and MAD are worse than WIT due to the conservative cross-validation tuning strategy and the non-ignorable bias of Lasso, respectively.

Further, for closer comparison, we present a replication of the simulation design considered in windmeijer2019confidence and its weak IVs variant:

We now vary sample size $n = 500$ to $1000$ and fix $\sigma_\eta = 1$ to strictly follow their design. Between them, Case 1(III) corresponds to the exact setting, while Case 1(IV) scales down the magnitude of $\boldsymbol{\gamma}^*$ to introduce small coefficients in the first-stage.

Table (ref) shows the results. In Case 1(III), CIIV outperforms TSHT because CIIV can utilize available information better windmeijer2019confidence. The WIT estimator performs similar to CIIV and approaches oracle-LIML. sisVIVE and Post-Alasso fail again due to a lack of majority rule. In Case 1(IV), scaling down the first-stage coefficients causes some problems for CIIV and TSHT, since the first-stage selection thresholding $\sigma_\eta\sqrt{2.01\log p/n} = 0.111<0.15$, which might break the plurality rule numerically. TSHT and CIIV perform poorly when $n = 500$ and improve when $n =1000$ when the issue of violating the plurality rule is mitigated. Among penalized methods, sisVIVE and Post-Alasso mistarget and perform like TSLS because an additional requirement for sisVIVE (ref) and majority rule fail simultaneously. Distinguished from them, the WIT estimator outperforms with acceptable MAD when $n=500$. The FPR and FNR improve when the sample size increases. Fig. (ref) presents all replications in Case 1(IV) when $n = 1000$. It shows that WIT is nearly oracle besides the mild number of incorrect selections. By contrast, CIIV and TSHT fail in selections more frequently.

Case 2: High dimension (many IVs)

To assess performance in many IVs, we consider the following examples:\\ { Case 2(I) ($p = 0.5 n$): \\$\boldsymbol{\gamma}^* = (\boldsymbol{2\sqrt{\log p/n^{1-\delta}}})_{p}^{\top}$; $\boldsymbol{\alpha}^* = (\boldsymbol{0}_{0.6p},\boldsymbol{2\sqrt{\log {p}/n^{1-\delta}}}_{0.4p})^{\top}$.\\ Case 2(II) ($p = 0.6n$):\\ $\boldsymbol{\gamma}^* = (\boldsymbol{2\sqrt{\log{p}/n^{1-\delta}}})_{p}^{\top}$; $\boldsymbol{\alpha}^* = (\boldsymbol{0}_{0.5p},\boldsymbol{-2\sqrt{\log{p}/n^{1-\delta}}}_{0.2p},\boldsymbol{4\sqrt{\log{p}/n^{1-\delta}}}_{0.3p}$).\\ Among Case 2(I) and (II), $\delta = 0.01$ is taken.} The numbers of valid and invalid IVs are growing with the sample size. To verify the theoretical result, we maintain the ratio of concentration parameter to sample size $n$ at a low constant level, i.e.\ $\mu_n/n = 0.5$, by adjusting $\sigma^2_\eta$. { We vary the sample size $n$ from $500$ to $1000$, and let the first-stage coefficients in $\boldsymbol{\gamma}$ and elements in $\boldsymbol{\alpha}_{\mathcal{V}^c}$ go to $0$ simultaneously, which aligns with the case in Proposition (ref).} Due to the computational burden in CIIV for many IVs case, we omit its results.

{ Table (ref) provides the detailed estimation results. Case 2(I) satisfies the majority rule, and only two groups are present. All weak IVs narrow available choices in TSHT. When $n =1000$, a low FPR but high FNR indicates that TSHT only selected a limited number of valid IVs. Both low FPR and FNR imply sisVIVE targets the correct solutions. Meanwhile, a high MAD shows that the TSLS-based method is biased in many IVs. The majority rule holds to ensure Post-Alasso is consistent. However, when $n=500$, Post-Alasso is severely biased with relatively high FPR. It may result from the sensitivity problem of the initial estimator of Post-Alasso with weak IVs. The WIT estimator behaves similarly to oracle-LIML and achieves the best performance in every measure, even when $n=500$. On the contrary, oracle-TSLS has a larger bias than LIML in many IVs cases, and the coverage rate is poor.

Case 2(II) has more invalid IVs with the sparsest rule. In terms of valid IV selection, TSHT has subpar performance in this case since the strong IVs-based plurality rule is unlikely to hold. sisVIVE correctly identifies solutions by chance in this example. However, its estimate suffers from the bias of Lasso and TSLS-based bias in many IVs. Without the majority rule, Post-Alasso targets the wrong solution as the initial median estimator is inconsistent. On the other hand, WIT works when the majority rule does not hold and circumvents the bias from TSLS. Thus, WIT outperforms the other methods in many IVs cases and approaches oracle-LIML when the sample size increases. The numerical results demonstrate the capacity that WIT can work in many IV cases under the conditions that $\alpha_j^*$ and $\gamma_j^*$ share the same rate going to zero. These simulation studies strongly support the conclusion in Proposition (ref)}.

The Effect of BMI on Diastolic Blood Pressure

This section illustrates the usefulness of the proposed WIT estimator for the method of Mendelian Randomization. We implement the WIT estimator to obtain an estimate of the causal effect of BMI on diastolic blood pressure and compare it to the other estimators, OLS, sisVIVE, Post-Alasso, TSHT and CIIV. This comparative analysis is designed to exemplify the efficiency and robustness of the WIT estimator in various research scenarios.

{

Data Description

In Mendelian Randomization genetic markers, called single nucleotide polymorphisms (SNPs), function as IVs for the identification and estimation of causal effects of modifiable phenotypes on outcomes von2016genetic. This research design utilizes the random distribution of alleles at conception locke2015genetic.

Nonetheless, SNPs can be invalid IVs. This is primarily due to pleiotropy, which refers to the potential of genetic variants to be associated with multiple phenotypes, but can also be due to influences such as linkage disequilibrium and population stratification. Unfortunately, these violations are typically unidentified prior to their selection. Furthermore, the correlation between SNPs and treatments is often weak, and so we are dealing with the problem of potentially many weak instruments.

In our study, we analyzed data from 105,276 individuals from the UK Biobank to investigate the effect of body mass index (BMI) on diastolic blood pressure (DBP). As suggested by windmeijer2019use,windmeijer2019confidence, we used 96 SNPs as potential IVs for BMI. To account for skewness, we applied log-transformations to both BMI and DBP. We consider the same model specification detailed in Section 8 of windmeijer2019confidence. The model also included age, age squared, and sex as explanatory variables, along with 15 principal components of the genetic relatedness matrix.

Result and Analysis

Table (ref) provides the estimation results of the effect of $\log$(BMI) on $\log$(DBP).

The OLS estimate of 0.206 is potentially severely biased due to endogeneity issues, such as reverse causality or latent confounders.

Regarding the IV methods, we first denote TSHT$^*$ and CIIV$^*$ to indicate the estimation results without first-stage thresholding. For TSLS and TSHT$^*$, both without first-stage thresholding, they yield identical estimates of 0.087. However, the near-zero Sargan test p-value strongly rejects the model, implying TSHT's inability to detect potentially invalid IVs among weak ones. Despite the slight improvement of TSHT in discerning invalid IVs through first-stage thresholding, it still leaves some possibly invalid instruments in the model, as indicated by the slightly increased but still very small Sargan test p-values.

Compared to TSHT, CIIV identified thirteen strong and invalid IVs, leading to an estimate 0f 0.174 and a Sargan p-value of 0.011. CIIV$^*$ detected the same (as CIIV) strong invalid IVs with a similar Sargan p-value. CIIV$^*$ adjusted the estimate to 0.140.

The sisVIVE approach yielded an estimate of 0.111 which is lower than that obtained with the CIIV method. First, we observe that sisVIVE picks 20 invalid IVs (the highest count among comparison methods). Second, a Sargan p-value of 0.064 indicates that it overly penalizes invalid IVs while missing some true targets. In contrast, the Post-Alasso estimate is equal to 0.163, with the minimum number of 11 identified invalid IVs, and a Sargan p-value of 0.013.

The WIT estimator produced an estimate of 0.123, markedly lower than that of CIIV (and a little lower than CIIV$^*$), accompanied by the highest Sargan test p-value of 0.140. The invalid IVs detected by WIT encompassed all the relevant and invalid IVs identified by CIIV without first-stage thresholding. Moreover, compared to CIIV$^*$, WIT effectively identified two additional individually weak and invalid IVs. It thus significantly improved the Sargan p-value and estimation while capturing all valid yet weak IVs' information. In hindsight, sisVIVE penalized too many weak and valid IVs, leading to a loss in efficiency. Notably, the minimum validity violation $|\hat{\boldsymbol{\alpha}}^{MCP}|$ is 0.0014, aligned with the average magnitude of first-stage coefficients. In summary, the WIT method minimizes the risk of including invalid IVs while fully utilizing all valid IVs.

To summarize, the WIT estimator serves as a powerful tool for estimating treatment effects in biomedical research using SNPs as potential IVs. Its robustness to invalid and weak IVs makes it highly suitable for Mendelian Randomization applications, where there are potentially many weak IVs with uncertain validity.

}

Conclusion

We extended the study of IV models with unknown invalid IVs to allow for many weak IVs. We provided a complete framework to investigate the identification issue of such models. Sticking to the sparsest rule, we proposed the surrogate sparsest penalty that fits the identification condition. We proposed a novel WIT estimator that addresses the issues that can lead to poor performance of sisVIVE and Post-Alasso, and can outperform the plurality rule-based TSHT and CIIV. Simulations and real data analysis support the theoretical findings and the advantages of the proposed method over existing approaches.

Acknowledgements