# Fighting Noise with Noise: Causal Inference with Many Candidate Instruments

Xinyi Zhang, Linbo Wang, Stanislav Volgushev and Dehan Kong
University of Toronto, Department of Statistical Sciences

###### Abstract

Instrumental variable methods provide useful tools for inferring causal
effects in the presence of unmeasured confounding. To apply these methods with
large-scale data sets, a major challenge is to find valid instruments from a
possibly large candidate set. In practice, most of the candidate instruments
are often not relevant for studying a particular exposure of interest.
Moreover, not all relevant candidate instruments are valid as they may
directly influence the outcome of interest. In this article, we propose a
data-driven method for causal inference with many candidate instruments that
addresses these two challenges simultaneously. A key component of our proposal
is a novel resampling method, which constructs pseudo variables to remove
irrelevant candidate instruments having spurious correlations with the
exposure. Synthetic data analyses show that the proposed method performs
favourably compared to existing methods. We apply our method to a Mendelian
randomization study estimating the effect of obesity on health-related quality
of life.

Keywords: A/B tests; Mendelian randomization; Selection bias; Spurious
Correlation

##  1 Introduction

The instrumental variable (IV) model is a workhorse in causal inference using
observational data. It is especially useful when standard adjustment methods
are biased due to unmeasured confounding. The key idea is to find an exogenous
variable to extract random variation in the exposure and use this random
variation to estimate the causal effect. This exogenous variable is known as
the instrumental variable or instrument if it is related to the exposure, but
does not affect the outcome except through the exposure.

In the past, selection of appropriate exogenous variables that can serve as
valid instruments was typically driven by expert knowledge (e.g. Angrist &
Keueger, 1991). In recent years, applied researchers have increasingly sought
to employ large-scale observational data sets to identify causal
relationships. As the size and complexity of data sets grow, it becomes
increasingly difficult to build a knowledge-based instrumental variable model
for modern data applications. A prominent example in genetic epidemiology is
Mendelian randomization (MR), which uses genetic variants, such as single-
nucleotide polymorphisms (SNPs), as instruments to assess the causal
relationship between a risk factor and an outcome (e.g. Burgess et al., 2015).
The candidate instrument set includes millions of SNPs in the human genome.
This challenge also arises in the literature using multi-condition experiments
as instruments (e.g. Eckles et al., 2016; Goldman & Rao, 2016; Peysakhovich &
Eckles, 2018). For example, in tech companies, it is common to have routine
product experimentation such as A/B tests. These large collections of
randomized experiments provide natural candidates for finding instrumental
variables.

Modern large-scale complex systems call for data-driven methods that find
valid instruments from a high-dimensional candidate set. Two main challenges
arise from such a pursuit. First, often most of the candidate instruments are
not relevant for studying a particular exposure of interest. To deal with this
problem, a common approach is to employ a screening procedure (e.g. Fan & Lv,
2008) to select variables in the candidate set that are associated with the
exposure. For example, in genetics, it is common to use genome-wide
association studies (GWAS) or variants thereof to identify genetic variants
that are relevant to a risk factor. Recent studies have suggested choosing a
less conservative selection threshold than the traditional $5\times 10^{-8}$
for MR applications (e.g. Zhao et al., 2020), in order to include more SNPs
that are associated with the risk factors and hence improve the efficiency of
the resulting causal effect estimate; see also Barber & Candes (2019) for a
similar recommendation in a related context. This is especially important for
datasets with small to moderate sample sizes, for which standard GWAS has
limited power so that only very few or even no genetic variants may be
selected; see our data application in Section 6 for an example. However, a
less conservative selection threshold will result in more false positives,
i.e. SNPs having spurious correlations with the risk factor of interest. Using
these SNPs as instruments will undermine the validity of downstream MR
procedures. A similar problem also arises in multi-condition experiments as
most of the experiments have very small effects on the exposure (Peysakhovich
& Eckles, 2018).

Moreover, some of the candidate instruments found to be associated with the
exposure might not be valid as they may directly influence the outcome of
interest. For example, in MR studies, this problem arises due to pleiotropy, a
phenomenon that one genetic variant influences multiple traits and possibly
through independent pathways. In response to this concern, recently there has
been a surging interest in methods that make valid causal inference in the
presence of invalid instruments. Kolesar et al. (2015) and Bowden et al.
(2015) provided inferential methods for treatment effects in the presence of
invalid instruments. Their methods assume that the direct effects of these
invalid instruments are uncorrelated with the associations between the
instruments and the exposure. Cheng & Liao (2015) studied the selection of
valid instruments by assuming that one knows a collection of valid IVs a
priori. Another popular assumption, known as the “majority rule” (Bowden et
al., 2016; Kang et al., 2016; Windmeijer et al., 2019; Hartford et al., 2021),
assumes that at least 50% of the candidate instruments are valid. Relaxing the
majority rule, Guo et al. (2018) and Windmeijer et al. (2021) assume that the
largest group of Wald ratios with equal values corresponds to valid
instruments, referred to as the “plurality rule.” They also developed two-
stage thresholding procedures to identify valid IVs and consistently estimate
causal effects. In the MR context, Morrison et al. (2020) proposed a method
based on genome-wide summary statistics that accounts for both uncorrelated
and correlated pleiotropic effects, while Zhao et al. (2020) proposed a
consistent causal effect estimator by adjusting the profile likelihood of the
summary data under a two-sample MR design where the exposure and the outcome
are observed from two independent samples.

The developments in the present paper are motivated by the somewhat surprising
finding that a naive combination of existing methods for solving the two
challenges discussed earlier leads to biased causal inference. More precisely,
we observe that estimates from this naive combination are close to the
ordinary least squares (OLS) estimate that does not account for any unmeasured
confounding.

Our contributions in the paper are two-fold. We first provide a theoretical
characterization of the observation described above. In particular, we show
that estimates from exogenous variables having spurious correlations with the
exposure are concentrated in a region that can be separated from the true
causal effect. Remarkably in our theoretical analysis, we explicitly take into
account the randomness in the first marginal screening step, while most of the
existing literature treats this as pre-processing and does not account for its
randomness.

We then introduce a novel strategy that addresses this challenge. A key
component of our proposal is a resampling method that constructs pseudo
variables to identify and remove spurious variables. The idea of pseudo
variables has been investigated before in other contexts, most notably in
forward selection (Wu et al., 2007). It is also similar in spirit to knockoffs
(Barber & Candes, 2015) although our construction of pseudo variables is much
simpler.

To facilitate inference, we also consider a variation of our proposal using
sample splitting, in which one portion of the data is used to detect spurious
variables, and the remaining data are used to identify the valid instruments
and estimate the causal effects.

The rest of this article is organized as follows. In Section 2, we provide
background information of IV analysis and our models. We illustrate the
challenges of identifying valid instruments from a large candidate set using a
simulated example and provide theoretical justifications for our findings in
Section 3. Section 4 details the proposed procedures. In Section 5, we
evaluate the performance of our procedures via simulations. In Section 6, we
apply the proposed methods to the Wisconsin Longitudinal Study data and
examine the causal effect of obesity on Health-Related Quality of Life (HRQL).
We conclude with discussions in Section 7.

##  2 Background

The instrumental variable (IV) approach is a popular method for inferring
causality from observational studies. It introduces exogenous variables called
instruments to control for unmeasured confounding $U$. Conditioning on
observed covariates, a variable $Z$ is called an IV for estimating the causal
effect of an exposure $D$ on an outcome $Y$ if it satisfies the three core
conditions (e.g. Didelez & Sheehan, 2007; Wang & Tchetgen Tchetgen, 2018)
below:

  * I1

Relevance: $Z$ is associated with the exposure;

  * I2

Exclusion restriction: $Z$ has no direct effect on the outcome;

  * I3

Independence: $Z$ is independent of unmeasured variables that affect the
exposure and the outcome.

Figure 1 gives a causal graph representation of the IV model.

$U$$D$$Y$$Z$$\alpha_{D}^{*}$$\alpha_{Y}^{*}$$\beta^{*}$$\alpha_{ZD}^{*}$$\bigtimes$$\alpha_{ZY}^{*}\qquad$$\qquad\bigtimes$
Figure 1: Graphical illustration of an IV model: $Z$ denotes a valid
instrument, $D,Y$ represent the exposure and outcome, respectively and $U$ is
the unobserved confounding. Parameters on the arrows represent effects on the
linear scale. Red dashed lines with $\times$ indicate absence of arrows.

Suppose we have $n$ independent samples $(Y_{i},D_{i},\bm{Z}_{i},\bm{X}_{i})$
from the distribution of $(Y,D,\bm{Z},\bm{X})$, where $Y$ is a continuous
outcome of interest, $D$ is a continuous exposure,
$\bm{Z}=(Z_{1},\ldots,Z_{p})^{\mathrm{\scriptstyle T}}\in\mathbb{R}^{p}$
denotes candidate instrumental variables, and
$\bm{X}=(X_{1},\ldots,X_{q})^{\mathrm{\scriptstyle T}}\in\mathbb{R}^{q}$
represents baseline covariates. We assume that $D$ and $Y$ are generated from
the following linear structural equation models:

$\displaystyle D$ $\displaystyle=\bm{Z}^{\mathrm{\scriptstyle
T}}\bm{\gamma}^{*}+\bm{X}^{\mathrm{\scriptstyle
T}}\bm{\psi}^{*}+\bm{U}^{\mathrm{\scriptstyle
T}}\bm{\alpha}_{D}^{*}+\epsilon_{D};$ (2.1) $\displaystyle Y$
$\displaystyle=\bm{Z}^{\mathrm{\scriptstyle
T}}\bm{\pi}^{*}+\beta^{*}D+\bm{X}^{\mathrm{\scriptstyle
T}}\bm{\phi}^{*}+\bm{U}^{\mathrm{\scriptstyle
T}}\bm{\alpha}_{Y}^{*}+\epsilon_{Y}.$ (2.2)

The parameter $\beta^{*}$ measures the causal effect of the exposure $D$ on
the outcome $Y$, which is the quantity of primary interest. Parameter
$\bm{\gamma}^{*}\in\mathbb{R}^{p}$ characterizes the partial correlations
between the instruments and the exposure. Parameter
$\bm{\pi}^{*}\in\mathbb{R}^{p}$ measures the degree of violation of the
exclusion restriction assumption I2. Similarly, parameters
$\bm{\psi}^{*},\bm{\phi}^{*}\in\mathbb{R}^{q}$ refer to the effects of
observed covariates on the exposure and the outcome respectively. Parameters
$\bm{\alpha}_{D}^{*}\in\mathbb{R}^{g}$ and
$\bm{\alpha}_{Y}^{*}\in\mathbb{R}^{g}$ represent the direct effects of the
multivariate unobserved confounders $\bm{U}\in\mathbb{R}^{g}$ on $D$ and $Y$,
respectively. The dimension $g$ is fixed. Unmeasured confounding $\bm{U}$ and
random errors $\epsilon_{D},\epsilon_{Y}$ are independently distributed with
mean zero and covariance matrix/variance
$\Sigma_{U},\sigma_{D}^{2},\sigma_{Y}^{2}$, respectively. We further assume
that $Z_{j}$ is independent of $\epsilon_{D},\epsilon_{Y}$ and $\bm{U}$ for
$j=1,\ldots,p$. This assumption is generally considered plausible in MR
studies after controlling for population stratification (Smith, 2007; Zhao et
al., 2020) and holds by design in multi-condition experiments. An instrument
$Z_{j}$ is called relevant if $\gamma_{j}^{*}\neq 0$, and is called irrelevant
otherwise. For relevant instruments, if $\pi_{j}^{*}=0$, then $Z_{j}$ is a
valid instrument; otherwise, it is called an invalid instrument.

For the rest of the paper, we suppress dependence on covariates $\bm{X}$ for
simplicity. With this, models (2.1) and (2.2) take the reduced forms

$\displaystyle D=\bm{Z}^{\mathrm{\scriptstyle
T}}\bm{\gamma}^{*}+\bm{U}^{\mathrm{\scriptstyle
T}}\bm{\alpha}_{D}^{*}+\epsilon_{D},$ (2.3) $\displaystyle
Y=\bm{Z}^{\mathrm{\scriptstyle T}}\bm{\Gamma}^{*}+\bm{U}^{\mathrm{\scriptstyle
T}}(\beta^{*}\bm{\alpha}_{D}^{*}+\bm{\alpha}_{Y}^{*})+e,$ (2.4)

where $\bm{\Gamma}^{*}=\bm{\pi}^{*}+\beta^{*}\bm{\gamma}^{*}$ and
$e=\beta^{*}\epsilon_{D}+\epsilon_{Y}$. Let
$\bm{\mathbb{Z}}\in\mathbb{R}^{n\times p},\bm{\mathrm{D}}\in\mathbb{R}^{n}$
and $\bm{\mathrm{Y}}\in\mathbb{R}^{n}$ denote the collection of observations
for $(\bm{Z},D,Y)$ from $n$ independent units, respectively. Here
$\bm{\mathbb{Z}}$ is also called the design matrix. Without loss of
generality, assume that $\bm{\mathbb{Z}},\bm{\mathrm{D}},\bm{\mathrm{Y}}$ are
all centered.

Given a valid instrument $Z_{j},$ the causal effect $\beta^{*}$ may be
identified as the Wald ratio (Lawlor et al., 2008). Alternatively, given a set
of valid instruments $\bm{Z}_{valid}$, a popular procedure to estimate the
causal effect under linear models is the two-stage least squares (2SLS)
estimator given by
$\widehat{\beta}_{\text{2SLS}}=(\bm{\mathrm{D}}^{\mathrm{\scriptstyle
T}}\bm{\mathrm{P}}\bm{\mathrm{D}})^{-1}(\bm{\mathrm{D}}^{\mathrm{\scriptstyle
T}}\bm{\mathrm{P}}\bm{\mathrm{Y}}),$ where
$\bm{\mathrm{P}}=\bm{\mathbb{Z}}_{valid}(\bm{\mathbb{Z}}_{valid}^{\mathrm{\scriptstyle
T}}\bm{\mathbb{Z}}_{valid})^{-1}\bm{\mathbb{Z}}_{valid}^{\mathrm{\scriptstyle
T}}$ is a projection matrix and $\bm{\mathbb{Z}}_{valid}$ is the design matrix
corresponding to the valid instruments.

In practice, prior knowledge of instrument validity is often not available.
Nevertheless, the causal effect can still be identified and estimated from
data under the majority rule that more than 50% of the candidate IVs are
valid. In this case, among the Wald ratios constructed using individual
candidate IVs, more than 50% are equal to the true causal effect $\beta^{*}$.
As a result, $\beta^{*}$ can be identified as the median of these Wald ratios
(e.g. Han, 2008; Bowden et al., 2016). The majority rule may be relaxed to the
so-called plurality rule, which assumes that the largest group of invalid
instruments with the same Wald ratio is smaller than the group of valid
instruments, whose Wald ratio equals the true causal effect. In this case, the
mode of the Wald ratios constructed using individual candidate IVs corresponds
to the true causal effect (e.g. Guo et al., 2018; Windmeijer et al., 2021).

##  3 Challenges for Causal Inference with Many Candidate Instruments

###  3.1 A naive combination of existing approaches with screening

In this section, we consider the causal effect estimation problem under models
(2.1) and (2.2), where the number of candidate instruments far exceeds the
sample size. For example, in the genetics application detailed in Section 6,
there are 3,683,868 candidate instruments measured on a data set of 3023
subjects. The ultra-high dimensional nature of the candidate instrument set
makes standard regularization methods like the lasso computationally
challenging. In this case, it is common to apply a marginal screening
procedure to quickly screen out variables that are unlikely to be relevant to
the exposure (e.g. Fan & Lv, 2008). For example, one may select the top $s$
candidate instruments by ranking the absolute marginal estimates
$|\sum_{i=1}^{n}D_{i}Z_{ij}|/(\sum_{i=1}^{n}Z_{ij}^{2})$ for $1\leq j\leq p$.
We then let $\widehat{\mathcal{S}}_{1}\subset\\{1,\ldots,p\\}$ be the index
set of the top $s$ candidate instruments, and
$\bm{Z}_{\widehat{\mathcal{S}}_{1}}$ be the corresponding candidate
instruments. We denote the design matrix of
$\bm{Z}_{\widehat{\mathcal{S}}_{1}}$ as
$\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}\in\mathbb{R}^{n\times s}$.

The marginal screening step may include candidate instruments that are
marginally correlated but conditionally uncorrelated with the exposure. To
further estimate the set of relevant instruments defined in a joint model, we
first evaluate joint relevance strength for candidates
$\bm{Z}_{\widehat{\mathcal{S}}_{1}}$ selected from marginal screening and
subsequently perform a hard thresholding step (Donoho & Johnstone, 1994). If
$\widehat{\mathcal{S}}_{1}$ includes all the relevant instruments, then models
(2.1) and (2.2) still hold with $\bm{Z}\in\mathbb{R}^{p}$ replaced by
$\bm{Z}_{\widehat{\mathcal{S}}_{1}}\in\mathbb{R}^{s}$, and parameters
$\bm{\gamma}^{*}$ and $\bm{\pi}^{*}$ replaced by
$\bm{\gamma}_{\widehat{\mathcal{S}}_{1}}^{*},\bm{\pi}_{\widehat{\mathcal{S}}_{1}}^{*}\in\mathbb{R}^{s}$,
respectively. The parameter $\bm{\Gamma}_{\widehat{\mathcal{S}}_{1}}^{*}$ in
the reduced model is defined similarly to $\bm{\Gamma}^{*}$ as
$\bm{\Gamma}_{\widehat{\mathcal{S}}_{1}}^{*}=\bm{\pi}_{\widehat{\mathcal{S}}_{1}}^{*}+\beta^{*}\bm{\gamma}_{\widehat{\mathcal{S}}_{1}}^{*}$.
To estimate the coefficients $\bm{\gamma}_{\widehat{\mathcal{S}}_{1}}^{*}$ and
$\bm{\Gamma}_{\widehat{\mathcal{S}}_{1}}^{*}$ in potentially high dimensional
settings, we proceed using the de-biased lasso estimator in van de Geer et al.
(2014), which is defined as follows,

$\displaystyle\widehat{\bm{\Gamma}}_{\widehat{\mathcal{S}}_{1}}$
$\displaystyle=\widetilde{\bm{\Gamma}}_{\widehat{\mathcal{S}}_{1}}+\frac{1}{n}\bm{M}\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{Y}}-\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}\widetilde{\bm{\Gamma}}_{\widehat{\mathcal{S}}_{1}}),$
$\displaystyle\widehat{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}}$
$\displaystyle=\widetilde{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}}+\frac{1}{n}\bm{M}\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{D}}-\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}\widetilde{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}}),$
(3.1)

with the lasso estimator
$\widetilde{\bm{\Gamma}}_{\widehat{\mathcal{S}}_{1}}=\mathop{\mathrm{argmin}}_{\bm{\Gamma}\in\mathbb{R}^{s}}\\{\|\bm{\mathrm{Y}}-\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}\bm{\Gamma}\|_{2}^{2}/n+2\lambda_{\Gamma}\|\bm{\Gamma}\|_{1}\\},$
and the lasso estimator $\widetilde{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}}$
is obtained similarly. Here the matrix $\bm{M}\in\mathbb{R}^{s\times s}$ is an
estimate of the precision matrix
$\bm{\Sigma}_{\widehat{\mathcal{S}}_{1}}^{-1}$; see Appendix G for more
details on the construction of $\bm{M}$. Following Donoho & Johnstone (1994),
the hard thresholding step proceeds by comparing
$|\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}|$ with a threshold
$\delta_{n}\times\text{SE}(\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l})$,
where $\delta_{n}=\sqrt{\omega\log(\max\\{n,s\\})}$ and the standard error of
the de-biased lasso estimate $\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}$
for $1\leq l\leq s$ is

$\text{SE}(\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l})=\sqrt{[\bm{M}\widehat{\bm{\Sigma}}_{\widehat{\mathcal{S}}_{1}}\bm{M}^{\mathrm{\scriptstyle
T}}/n]_{ll}\|\bm{\mathrm{D}}-\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}\widetilde{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}}\|^{2}/n}$
(3.2)

with $[]_{ll}$ representing the $l$th diagonal element of the matrix. Here
$\omega$ in the threshold is a tuning parameter that controls the number of
candidates estimated as relevant to the exposure. We use
$\widehat{\mathcal{S}}_{2}$ to denote the set of candidate instruments that
pass this joint thresholding step, that is,
$\widehat{\mathcal{S}}_{2}=\\{l:1\leq l\leq
s,|\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}|\geq\delta_{n}\times\text{SE}(\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l})\\}$.

Now suppose that the plurality rule holds for $\widehat{\mathcal{S}}_{2}$,
i.e. the number of invalid instruments with the same Wald ratio
$\Gamma_{\widehat{\mathcal{S}}_{1},l}^{*}/\gamma_{\widehat{\mathcal{S}}_{1},l}^{*}$
for $l\in\widehat{\mathcal{S}}_{2}$ is smaller than the number of valid
instruments in $\widehat{\mathcal{S}}_{2}$. To estimate valid instruments
under this plurality rule, one may apply the voting framework by Guo et al.
(2018). Specifically, their voting procedure uses pairwise comparisons to test
whether two candidates have the same Wald ratios by testing
$H_{0}^{(j,l)}:\Gamma_{\widehat{\mathcal{S}}_{1},j}^{*}-\gamma_{\widehat{\mathcal{S}}_{1},j}^{*}(\Gamma_{\widehat{\mathcal{S}}_{1},l}^{*}/\gamma_{\widehat{\mathcal{S}}_{1},l}^{*})=0$.
For each voter $l\in\widehat{\mathcal{S}}_{2}$, the collection
$\big{\\{}j\in\widehat{\mathcal{S}}_{2}:H_{0}^{(j,l)}\text{ not
rejected}\big{\\}}$ would be candidates that $l$ is in favor of. Once the
polls close, by counting the number of votes for each candidate, the
candidates that win the most votes are estimated as valid instruments. See Guo
et al. (2018) for more detailed discussions.

We summarize the aforementioned procedure in Algorithm 1, which simply
combines marginal screening and existing work that identifies valid
instruments under the plurality rule.

INPUT: Design matrix $\bm{\mathbb{Z}}\in\mathbb{R}^{n\times p}$, observed
exposure $\bm{\mathrm{D}}\in\mathbb{R}^{n}$, observed outcome
$\bm{\mathrm{Y}}\in\mathbb{R}^{n}$.

  * 1.

Marginal screening (Fan & Lv, 2008):

For $j=1,\ldots,p$, compute
$\widehat{\rho}_{j}=|\sum_{i=1}^{n}Z_{ij}D_{i}|/\sum_{i=1}^{n}Z_{ij}^{2}$.

Select the $s$ variables with the largest $\widehat{\rho}_{j}$ and denote them
as $\widehat{\mathcal{S}}_{1}$.

  * 2.

Joint thresholding (Donoho & Johnstone, 1994; van de Geer et al., 2014):

Fit a linear model for
$\bm{\mathrm{D}}\sim\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}$: let
$\widehat{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}}$ be the debiased lasso
estimate defined in (3.1); here $\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}$
is formed by columns of $\bm{\mathbb{Z}}$ indexed by
$\widehat{\mathcal{S}}_{1}$. Also obtain their standard error estimate
$\text{SE}(\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l})$ in (3.2).

Fit a linear model for
$\bm{\mathrm{Y}}\sim\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}$ and obtain
the debiased lasso estimate
$\widehat{\bm{\Gamma}}_{\widehat{\mathcal{S}}_{1}}$ in a similar fashion.

Estimate the relevant variables $\widehat{\mathcal{S}}_{2}=\big{\\{}1\leq
l\leq
s:|\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}|\geq\delta_{n}\times\text{SE}(\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l})\big{\\}},$
where $\delta_{n}=\sqrt{\omega\log(\max\\{n,s\\})}$ and $\omega$ is a tuning
parameter.

  * 3.

Voting (Guo et al., 2018):

For $j,l\in\widehat{\mathcal{S}}_{2}$, test
$H_{0}^{(j,l)}:\Gamma_{\widehat{\mathcal{S}}_{1},j}^{*}-\gamma_{\widehat{\mathcal{S}}_{1},j}^{*}\Big{(}\frac{\Gamma_{\widehat{\mathcal{S}}_{1},l}^{*}}{\gamma_{\widehat{\mathcal{S}}_{1},l}^{*}}\Big{)}=0.$

For $l\in\widehat{\mathcal{S}}_{2}$, find candidates that $l$ votes for:
$\widehat{\mathcal{V}}_{l}=\big{\\{}j\in\widehat{\mathcal{S}}_{2}:H_{0}^{(j,l)}\text{
not rejected}\big{\\}}.$

Count votes for $j\in\widehat{\mathcal{S}}_{2}$:
$\text{vc}_{j}=\sum_{l\in\widehat{\mathcal{S}}_{2}}\mathbbm{1}\big{(}j\in\widehat{\mathcal{V}}_{l}\big{)}$.

Find valid IVs:
$\widehat{\mathcal{S}}_{4}=\big{\\{}\mathop{\mathrm{argmax}}_{j\in\widehat{\mathcal{S}}_{2}}\text{vc}_{j}\big{\\}}.$

  * 4.

Causal effect estimation: Take
$\widehat{\mathcal{S}}_{4},\bm{\mathrm{D}},\bm{\mathrm{Y}}$ as input into 2SLS
to estimate causal effect and construct a confidence interval.

Algorithm 1 Naive algorithm: causal effect estimation with many candidate
instruments

###  3.2 A numerical example

In this section we analyze the naive approach from the previous section in a
simulated example. We generate 1000 data sets of size $n=500$ with $p=50,000$
candidate instruments. The true causal effect $\beta^{*}$ is set to be 2.
Candidate instruments $Z_{1},Z_{2}$ are invalid, $Z_{3},\ldots,Z_{9}$ are
valid, and $Z_{10},\ldots,Z_{50,000}$ are irrelevant. The valid instruments
are distributed according to $\bm{Z}_{3:9}\sim N(\bm{0},\Sigma_{v})$ with
$[\Sigma_{v}]_{jk}=0.25^{|j-k|}$. The other candidate instruments
$Z_{j},j=1,2,10,\ldots,50,000$, observed covariates $X_{1},X_{2}$, random
error $\epsilon_{Y}$, and a univariate unobserved confounding $U$ are
generated from independent standard normal distributions. The treatment $D$
and the outcome $Y$ follow models (2.1) and (2.2), where
$(\alpha_{D}^{*},\alpha_{Y}^{*})=(4,-3)$, the relevance strength
$\bm{\gamma}_{1:9}^{*}=3$ and the degree of violation to the exclusion
restriction assumption is characterized by
$\bm{\pi}_{1:2}^{*}=(-3.5,3.5)^{\mathrm{\scriptstyle T}}$; the effects of the
observed covariates on the exposure and the outcome are given by
$\bm{\psi}^{*}=(1.5,2)^{\mathrm{\scriptstyle T}}$ and
$\bm{\phi}^{*}=(1.2,1.5)^{\mathrm{\scriptstyle T}}$, respectively. For
simplicity, we let $\epsilon_{D}=0$.

To reduce the dimension of candidate instruments and exclude variables that
are unlikely to be informative to the exposure, following Step 1 of Algorithm
1, we first marginally regress the exposure on each candidate instrument and
keep the top 500 with the largest marginal associations in absolute value. We
then apply Steps 2 – 4 of Algorithm 1 to this remaining set for identifying
valid IVs and estimating the causal effect. Empirical results suggest that
most of the valid instruments identified by Algorithm 1 (Naive) are in fact
irrelevant. On average across 1000 Monte Carlo runs, the “Naive” algorithm
identifies 15.27 variables as valid IVs, among which 15.25 are in fact
irrelevant. The resulting causal effect estimate is 1.78 on average and the
empirical coverage probability is 1%, compared to a nominal level of 95%.

To provide insights into this problem, we study the causal effect estimates
from $Z_{j}^{\prime}s$ with $j\in\widehat{\mathcal{S}}_{2}.$ In Figure 2(a),
we plot the histograms of causal effect estimates aggregated over the 1000
Monte Carlo runs. We use different colors to distinguish the estimates
produced by the valid, invalid IVs and irrelevant variables. One can see that
the plurality rule is violated empirically even after joint thresholding: the
largest number of causal effect estimates with similar values corresponds to
irrelevant variables. The latter variables shall be referred to as “spurious
instruments.” The plurality rule is violated here because there are much more
irrelevant candidate instruments to begin with, and the sample size is small
relative to the number of candidate instruments. Somewhat surprisingly, all
spurious instruments lead to similar causal effect estimates. Moreover, these
effect estimates lie in a region that is separated from those estimates
regarding the valid IVs. As a result, Algorithm 1 that assumes the plurality
rule misidentifies the spurious instruments as valid IVs, as shown in Figure
2(b). A theoretical explanation of those findings is given in the following
subsection.

(a) $\widehat{\mathcal{S}}_{2}$ (b) $\widehat{\mathcal{S}}_{4}$

Figure 2: Histograms of causal effect estimates aggregated over the 1000 Monte
Carlo runs: (a) plots the estimates by candidate instruments that pass the
joint thresholding step ($\widehat{\mathcal{S}}_{2}$) of Algorithm 1, and (b)
plots the estimates by candidate instruments identified as valid instruments
($\widehat{\mathcal{S}}_{4}$) by Algorithm 1. The blue pile represents causal
effect estimates by valid IVs; the green piles represent those by invalid IVs;
the orange pile represents those by irrelevant variables. The red dashed line
is the true causal effect $\beta^{*}=2$.

###  3.3 Results for causal estimates from variables passing joint
thresholding

In this section, we theoretically characterize the behaviour of causal effect
estimates from spurious, valid, and invalid instruments that pass joint
thresholding. Before introducing our concentration bounds, we define the sets
of relevant and irrelevant variables in $\widehat{\mathcal{S}}_{1}$, the set
of variables selected by marginal screening

$\displaystyle\check{\mathcal{R}}=\\{l:1\leq l\leq
s,\gamma_{\widehat{\mathcal{S}}_{1},l}^{*}\neq
0\\},\quad\check{\mathcal{I}}=\\{l:1\leq l\leq
s,\gamma_{\widehat{\mathcal{S}}_{1},l}^{*}=0\\}.$

The assumptions imposed to establish the theoretical results are presented in
the following:

  * A1. Unmeasured confounders $\bm{U}_{i}$’s are $i.i.d.\text{ }N(\bm{0},\Sigma_{U})$ where $\Sigma_{U}\in\mathbb{R}^{g\times g}$ and $g$ is fixed. Random errors $\epsilon_{iD}$’s are $i.i.d.\text{ }N(0,\sigma_{D}^{2})$ and $\epsilon_{iY}$’s are $iid$ $N(0,\sigma_{Y}^{2})$. Further $\epsilon_{iD},\epsilon_{iY}$ and $\bm{U}_{i}$ are mutually independent for $i=1,\ldots,n$ and are all independent of $Z_{ij}$ for $1\leq i\leq n$ and $1\leq j\leq p$.

  * A2. Let $\bm{\Sigma}\in\mathbb{R}^{p\times p}$ be the covariance matrix of $\bm{Z}\in\mathbb{R}^{p}$ . There exist constants $C_{L},C_{R},B>0$ such that for all $p$, $\lambda_{\min}(\bm{\Sigma})\geq C_{L}>0$, $\lambda_{\max}(\bm{\Sigma})\leq C_{R}<\infty$, and $\max_{1\leq j\leq p}\bm{\Sigma}_{jj}\leq B$. Further assume $\bm{Z}_{i}$ are $i.i.d$ sub-gaussian with variance proxy $\sigma^{2}$ and zero mean, i.e., for any $\lambda\in\mathbb{R}$ and unit vector $\bm{u}\in\mathcal{S}^{p-1}$, $E\left\\{e^{\lambda\bm{u}^{\mathrm{\scriptstyle T}}(\bm{\Sigma}^{-1/2}\bm{Z}_{1})}\right\\}\leq e^{\frac{\sigma^{2}\lambda^{2}}{2}}$ . Denote by $\kappa=\|\bm{\Sigma}^{-1/2}\bm{Z}_{1}\|_{\psi_{2}}$ the sub-Gaussian norm of $\bm{\Sigma}^{-1/2}\bm{Z}_{i}$ for $i=1,\ldots,n$, where $\kappa$ is a finite positive constant.

  * A3. Variables irrelevant the exposure are also irrelevant to the outcome, i.e. $\\{k:\gamma_{k}^{*}=0\\}\subseteq\\{k:\pi_{k}^{*}=0\\}$.

  * A4. The number of variables relevant to the exposure $D$, denoted by $s_{1}^{*}$, is fixed. The number of invalid instruments that violate the exclusion restriction assumption (I2) is also fixed and denoted by $s_{2}^{*}$.

  * A5. There exist positive $c_{0},C_{0},\breve{c}_{0},\breve{C}_{0}$ such that $c_{0}\leq|\gamma_{k}^{*}|\leq C_{0}$ for all relevant variables $k$ and for all invalid instrumenst $k$, the degree of violation satisfies $\breve{c}_{0}\leq|\pi_{k}^{*}|\leq\breve{C}_{0}$. There exists $c>0$ such that $\min_{k:\gamma_{k}^{*}\neq 0}|\mathrm{cov}(D/\gamma_{k}^{*},Z_{k})|\geq c$.

  * A6. Let $\eta_{max}=\max(1,B\sup_{\mathcal{S}:|\mathcal{S}|=s}\max_{1\leq j\leq s}\|[\bm{\Sigma}_{\mathcal{S}}^{-1}]_{-j,j}\|_{1})$, where $\bm{\Sigma}_{\mathcal{S}}^{-1}$ denotes the inverse of $\bm{\Sigma}_{\mathcal{S}}$ formed by rows and columns of $\bm{\Sigma}$ indexed by $\mathcal{S}$ and $[]_{-j,j}$ represents the $j$th column without the $j$th component of the matrix. Assume $\eta_{\max}=o\big{(}(n/\log p)^{1/8}\big{)}$.

  * A7. Let $\lambda_{\gamma}=C_{\gamma}\sqrt{(\log p)/n}$ with $C_{\gamma}\geq 5.7\sqrt{B(\bm{\alpha}_{D}^{*\mathrm{\scriptstyle T}}\Sigma_{U}\bm{\alpha}_{D}^{*}+\sigma_{D}^{2})}$ and $\lambda_{\Gamma}=C_{\Gamma}\sqrt{(\log p)/n}$ with $C_{\Gamma}\geq 5.7\sqrt{B\\{(\beta^{*}\bm{\alpha}_{D}^{*}+\bm{\alpha}_{Y}^{*})^{\mathrm{\scriptstyle T}}\Sigma_{U}(\beta^{*}\bm{\alpha}_{D}^{*}+\bm{\alpha}_{Y}^{*})+\beta^{*2}\sigma_{D}^{2}+\sigma_{Y}^{2}\\}}$, where $\lambda_{\gamma},\lambda_{\Gamma}$ are tuning parameters in the lasso problems defined after (3.1). Moreover, the tuning parameters in the nodewise lasso problems (G.1) satisfy $\lambda_{j}\equiv 4a_{2}\eta_{\max}\sqrt{(\log p)/n},j=1,\ldots,s$ and a constant $a_{2}>4\sqrt{3}eB\kappa^{2}$.

Assumptions A1 and A2 are distribution assumptions for candidate instruments,
unmeasured confounding and random errors in the models (2.1) and (2.2). In A3
we for simplicity assume that candidates that are irrelevant to the exposure
are also irrelevant to the outcome. This can be relaxed at the cost of
additional technicalities in the arguments. Similarly, assumption A4 is made
for simplicity and can be relaxed at the cost of more technical arguments. The
relevance strength for valid and invalid instruments is provided in A5; A5
also rules out candidates that are marginally uncorrelated but jointly
correlated with the exposure. We can allow the relevance strength to decay at
rates $c_{0}\gtrsim\eta_{\max}\sqrt{(\log p)/n}$ at the cost of introducing
additional technicalities in the proofs. Assumption A6 gives conditions on
maximal row-wise $\ell_{1}$ norm for precision matrix
$\bm{\Sigma}_{\widehat{\mathcal{S}}_{1}}^{-1}$ uniformly over all
$\widehat{\mathcal{S}}_{1}$ of cardinality $s$. Restrictions on this qunatity
are also considered in van de Geer et al. (2014). Conditions for lasso tuning
parameters are collected in A7.

####  3.3.1 Concentration of estimates from spurious instruments

We now consider concentration of Wald ratios from spurious instruments that
pass joint thresholding.

######  Theorem 3.1 (Concentration results for spurious instruments).

Suppose assumptions A1–A7 hold. Further assume $\log p=O(n^{\tau_{1}})$ for
$0<\tau_{1}<1$, $n^{\tau_{2}}=O(s)$ for $0<\tau_{2}<1$ and $(\log
p)^{3/4}=o(n^{1/4}\sqrt{\log s})$. Let

$\widetilde{C}:=\frac{8\Big{\\{}\mathrm{var}(\bm{\alpha}_{Y}^{*\mathrm{\scriptstyle
T}}\bm{U})-\frac{\mathrm{cov}^{2}(\bm{\alpha}_{Y}^{*\mathrm{\scriptstyle
T}}\bm{U},\bm{\alpha}_{D}^{*\mathrm{\scriptstyle
T}}\bm{U})}{\mathrm{var}(\bm{\alpha}_{D}^{*\mathrm{\scriptstyle
T}}\bm{U}+\epsilon_{D})}+\sigma_{Y}^{2}\Big{\\}}^{1/2}}{\sqrt{\mathrm{var}(\bm{\alpha}_{D}^{*\mathrm{\scriptstyle
T}}\bm{U}+\epsilon_{D})}}$

There exists a constant $C_{1}>0$ independent of $\omega$ such that for $n,p$
sufficiently large

$P\bigg{(}\cup_{l\in\check{\mathcal{I}}\cap\widehat{\mathcal{S}}_{2}}\bigg{\\{}\bigg{|}\frac{\widehat{\Gamma}_{\widehat{\mathcal{S}}_{1},l}}{\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}}-\beta^{*}-C_{*}\bigg{|}>\widetilde{C}/\sqrt{\omega}\bigg{\\}}\bigg{)}\leq
2s^{-3}+e^{-C_{1}\log(n\wedge p)}.$

Theorem 3.1 formally verifies the empirical observation that estimates from
spurious instruments concentrate in a region around $\beta^{*}+C_{*}$. The
width of this region shrinks for higher threshold parameters $\omega$ and can
be made small if $\omega$ is chosen sufficiently large.

The centering $\beta^{*}+C_{*}$ coincides with the ordinary least squares
estimate, and also appears in the weak IV literature (Nelson & Startz, 1990b;
Bound et al., 1995; Stock et al., 2002) which suggest that as the instrumental
strength approaches zero, the causal effect estimate approaches
$\beta^{*}+C_{*}$. These results, however, assume that the number of weak
instruments is fixed or increases at a slower rate than the sample size.
Moreover, they place additional constraints on the convergence rate of the
relevance strength of weak IVs. In contrast, we allow the number of irrelevant
candidate instruments before screening to grow at an exponential rate and
explicitly take into account randomness from both – first stage screening and
hard thresholding based on the debiased lasso. This makes our theoretical
analysis substantially more involved.

####  3.3.2 Concentration of estimates by relevant instruments

Next, consider the behaviour of causal effect estimates from relevant
instruments that pass joint thresholding.

######  Theorem 3.2 (Concentration results for relevant (valid and invalid)
instruments).

Suppose assumptions A1, A2 and A4–A7 hold. Further assume $\log
p=O(n^{\tau_{1}})$ for $0<\tau_{1}<1$ and $n^{\tau_{2}}=O(s)$ for
$0<\tau_{2}<1$. Then the set $\widehat{\mathcal{S}}_{2}$ includes all relevant
variables with probability going to one. Moreover, there exist constants
$C_{1},C_{2}>0$ such that the following holds when $n,p$ are sufficiently
large

$\displaystyle\quad
P\bigg{(}\max_{l\in\check{\mathcal{R}}\cap\widehat{\mathcal{S}}_{2}}\bigg{|}\frac{\widehat{\Gamma}_{\widehat{\mathcal{S}}_{1},l}}{\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}}-\bigg{(}\beta^{*}+\frac{\pi_{\widehat{\mathcal{S}}_{1},l}^{*}}{\gamma_{\widehat{\mathcal{S}}_{1},l}^{*}}\bigg{)}\bigg{|}>C_{2}\sqrt{\frac{\log
p}{n}}\frac{\eta_{max}}{C_{L}}\bigg{)}\leq e^{-C_{1}\log(n\wedge p)}.$

Theorem 3.2 confirms that estimates from valid IV’s concentrate around the
true causal effect and causal effect estimates from invalid IVs which violate
the exclusion restriction assumption (I2) concentrate around
$\beta^{*}+\pi_{\widehat{\mathcal{S}}_{1},l}^{*}/\gamma_{\widehat{\mathcal{S}}_{1},l}^{*}$,
even after marginal screening and joint threholding via debiased lasso
estimates.

Combining the findings in Theorem 3.1 and Theorem 3.2 shows that the estimates
from valid instruments can be separated from those from spurious instruments
with probability going to one for $\omega$ satisfying

$\frac{1}{\sqrt{\omega}}<\frac{|\mathrm{cov}(\bm{\alpha}_{D}^{*\mathrm{\scriptstyle
T}}\bm{U},\bm{\alpha}_{Y}^{*\mathrm{\scriptstyle
T}}\bm{U})|}{\widetilde{C}\mathrm{var}(\bm{\alpha}_{D}^{*\mathrm{\scriptstyle
T}}\bm{U}+\epsilon_{D})}$ (3.3)

Indeed, from Theorem 3.2 we see that causal estimates from all valid
instruments will concentrate around the true effect $\beta^{*}$. In contrast,
estimates from spurious instruments concentrate within the interval
$[\beta^{*}+C_{*}-\widetilde{C}/\sqrt{\omega},\beta^{*}+C_{*}+\widetilde{C}/\sqrt{\omega}]$
as shown in Theorem 3.1. This interval does not contain $\beta^{*}$ provided
that $\widetilde{C}/\sqrt{\omega}<|C_{*}|$, which can always be achieved by
choosing $\omega$ sufficiently large.

Results in this section formalize the intuition for why a naive application of
the procedure to identify valid instruments under plurality after applying
marginal screening can fail – all estimates from spurious instruments
concentrate around the same value, leading to a violation of the empirical
plurality rule. Our results also suggest a way out: since valid and invalid
instruments still concentrate around the same values as without marginal
screening, and since all estimates from spurious instruments are similar and
separated from valid instruments, all that we need is to identify estimates
from spurious instruments. This motivates the methodology in the next section.

##  4 Fighting Noise with Noise

As established in the previous sections, a major challenge for causal
inference with many candidate instruments is the presence of spurious
instruments and distinguishing the valid instruments from the spurious ones.
To solve this problem, the key idea in our developments is to add independent
noises, known as pseudo variables, to the initial candidate instrument set.
These pseudo variables mimic the behaviour of the irrelevant variables in the
original set, as both of them are independent of the exposure $D$.

Specifically, since the vast majority of candidate instruments in most
applications are irrelevant, we propose to generate $p$ pseudo variables
$Z_{p+1},\ldots,Z_{2p}$ by randomly permuting the rows of the design matrix
$\bm{\mathbb{Z}}\in\mathbb{R}^{n\times p}$. The resulting matrix
$\bar{\bm{\mathbb{Z}}}\in\mathbb{R}^{n\times p}$ is independent of $D$ and has
the same correlation structure as the original candidate instruments $\bm{Z}$.
The expanded candidate instrument matrix
$\widetilde{\bm{\mathbb{Z}}}\in\mathbb{R}^{n\times 2p}$ is defined as the
concatenation of $\bm{\mathbb{Z}}$ and $\bar{\bm{\mathbb{Z}}}$; see step 0 of
Algorithm 2.

Steps 1–2 of Algorithm 2 are direct applications of Steps 1–2 of Algorithm 1
with inputs $\bm{\mathrm{D}},\bm{\mathrm{Y}}$ and
$\widetilde{\bm{\mathbb{Z}}}$. To be more specific, we use marginal screening
and joint thresholding to estimate relevant variables to the exposure within
the expanded set of $2p$ candidates. With slight abuse of notation, we use
$\check{\mathcal{S}}_{1}\subset\\{1,\ldots,2p\\}$ to denote the index set of
the $s$ candidate instruments selected from marginal screening. Let
$\widehat{\Gamma}_{\check{\mathcal{S}}_{1}}$ and
$\widehat{\gamma}_{\check{\mathcal{S}}_{1}}$ denote the debiased lasso
estimators based on the candidate instruments in $\check{\mathcal{S}}_{2}$.
The resulting estimated set of relevant instruments is
$\check{\mathcal{S}}_{2}=\big{\\{}1\leq l\leq
s:|\widehat{\gamma}_{\check{\mathcal{S}}_{1},l}|\geq\delta_{n}\times\text{SE}(\widehat{\gamma}_{\check{\mathcal{S}}_{1},l})\big{\\}},$
where $\delta_{n}=\sqrt{\omega\log(\max\\{n,s\\})}$ and the standard error
$\text{SE}(\widehat{\gamma}_{\check{\mathcal{S}}_{1},l})$ is obtained
similarly to (3.2).

######  Remark 4.1.

The tuning parameter $\omega$ plays an important role in our procedure. As we
demonstrated theoretically, larger values for $\omega$ lead to better
concentration of effect estimates from spurious instruments passing joint
thresholding. At the same time, large values of $\omega$ mean that the joint
thresholding step filters out more variables. In finite samples this can lead
to discarding some valid instruments, leading to efficiency loss in downstream
analyses.

In Step 3 of Algorithm 2, we identify the spurious instruments and remove them
from the candidate set, together with all the pseudo variables. Specifically,
we calculate the ratio estimate
$\widehat{\beta}^{(l)}=\widehat{\Gamma}_{\check{\mathcal{S}}_{1},l}/\widehat{\gamma}_{\check{\mathcal{S}}_{1},l}$
for each $l\in\check{\mathcal{S}}_{2}$. We then remove all the pseudo
variables from $\check{\mathcal{S}}_{2}$, together with all the candidate
instruments whose corresponding causal effect estimate $\widehat{\beta}^{(l)}$
falls into the range of causal effect estimates corresponding to the pseudo
variables. The resulting set of candidate instruments is denoted as
$\check{\mathcal{S}}_{3}$.

In Step 4 of Algorithm 2, we develop a mode-finding procedure to distinguish
the valid instruments from the invalid and remaining spurious instruments.
Similar to the literature on invalid IVs, we assume the plurality rule holds
in $\check{\mathcal{S}}_{3}$, i.e. the largest group of invalid or remaining
spurious instruments with the same Wald ratio is smaller than the group of
valid instruments. Recall that the population mode of
$\\{\beta^{(j)}:j\in\check{\mathcal{S}}_{3}\\}$ is the value that appears most
frequently, i.e.

$\text{mode}\\{{\beta}^{(l)}:l\in\check{\mathcal{S}}_{3}\\}=\mathop{\mathrm{argmax}}\limits_{{\beta}^{(l)}:l\in\check{\mathcal{S}}_{3}}|\\{j\in\check{\mathcal{S}}_{3}:{\beta}^{(j)}=\beta^{(l)}\\}|,$
(4.1)

where $|\cdot|$ denotes cardinality of a set and we define
${\beta}^{(j)}=\beta^{*}+C_{*}$ for indices $j$ corresponding to irrelevant
instruments. Under the population plurality rule the set $\mathcal{V}$ of
valid instruments in $\check{\mathcal{S}}_{3}$ is given by
$\mathcal{V}=\\{j\in\check{\mathcal{S}}_{3}:\beta^{(j)}=\text{mode}\\{{\beta}^{(l)}:l\in\check{\mathcal{S}}_{3}\\}\\}$.

To account for statistical uncertainty due to estimating ${\beta}^{(l)}$ in
(4.1), we say two estimates $\widehat{\beta}^{(j)}$ and
$\widehat{\beta}^{(l)}$ are close if
$|\widehat{b}^{(j,l)}|/\text{SE}(\widehat{b}^{(j,l)})$ is below a certain
threshold; here
$\widehat{b}^{(j,l)}=\widehat{\beta}^{(l)}-\widehat{\beta}^{(j)}$ is the
distance between the estimates $\widehat{\beta}^{(j)}$ and
$\widehat{\beta}^{(l)}$, and $\text{SE}(\widehat{b}^{(j,l)})$ is given in
(H.1) of the appendix. Following Donoho & Johnstone (1994), we recommend
$\sqrt{\omega^{2}\log(\max\\{n,s\\})}$ as the threshold. Here we use
$\omega^{2}$ to account for the number of pairwise comparisons. The estimated
set $\check{\mathcal{S}}_{4}$ of valid instruments consists of all candidate
instruments whose corresponding causal effect estimates is close to the
maximal number of other causal effect estimates as formalized in Step 4 of
Algorithm 2.

Once we have identified the set of valid IVs, it is straightforward to
estimate the causal effect using two-stage least squares. This is summarized
in Step 5 of Algorithm 2.

######  Remark 4.2.

The voting algorithm implemented in Section 3.2, originally developed by Guo
et al. (2018), uses a different approach than our Step 4 in Algorithm 2. It
treats each of the candidate instrument $Z_{l},l\in\check{\mathcal{S}}_{3}$ as
an expert, and each expert casts votes for other
$Z_{j},j\in\check{\mathcal{S}}_{3}.$ Importantly, this voting algorithm is not
symmetric in that expert $l$ casting vote for $j$ does not imply expert $j$
casting vote for $l$. In fact, under their algorithm, an effect estimate with
a larger variance is likely to receive more votes. This is undesirable as the
effect estimates with larger variances likely come from weak or spurious IVs.
The procedure we propose solves this problem under our setting.

INPUT: Design matrix $\bm{\mathbb{Z}}\in\mathbb{R}^{n\times p}$, observed
exposure $\bm{\mathrm{D}}\in\mathbb{R}^{n}$, observed outcome
$\bm{\mathrm{Y}}\in\mathbb{R}^{n}$.

  * 0.

Pseudo variables generation:

Generate $p$ pseudo variables $Z_{p+1},\ldots,Z_{2p}$ by randomly permuting
the rows of $\bm{\mathbb{Z}}$;

Denote the resulting design matrix as $\bar{\bm{\mathbb{Z}}}$, and the
concatenation of $\bm{\mathbb{Z}}$ and $\bar{\bm{\mathbb{Z}}}$ as
$\widetilde{\bm{\mathbb{Z}}}$.

  * 1+2

Marginal screening and joint thresholding:

Apply Steps 1 and 2 of Algorithm 1 with $\bm{\mathbb{Z}}\in\mathbb{R}^{n\times
p}$ replaced by $\widetilde{\bm{\mathbb{Z}}}\in\mathbb{R}^{n\times 2p}$. The
resulting set of variables that pass marginal and joint thresholding is
denoted by $\check{\mathcal{S}}_{2}$. Denote pseudos in
$\check{\mathcal{S}}_{2}$ by $\widehat{\mathcal{P}}$.

  * 3.

Irrelevant variables removal:

For each $l\in\check{\mathcal{S}}_{2}$, calculate a causal effect estimate
$\widehat{\beta}^{(l)}={\widehat{\Gamma}_{\check{\mathcal{S}}_{1},l}}/{\widehat{\gamma}_{\check{\mathcal{S}}_{1},l}};$

Remove variables similar to pseudo variables in their causal effect estimate
$\widehat{\beta}^{(l)}$ from $\check{\mathcal{S}}_{2}$:

$\check{\mathcal{S}}_{3}=\check{\mathcal{S}}_{2}\setminus\left\\{l\in\check{\mathcal{S}}_{2}:\min_{j\in\widehat{\mathcal{P}}}\widehat{\beta}^{(j)}\leq\widehat{\beta}^{(l)}\leq\max_{j\in\widehat{\mathcal{P}}}\widehat{\beta}^{(j)}\right\\}.$

  * 4.

Mode finding:

For $j,l\in\check{\mathcal{S}}_{3}$, compute
$\widehat{b}^{(j,l)}=\widehat{\beta}^{(l)}-\widehat{\beta}^{(j)}$ and
$\text{SE}(\widehat{b}^{(j,l)})$;

For $j\in\check{\mathcal{S}}_{3},$ find the number of candidates leading to a
similar causal effect estimate with $Z_{j}$:

$V_{j}=\Big{\lvert}\Big{\\{}l\in\check{\mathcal{S}}_{3}:|\widehat{b}^{(j,l)}|/\text{SE}(\widehat{b}^{(j,l)})\leq\sqrt{\omega^{2}\log(\max\\{n,s\\})}\Big{\\}}\Big{\rvert};$

Find valid IVs using the plurality rule:

$\check{\mathcal{S}}_{4}=\bigg{\\{}j:V_{j}=\max_{l\in\check{\mathcal{S}}_{3}}V_{l}\bigg{\\}}.$

  * 5.

Causal effect estimation: Take
$\widehat{\mathcal{S}}_{4},\bm{\mathrm{D}},\bm{\mathrm{Y}}$ as input into 2SLS
to estimate causal effect and construct confidence interval.

Algorithm 2 Fighting noise with noise

Algorithm 2 plays a key role in removing most irrelevant variables, but can
introduce selection bias into causal effect estimate. A classical approach to
deal with selection bias is sample splitting. Below, we discuss how sample
splitting can be implemented in the present context. Divide data into two
halves: $H_{1}$ of size $n_{1}$ and $H_{2}$ of size $n_{2}$. The first
subsample is used to identify and remove the irrelevant variables,
corresponding to Steps 0–3 of Algorithm 2. On the second subsample, we apply
an additional joint thresholding step to further estimate relevant variables
to the exposure. Then we use the proposed mode finding algorithm in Step 4 of
Algorithm 2 to identify the valid instruments and estimate the causal effect.
These steps are summarized in Algorithm 1, with details discussed in Section
A, both in the supplementary materials.

Although sample splitting can deal with selection bias in theory, it results
in smaller samples used for identifying valid instruments and in the final
estimation step. This can lead to reduced efficiency and imprecise inference
in smaller samples, which may outweigh the benefits of removing selection
bias. We observe this effect in some of our simulations reported in Section 5.

##  5 Simulation Studies

In this section, we run numerical experiments to evaluate the performance of
the proposed methods. The simulation setting here is the same as in the
motivating example in Section 3.2, except that here we consider various values
of the random error variance in the exposure $\sigma_{D}^{2}=0,2,4$. We
include the following methods for comparison:

  * •

“Proposed (full)” represents our proposed method in Algorithm 2. The lasso
penalization parameter $\lambda_{\gamma},\lambda_{\Gamma}$ are chosen equal to
$\sqrt{(\log p)/n}$, which is approximately 0.15. For the node-wise lasso, we
follow van de Geer et al. (2014) and choose the same tuning parameter
$\lambda_{j}$ using 10-fold cross-validation among all node-wise regressions.
We take $\omega=2.01$ (Donoho & Johnstone, 1994; Donoho, 1995).

  * •

“Proposed (split)” refers to the sample splitting method outlined at the end
of Section 4 and detailed in Algorithm 1 of the supplementary materials. We
choose $n_{1}=300$ and $n_{2}=200$. The tuning parameters are chosen similarly
to the “Proposed (full)”. To be more specific, on the first sample, we use
$\sqrt{(\log s)/n_{1}}$ in the joint threshold to estimate relevant variables,
while on the second sample, we use $\sqrt{(\log q)/n_{2}}$ in joint
thresholding, where $q$ denotes the number of remaining candidates after
removing spurious IVs.

  * •

“Naive” is the naive Algorithm 1. The tuning parameters
$\lambda_{\gamma},\lambda_{\Gamma}$ are chosen in the same way as “Proposed
(full)”.

  * •

“Oracle” is the 2SLS estimator based on valid instruments only, which provides
a benchmark for comparison.

Motivated by the data application in Section 6, for the implementation of the
methods “Proposed (full)”, “Proposed (split)” and “Naive”, we take $s=500$ in
the marginal screening step. We measure the performance through 1000 Monte
Carlo simulations. Simulation results regarding estimated bias, root-mean-
square error (RMSE) and empirical coverage are summarized in Table 1.

Table 1: Performance summary for various methods and $\sigma_{D}^{2}$ across
1000 Monte Carlo runs with $n=500$ and $p=50,000$. We make the numbers bold
for the method with the best performance except for the oracle estimate.
Standard errors are presented in the bracket if applicable. For simplicity, we
write $\text{SE}=0.00$ if $\text{SE}<0.005$

| Method |  $\text{Bias}\times 10$ ($\text{SE}\times 10$) | $\text{RMSE}\times 10$ | Coverage (Nominal = 95%)
---|---|---|---|---
| Proposed (full) |  $-$0.10(0.02) | 0.66 | 0.92
$\sigma_{D}^{2}=0$ | Proposed (split) |  $-$0.26(0.07) | 2.16 | 0.90
| Naive |  $-$2.23(0.02) | 2.28 | 0.01
| Oracle |  $-$0.02(0.01) | 0.28 | 0.94
| Proposed (full) |  $-$0.08(0.02) | 0.56 | 0.92
$\sigma_{D}^{2}=2$ | Proposed (split) |  $-$0.18(0.08) | 2.38 | 0.87
| Naive |  $-$1.82(0.02) | 1.94 | 0.06
| Oracle |  $-$0.02(0.01) | 0.27 | 0.94
| Proposed (full) |  $-$0.16(0.04) | 1.24 | 0.90
$\sigma_{D}^{2}=4$ | Proposed (split) |  $-$0.09(0.09) | 2.86 | 0.82
| Naive |  $-$1.09(0.02) | 1.28 | 0.21
| Oracle |  $-$0.02(0.01) | 0.27 | 0.94

The simulation results suggest that the proposed procedure using full data
outperforms competing methods in terms of bias, efficiency and validity of
inference. Similar to our findings in Section 3.2, estimates by “Naive” are
biased. Although sample splitting gives acceptable empirical coverage when
$\sigma_{D}^{2}$ is small, the bias and RMSE are much larger compared to the
ones obtained by the proposed method using full data. When the variance of the
random error in the exposure $\epsilon_{D}$ increases, the performance of
“Proposed (full)” remains similar, which indicates that our approach using
full data works even when the variance of unobserved confounding is relatively
small compared to that of the random error $\epsilon_{D}$. For the sample
splitting method, one can see the RMSE increases and the coverage probability
decreases further below the nominal level when $\epsilon_{D}$ increases.

To better understand the performance of various methods as summarized in Table
1, we report the average number of valid, invalid and spurious instruments in
$\check{\mathcal{S}}_{2},\check{\mathcal{S}}_{3}$ and
$\check{\mathcal{S}}_{4}$ respectively in Table 2.

Table 2:  Average number of invalid, valid IVs and irrelevant variables in
$\check{\mathcal{S}}_{2}$ (the candidate set after joint thresholding),
$\check{\mathcal{S}}_{3}$ (the candidate set after removing spurious
instruments) and $\check{\mathcal{S}}_{4}$ (the estimated set of valid IVs),
respectively. Note that “Naive” is applied without including pseudo variables.
The results are based on $n=500,$ $p=50,000$ and 1000 Monte Carlo runs. “$-$”
stands for “inapplicable”

|  | $\check{\mathcal{S}}_{2}$ | $\check{\mathcal{S}}_{3}$ | $\check{\mathcal{S}}_{4}$
---|---|---|---|---
| Method | Invalid | Valid | Irrelevant |  | Invalid | Valid | Irrelevant |  | Invalid | Valid | Irrelevant
| Proposed (full) | 2.00 | 7.00 | 30.53 |  | 1.99 | 6.99 | 1.87 |  | 0.00 | 6.69 | 0.51
$\sigma_{D}^{2}=0$ | Proposed (split) | 1.93 | 7.00 | 48.12 |  | 1.84 | 6.74 | 2.08 |  | 0.19 | 6.73 | 0.07
| Naive | 2.00 | 7.00 | 15.26 |  | $-$ | $-$ | $-$ |  | 0.00 | 0.02 | 15.25
| Proposed (full) | 2.00 | 7.00 | 35.19 |  | 1.97 | 6.99 | 1.95 |  | 0.01 | 6.88 | 1.04
$\sigma_{D}^{2}=2$ | Proposed (split) | 1.92 | 7.00 | 53.37 |  | 1.80 | 6.46 | 2.01 |  | 0.24 | 6.42 | 0.05
| Naive | 2.00 | 7.00 | 17.67 |  | $-$ | $-$ | $-$ |  | 0.00 | 0.60 | 17.66
| Proposed (full) | 2.00 | 7.00 | 39.55 |  | 1.98 | 6.61 | 2.00 |  | 0.04 | 6.50 | 1.16
$\sigma_{D}^{2}=4$ | Proposed (split) | 1.92 | 7.00 | 58.00 |  | 1.81 | 5.29 | 2.01 |  | 0.33 | 5.19 | 0.12
| Naive | 2.00 | 7.00 | 20.45 |  | $-$ | $-$ | $-$ |  | 0.00 | 2.93 | 20.45

The results show that “Naive” misidentifies spurious instruments as valid and
misses many or almost all the valid instruments. The sample splitting
procedure can find most valid instruments when the variance $\sigma_{D}^{2}$
is small. However, the estimated set of valid instruments will include some
invalid ones since we only use part of the data and the estimates have larger
variation than those obtained using the full sample. Consequently, the
resulting causal effect estimate has moderate bias and the empirical coverage
deviates from the nominal level. In comparison, using our proposed method on
the full sample can locate spurious instruments with the help of pseudo
variables and correctly identify most valid instruments even if
$\sigma_{D}^{2}$ is large.

To provide further insights, we plot the distribution of
$\widehat{\beta}^{(j)}$ collected from all Monte Carlo runs for valid, invalid
and spurious instruments in $\check{\mathcal{S}}_{2}$,
$\check{\mathcal{S}}_{3}$ and $\check{\mathcal{S}}_{4}$ respectively in
(a)-(c) of Figure 3, obtained from “Proposed (full)” under $\sigma_{D}^{2}=0$.
We can see that causal effect estimates from spurious instruments and valid
instruments are separable. Most spurious instruments lie in the region formed
by causal effect estimates from pseudo variables, and thus are removed by
following the proposed procedure. Thereafter the collection of valid
instruments satisfies the plurality rule empirically, so the mode-finding step
is able to identify the valid instruments.

(a) $\check{\mathcal{S}}_{2}$ (b) $\check{\mathcal{S}}_{3}$: Remove spurious
IVs (c) $\check{\mathcal{S}}_{4}$ (d) $\check{\mathcal{S}}_{4}$ (refit)

Figure 3: Plots (a)-(d) correspond to the simulation in Section 5: Plots (a)
to (c) are histograms of causal effect estimates from valid, invalid and
spurious instruments in $\check{\mathcal{S}}_{2}$, $\check{\mathcal{S}}_{3}$
and $\check{\mathcal{S}}_{4}$ respectively obtained from Algorithm 2 across
1000 Monte Carlo runs. The blue pile represents causal effect estimates by
valid IVs; the green piles represent those by invalid IVs; the orange pile
represents those by irrelevant variables in the initial candidate set; the
black pile represents those by pseudo variables. The red dashed line is the
true causal effect $\beta^{*}=2$. The results are based on
$(n,p,\sigma_{D}^{2})=(500,50,000,0)$. Plot (d) gives the distribution of
effect estimates from refitted model.

To further evaluate the performance of the proposed method, we consider a
variety of additional settings in Section F of the supplementary materials. In
particular, under the same data generating setting here, we evaluate our
method when there are much more noise variables to start with. Additionally,
we consider scenarios where irrelevant variables are correlated with valid
instruments, some irrelevant variables are correlated with each other, and
simulations under multivariate unmeasured confounders. We also carry out
numerical experiments, in which the “independence” IV assumption (I3) is
violated. Furthermore, since irrelevant variable can be seen as an extreme
case of the weak instrument, we compare with competing methods (Donald &
Newey, 2001; Carrasco, 2012; Carrasco & Tchuente, 2015) in the econometric
literature that deal with many IVs where some are possibly weak.

##  6 Real Data Application

We apply the proposed procedure to evaluate the effect of obesity on Health-
Related Quality of Life (HRQL). Obesity increases the risks of many health
problems, such as heart disease, diabetes, high blood pressure, and certain
cancers. Mechanical stress resulting from obesity can also impair mental well-
being (Myers & Rosen, 1999; Katz et al., 2000). Health-related quality of life
(HRQL) is emerging as an important outcome in obesity studies, which refers to
the overall effects of medical conditions on physical and mental functioning
and well-being as subjectively evaluated and reported by the patient (Fontaine
et al., 1996). The development of standardized measures of HRQL helps us
better understand the impact of a given illness on physical and mental health,
and further influences public health policies.

To examine the causal relationship between obesity and HRQL, we use data from
the Wisconsin Longitudinal Study (WLS). The WLS genetic data is sponsored by
the National Institute on Aging (grant numbers R01AG009775, R01AG033285, and
R01AG041868) and was conducted by the University of Wisconsin. More
information can be found at https://www.ssc.wisc.edu/wlsresearch/. Study
participants were graduates from Wisconsin high schools in 1957 and their
siblings (Herd et al., 2014). Our analysis focuses on unrelated graduates
reinterviewed in 2011. We take BMI as the exposure, a measure commonly used to
define obesity. It ranges from 12 to 64 in the data set. The healthy range for
BMI is between 18.5 and 24.9. As shown in Figure 6 in the supplementary
materials, for underweight and normal individuals, HRQL is positively
correlated with BMI, while for overweight or obese ones, HRQL is negatively
correlated with BMI. Motivated by this observation, We shall restrict our
analysis to participants whose BMI is above 25. There are 3023 subjects in
total, 51% of which are females. In our analysis, the average BMI is 30.6 with
standard deviation 4.93. To measure HRQL, we use the Health Utility Index Mark
3 (HUI-3) which takes values between $-0.22$ and 1. HUI-3 measured on the
subjects in our analysis has mean 0.786 and standard deviation 0.227. We
adjust for observed covariates including age, gender and year of education.
Among participants considered in our analysis, the average age is 71.2 years
old with standard deviation 0.9; the average education length is 13.8 years
with standard deviation 2.38. We also adjust for the top six principal
components to account for population stratification (Patterson et al., 2006;
Sanderson et al., 2021), which are obtained by performing principal components
analysis (PCA) on unrelated study participants.

We begin with a crude analysis that ignores the unmeasured confounding. Simply
regressing the HUI-3 on BMI gives the OLS estimate $-0.011$ with 95% CI
$[-0.013,-0.009]$. However, unmeasured factors such as smoking status and the
use of alcohol may confound the relationship between BMI and HUI-3. In
response to this potential problem, we use the MR method to investigate the
causal relationship between these two. After quality control (QC) filtering of
SNPs, we get a candidate set of 3,683,868 genetic variants in our MR analysis.
Details of the QC steps are discussed in Section B.2 of the supplementary
materials. Following Algorithm 2, we generate pseudo variables by randomly
permuting the SNP design matrix by rows, resulting in a candidate set of
7,367,736 variants in total.

We first apply a marginal screening step to quickly reduce the dimension of
initial candidate instruments via ranking the marginal correlation in absolute
value between each candidate instrument and the exposure. A common choice is
to keep $\lfloor n/\log(n)\rfloor$ candidates (Fan & Lv, 2008), which is
approximately 460 in our analysis. We round it to the nearest multiple of 100.
Specifically, in the marginal screening step, by ranking the absolute marginal
strength of association with BMI, we keep the top 500 candidate instruments,
denoted by $\check{\mathcal{S}}_{1}(\text{Proposed})$. There were 268 true
variants and 232 pseudo ones included in this set. To further identify genetic
variants related to BMI, we fit a joint model for the 500 candidate
instruments.

We visualize in Figure 4(a) the result from our proposed method. Specifically,
a total of 75 SNPs ($\check{\mathcal{S}}_{2}(\text{Proposed})$) pass the joint
thresholding, corresponding to all the points shown in Figure 4(a). The 34
points on the left of the orange dashed line correspond to true genetic
variants, while the 37 red dots on the right correspond to pseudo variants. In
$\check{\mathcal{S}}_{2}(\text{Proposed})$, causal effect estimates from 71
SNPs fall into the region formed by pseudo variants, which is represented by
the band between the two black dashed lines. These SNPs are estimated as
spurious and are removed from the candidate set. The remaining 4 points
outside the region are genetic variants that are estimated to be relevant to
BMI ($\check{\mathcal{S}}_{3}(\text{Proposed})$). An application of the
proposed mode-finding procedure on $\check{\mathcal{S}}_{3}(\text{Proposed})$
gives us 3 valid instruments ($\check{\mathcal{S}}_{4}(\text{Proposed})$),
corresponding to the blue squares in Figure 4(a): “rs144468788” on chromosome
2 and in the MIR4432 host gene (MIR4432HG); “rs113116137” and “rs2546230” on
chromosome 8, and “rs113116137” is in the CUB And Sushi Multiple Domains 1
(CSMD1) gene. The purple point at the bottom represents the estimated invalid
instrument. Finally, we fit two-stage least squares using the three identified
valid IVs and obtain the causal effect estimate $\widehat{\beta}=-0.039$ with
95% confidence interval $[-0.052,-0.025]$. This suggests that in the
overweight or obese population, one unit increase in BMI will result in 0.039
unit decrease of HUI-3. In other words, one standard deviation increase in BMI
leads to roughly one standard deviation decrease in health-related quality of
life.

(a) Proposed Algorithm 2 (b) Naive Algorithm 1 (c) Naive Algorithm 1: pseudo
included

Figure 4: Causal effect estimates from estimated relevant SNPs versus SNP
index. (a): Causal effect estimates associated with genetic variants in
$\check{\mathcal{S}}_{2}(\text{Proposed})$. The orange dot-dashed line
separates pseudo and real genetic variants. Red dots represent pseudo SNPs and
all the others are true genetic variants. The two black dashed lines represent
the estimated range. The three blue squares represent valid instruments
identified through our procedure, and the gray diamond represents the invalid
instrument we identified. (b): Causal effect estimates from relevant SNPs
estimated by Algorithm 1 (“Naive”) which is applied to true variants only. The
blue squares represent valid IVs identified by “Naive” and green triangles
represent estimated invalid IVs. (c) Results from Algorithm 1 with pseudos
included. The orange dot-dashed line separates pseudo and real genetic
variants. Green triangles refer to SNPs identified as invalid. The blue
squares are SNPs identified as valid instruments.

To explore that the standard GWAS threshold might be too stringent to retain
genetic variants associated with the trait of interest in our data
application, we apply the genome-wide significance $p$-value threshold of
$5\times 10^{-8}$ to the true genetic variants to identify significant
associations with BMI. This gives us only one SNP “rs148605555” on chromosome
9. Although the SNP is estimated as relevant to BMI using our procedure, that
is, the SNP is in the set $\check{\mathcal{S}}_{2}(\text{Proposed})$, the
associated effect estimate given by $-0.008$ lies in the middle of the
estimated range by pseudos and thus the SNP is identified as spurious in the
proposed method. In addition, the effect estimate of the SNP that passes the
GWAS threshold is close to the OLS estimate, an indication that this might be
a spurious one.

For comparison, we also apply Algorithm 1, denoted by “Naive”, to this data
set. The results are visualized in Figure 4(b). Specifically, we apply
marginal screening to the true genetic variants only and keep the top 500
($\widehat{\mathcal{S}}_{1}(\text{Naive})$) by ranking absolute marginal
strength concerning the exposure. The joint thresholding step keeps 44
variants $(\widehat{\mathcal{S}}_{2}(\text{Naive}))$, and 29 of them are
estimated as valid IVs $(\widehat{\mathcal{S}}_{4}(\text{Naive}))$. Comparing
the candidate instruments retained after the joint thresholding, we note that
$\check{\mathcal{S}}_{2}(\text{Proposed})$ and
$\widehat{\mathcal{S}}_{2}(\text{Naive})$ have 15 variants in common, and only
one common variant “rs144468788”, corresponding to the point at the bottom
left corner in Figure 4(b), is also estimated as relevant by “Proposed”, i.e.
$\\{\text{rs144468788}\\}\subset\check{\mathcal{S}}_{3}(\text{Proposed})$. The
“rs144468788” is estimated as invalid by “Naive”, but valid by “Proposed”,
i.e.,
$\widehat{\mathcal{S}}_{4}(\text{Naive})\cap\check{\mathcal{S}}_{4}(\text{Proposed})=\emptyset$.
Interestingly, the resulting causal effect estimate by “Naive”, that is
$-0.010$ (95% CI $[-0.015,-0.005]$), is roughly the same as the OLS estimate,
equivalently saying that one standard deviation increase in BMI leads to
approximately 0.2 standard deviation decrease in health-related quality of
life.

To further understand the behaviour of the naive Algorithm 1 with this data
set, we apply this method to the candidate set
$\check{\mathcal{S}}_{1}(\text{Proposed})$ obtained after marginal screening
in our procedure. In doing this, we pretend that the analyst does not know
which variables in $\check{\mathcal{S}}_{1}(\text{Proposed})$ are pseudo
variables. Analysis results in Figure 4(c) show that as expected, “Naive”
misidentifies many pseudo variables as valid instruments. The resulting causal
effect estimate, $-0.012$ (95% CI $[-0.016,-0.008]$), which is close to the
OLS estimate. That is to say, one standard deviation increase in BMI roughly
leads to 0.2 standard deviation decrease in health-related quality of life.

In our data application, we did not implement the sample splitting Algorithm 1
for two main reasons: for one thing, the sample splitting method suffers from
power loss, which is also demonstrated in our numerical experiments. And for
another, the SNP matrix is encoded as 0,1,2, which represents possible
genotypes. Splitting the data might dramatically change the genotype
distribution since possibly some genotypes appear more often in one portion of
the data but not the other (Barber & Candes, 2019).

We close this section by performing additional numerical studies to explore
the uncertainty of the aforementioned results, due to the randomness induced
by generating pseudo variants. We run the procedure 50 times by repeatedly
permuting the true SNP design matrix to generate pseudos. The median of causal
effect estimates is $-0.037$. For comparison, we also repeat Algorithm 1 to
the same candidate set used in our procedure 50 times. The median of the
resulting effect estimates is $-0.010$.

##  7 Discussion

To make causal inference with many candidate instruments, it is common to
first reduce the dimension of candidate instruments with screening procedures.
If the number of candidate instruments is large, and/or the sample size is
small to moderate, there may be more spurious instruments than valid
instruments in the estimated set of variables relevant to the exposure. In
this paper, we show that this can lead to biased causal inference. In
particular, the resulting causal effect estimate may be close to the OLS
estimate that does not account for any unmeasured confounding. We also develop
a novel strategy that constructs pseudo variables to find the spurious
instruments.

We note that spurious instruments may not be a practical problem if the number
of candidate instruments is moderate, or if the sample size is very large. For
example, this may not be a major concern in large-scale GWAS with hundreds of
thousands of samples.

Our work is related to a strand of research in econometrics on consistent
estimation in instrumental variables regression with many (possibly weak)
instruments (e.g. Chao & Swanson, 2005; Hansen et al., 2008; Carrasco &
Tchuente, 2015). The setup they consider only involves valid instruments,
while our method deals with the potential irrelevant and invalid instruments
in the candidate instruments set. Our method may also be extended to deal with
weak instruments. On a high level, a weak instrument behaves like a valid,
invalid or spurious instrument depending on the rate of convergence of the
instrumental strength. Under the simple setting as displayed in Figure 5,
where the instrument $Z$ is weakly correlated with the exposure $D$ and is
independent of the unobserved confounding $U$, Nelson & Startz (1990b) has
shown that the Wald estimate satisfies the expansion below

$\widehat{\beta}-\beta^{*}=\frac{\widehat{\mathrm{cov}}(Z,U)}{\frac{1}{\alpha^{*}}\widehat{\mathrm{cov}}(Z,U)+\gamma^{*}}.$
(7.1)

Here $\widehat{\mathrm{cov}}$ denotes the sample covariance. The limit of the
right-hand side of (7.1) depends on the rate at which the relevance parameter
$\gamma^{*}$ converges to 0. Specifically, if the covariance term
$\widehat{\mathrm{cov}}(Z,U)$ converges in probability to zero at $n^{-1/2}$
rate, then the bias of IV estimate concentrates at $\alpha^{*}$ when
$\gamma^{*}=0$ or $\gamma^{*}\to 0$ faster than $n^{-1/2}$. If $\gamma^{*}\to
0$ slower than $n^{-1/2}$, then the bias is asymptotically zero and $Z$
behaves like a valid instrument. When $\gamma^{*}\to 0$ at the $n^{-1/2}$
rate, or mathematically, $\gamma^{*}=cn^{-1/2}$ for some constant $c>0$, $Z$
behaves like an invalid instrument with the degree of of violation determined
by $c$. Extending this result for effect estimation after screening and
thresholding is an interesting topic for future research.

$U$$D$$Y$$Z$$1/\alpha^{*}$1$\beta^{*}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\gamma^{*}}$
Figure 5: Graphical illustration of an instrumental variable model: $Z$
denotes a weak instrument, $D,Y$ represent the exposure and outcome,
respectively and $U$ is the unobserved confounding. Parameters on the arrows
represent effects on the linear scale. The red dashed line refers to sample
correlation.

Our framework is also applicable under the two-sample design where the
exposure and the outcome are collected from two independent samples, say
$\Lambda_{1}$ and $\Lambda_{2}$ of sizes $n_{1}$ and $n_{2}$ respectively.
Here effect estimates from spurious instruments are still separable from valid
instruments. But different from the one-sample design, now effect estimates
from selected irrelevant variables are centered around zero since these
variables exhibit spurious correlation with the exposure only on $\Lambda_{1}$
but not with the outcome on $\Lambda_{2}$. For better illustration, we
consider a simple setting where $D=\alpha_{D}^{*}U$ and
$Y=\beta^{*}D+\alpha_{Y}^{*}U$ and assume $n_{1},n_{2}$ are of the same order.
Suppose only one variable $Z_{j_{0}}$ survives the marginal screening, which
implies that the spurious correlation between $Z_{j_{0}}$ and
$D_{\Lambda_{1}}$ is the maximal correlation between $D_{\Lambda_{1}}$ and
candidate $Z_{j}$ for all $1\leq j\leq p$. One can show that the effect of
$Z_{j_{0}}$ on the exposure is of order $\sqrt{(\log p)/n_{1}}$, while the
effect of $Z_{j_{0}}$ on the outcome converges to zero faster at rate
$1/\sqrt{n_{2}}$. We leave this extension for future work.

In addition, our current framework only considers invalid instruments that
have a direct effect on the outcome. Empirical evidence (see Section F.6 of
the supplementary materials for more details) demonstrates that our method is
robust to the violation of independence assumption (I3), it would be of
interest to extend our theory and validate the proposed work in more general
contexts. It is also worth exploring whether our findings can be adapted to
models with multiple treatments.

## References

  * WLS (????) Wisconsin Longitudinal Study. Produced and distributed by the University of Wisconsin with funding from the National Institute on Aging (grant numbers R01AG009775, R01AG033285, and R01AG041868), Madison, WI.
  * Angrist & Keueger (1991) Angrist, J. D. & Keueger, A. B. (1991).  Does compulsory school attendance affect schooling and earnings?  The Quarterly Journal of Economics 106, 979–1014.
  * Barber & Candes (2015) Barber, R. F. & Candes, E. J. (2015).  Controlling the false discovery rate via knockoffs.  Ann. Statist. 43, 2055–2085.
  * Barber & Candes (2019) Barber, R. F. & Candes, E. J. (2019).  A knockoff filter for high-dimensional selective inference.  Ann. Statist. 47, 2504–2537.
  * Bound et al. (1995) Bound, J., Jaeger, D. A. & Baker, R. M. (1995).  Problems with instrumental variables estimation when the correlation between the instruments and the endogeneous explanatory variable is weak.  J Am Statist Assoc. 90, 443–450.
  * Bowden et al. (2015) Bowden, J., Smith, G. D. & Burgess, S. (2015).  Mendelian randomization with invalid instruments: effect estimation and bias detection through Egger regression.  Int. J. Epidem. 44, 512–525.
  * Bowden et al. (2016) Bowden, J., Smith, G. D., Haycock, P. C. & Burgess, S. (2016).  Consistent estimation in Mendelian randomization with some invalid instruments using a weighted median estimator.  Genet. Epidem. 40, 304–314.
  * Browning (2008) Browning, S. (2008).  Missing data imputation and haplotype phase inference for genome-wide association studies.  Hum Genet. 124, 439–450.
  * Bühlmann & van de Geer (2011) Bühlmann, P. & van de Geer, S. (2011).  Statistics for High-Dimensional Data.  Springer Series in Statistics. Springer.
  * Burgess et al. (2015) Burgess, S., Timpson, N. J., Ebrahim, S. & Smith, G. D. (2015).  Mendelian randomization: where are we now and where are we going?  Int J Epidemiol. 44, 379–388.
  * Carrasco (2012) Carrasco, M. (2012).  A regularization approach to the many instruments problem.  J. Econom. 170, 383–398.
  * Carrasco & Tchuente (2015) Carrasco, M. & Tchuente, G. (2015).  Regularized LIML for many instruments.  J. Econom. 186, 427–442.
  * Chao & Swanson (2005) Chao, J. C. & Swanson, N. R. (2005).  Consistent estimation with a large number of weak instruments.  Econometrica 73, 1673–1692.
  * Cheng & Liao (2015) Cheng, X. & Liao, Z. (2015).  Select the valid and relevant moments: An information-based LASSO for GMM with many moments.  J. Econom. 186, 443–464.
  * de Bakker et al. (2008) de Bakker, P. I., Ferreira, M. A., Jia, X., Neale, B. M., Raychaudhuri, S. & Voight, B. F. (2008).  Practical aspects of imputation-driven meta-analysis of genome-wide association studies.  Hum Mol Genet. 17, R122–R128.
  * Didelez & Sheehan (2007) Didelez, V. & Sheehan, N. (2007).  Mendelian randomization as an instrumental variable approach to causal inference.  Stat Methods Med Res. 16, 309–330.
  * Donald & Newey (2001) Donald, S. G. & Newey, W. K. (2001).  Choosing the number of instruments.  Econometrica. 69, 1161–1191.
  * Donoho (1995) Donoho, D. (1995).  De-noising by soft-thresholding.  IEEE Trans. Inform. Theory 41, 613–627.
  * Donoho & Johnstone (1994) Donoho, D. L. & Johnstone, I. M. (1994).  Ideal spatial adaptation by wavelet shrinkage.  Biometrika. 81, 425–455.
  * Eckles et al. (2016) Eckles, D., Kizilcec, R. F. & Bakshy, E. (2016).  Estimating peer effects in networks with peer encouragement designs.  Proceedings of the National Academy of Sciences 113, 7316–7322.
  * Fan & Lv (2008) Fan, J. & Lv, J. (2008).  Sure independence screening for ultrahighdimensional feature space.  J. R. Statist. Soc. B. 70, 849–911.
  * Fontaine et al. (1996) Fontaine, K. R., Cheskin, L. J. & Barofsky, I. (1996).  Health-related quality of life in obese persons seeking treatment.  J Fam Pract. 43, 265–270.
  * Goldman & Rao (2016) Goldman, M. & Rao, J. (2016).  Experiments as instruments: Heterogeneous position effects in sponsored search auctions.  EAI Endorsed Trans. Serious Games 3, e2.
  * Guo et al. (2018) Guo, Z., Kang, H., Cai, T. T. & Small, D. S. (2018).  Confidence intervals for causal effects with invalid instruments by using two‐stage hard thresholding with voting.  J. R. Statist. Soc. B. 80, 793–815.
  * Han (2008) Han, C. (2008).  Detecting invalid instruments using $l_{1}$-GMM.  Econ. Lett. 101, 285–287.
  * Hansen et al. (2008) Hansen, C., Hausman, J. & Newey, W. (2008).  Estimation with many instrumental variables.  J. Bus. Econ. Statist. 26, 398–422.
  * Hartford et al. (2021) Hartford, J. S., Veitch, V., Sridhar, D. & Leyton-Brown, K. (2021).  Valid causal inference with (some) invalid instruments.  ICML. 139, 4096–4106.
  * Herd et al. (2014) Herd, P., Carr, D. & Roan, C. (2014).  Cohort Profile: Wisconsin Longitudinal Study (WLS).  Int. J. Epidemiol. 43, 34–41.
  * Howie et al. (2011) Howie, B., Marchini, J. & Stephens, M. (2011).  Genotype imputation with thousands of genomes.  G3 (Bethesda). 1, 457–470.
  * Howie et al. (2009) Howie, B. N., Donnelly, P. & Marchini, J. (2009).  A flexible and accurate genotype imputation method for the next generation of genome-wide association studies.  PLoS Genet 5, e1000529.
  * Javanmard & Montanari (2014) Javanmard, A. & Montanari, A. (2014).  Confidence intervals and hypothesis testing for high-dimensional regression.  J. Mach. Learn. Res. 15, 2869–2909.
  * Kang et al. (2016) Kang, H., Zhang, A., Cai, T. T. & Small, D. S. (2016).  Instrumental variables estimation with some invalid instruments and its application to Mendelian randomization.  J. Am. Statist. Ass. 111, 132–144.
  * Katz et al. (2000) Katz, D. A., McHorney, C. A. & Atkinson, R. L. (2000).  Impact of obesity on health-related quality of life in patients with chronic illness.  J Gen Intern Med. 15, 789–796.
  * Kolesar et al. (2015) Kolesar, M., Chetty, R., Friedman, J., Glaeser, E. & Imbens, G. W. (2015).  Identification and inference with many invalid instruments.  J. Bus. Econ. Statist. 33, 474–484.
  * Laurie et al. (2010) Laurie, C. C., Doheny, K. F., Mirel, D. B. et al. (2010).  Quality control and quality assurance in genotypic data for genome-wide association studies.  Genet Epidemiol. 34, 591–602.
  * Lawlor et al. (2008) Lawlor, D. A., Harbord, R. M., Sterne, J. A. C., Timpson, N. & Smith, G. D. (2008).  Mendelian randomization: Using genes as instruments for making causal inferences in epidemiology.  Statist. Med. 27, 1133–1163.
  * Li et al. (2009) Li, Y., Willer, C., Sanna, S. & Abecasis, G. (2009).  Genotype imputation.  Annu Rev Genomics Hum Genet. 10, 387–406.
  * Li et al. (2010) Li, Y., Willer, C. J., Ding, J., Scheet, P. & Abecasis, G. R. (2010).  MaCH: Using sequence and genotype data to estimate haplotypes and unobserved genotypes.  Genet Epidemiol. 34, 816–834.
  * Marchini et al. (2007) Marchini, J., Howie, B., Myers, S., McVean, G. & Donnelly, P. (2007).  A new multipoint method for genome-wide association studies by imputation of genotypes.  Nat. Genet. 39, 906–913.
  * Morrison et al. (2020) Morrison, J., Knoblauch, N., Marcus, J. H. et al. (2020).  Mendelian randomization accounting for correlated and uncorrelated pleiotropic effects using genome-wide summary statistics.  Nat. Genet. 52, 740—747.
  * Myers & Rosen (1999) Myers, A. & Rosen, J. C. (1999).  Obesity stigmatization and coping: relation to mental health symptoms, body image, and self-esteem.  Int J Obes Relat Metab Disord. 23, 221–230.
  * Nelson & Startz (1990a) Nelson, C. R. & Startz, R. (1990a).  The distribution of the instrumental variables estimator and its t-ratio when the instrument is a poor one.  J. Bus. 63, S125–S140.
  * Nelson & Startz (1990b) Nelson, C. R. & Startz, R. (1990b).  Some further results on the exact small sample properties of the instrumental variable estimator.  Econometrica. 58, 967–976.
  * Patterson et al. (2006) Patterson, N., Price, A. L. & Reich, D. (2006).  Population structure and eigenanalysis.  PLoS Genet. 2, e190.
  * Peysakhovich & Eckles (2018) Peysakhovich, A. & Eckles, D. (2018).  Learning causal effects from many randomized experiments using regularized instrumental variables.  In Proceedings of the 2018 World Wide Web Conference.
  * Sanderson et al. (2021) Sanderson, E., Richardson, T. G., Hemani, G. & Smith, G. D. (2021).  The use of negative control outcomes in Mendelian randomization to detect potential population stratification.  Int J Epidemiol. 50, 1350–1361.
  * Scheet & Stephens (2006) Scheet, P. & Stephens, M. (2006).  A fast and flexible statistical model for large-scale population genotype data: Applications to inferring missing genotypes and haplotypic phase.  Am. J. Hum. Genet. 78, 629–644.
  * Smith (2007) Smith, G. D. (2007).  Capitalizing on Mendelian randomization to assess the effects of treatments.  J R Soc Med. 100, 432–435.
  * Stock et al. (2002) Stock, J. H., Wright, J. H. & Yogo, M. (2002).  A survey of weak instruments and weak identification in generalized method of moments.  J. Bus. Econ. Statist. 20, 518–529.
  * Stock & Yogo (2005) Stock, J. H. & Yogo, M. (2005).  Testing for Weak Instruments in Linear IV Regression, chap. 5.  Cambridge University Press, pp. 80–108.
  * van de Geer et al. (2014) van de Geer, S., Bühlmann, P., Ritov, Y. & Dezeure, R. (2014).  On asymptotically optimal confidence regions and tests for high-dimensional models.  Ann. Statist. 42, 1166–1202.
  * Vershynin (2012) Vershynin, R. (2012).  Introduction to the non-asymptotic analysis of random matrices, chap. 5.  Cambridge Univ. Press.
  * Vershynin (2018) Vershynin, R. (2018).  High-Dimensional Probability: An Introduction with Applications in Data Science High-Dimensional Probability: An Introduction with Applications in Data Science.
  * Wainwright (2019) Wainwright, M. J. (2019).  High-Dimensional Statistics: A Non-Asymptotic Viewpoint.  Cambridge University Press.
  * Wang & Tchetgen Tchetgen (2018) Wang, L. & Tchetgen Tchetgen, E. (2018).  Bounded, efficient and multiply robust estimation of average treatment effects using instrumental variables.  J. R. Statist. Soc. B. 80, 531–550.
  * Windmeijer et al. (2019) Windmeijer, F., Farbmacher, H., Davies, N. & Smith, G. D. (2019).  On the use of the lasso for instrumental variables estimation with some invalid instruments.  J. Am. Statist. Ass. 114, 1339–1350.
  * Windmeijer et al. (2021) Windmeijer, F. A., Liang, X., Hartwig, F. P. & Bowden, J. (2021).  The confidence interval method for selecting valid instrumental variables.  J. R. Statist. Soc. B. 83, 752–776.
  * Wu et al. (2007) Wu, Y., Boos, D. D. & Stefanski, L. A. (2007).  Controlling variable selection by the addition of pseudovariables.  J. Am. Statist. Ass. 102, 235–243.
  * Zhao et al. (2020) Zhao, Q., Wang, J., Hemani, G., Bowden, J. & Small, D. S. (2020).  Statistical inference in two-sample summary-data Mendelian randomization using robust adjusted profile score.  Ann. Statist. 48, 1742–1769.

The supplementary file is organized as follows. We introduce the sample
splitting algorithm in Section A. More details on our real data application
are provided in Section B. We introduce the notations used throughout the
paper in Section C. Technical details are collected in Section D. We discuss
in Section E the derivation of standard error for bias estimate. Additional
simulations to demonstrate the robustness of the proposed method are presented
in Section F.

##  Appendix A Alternative Algorithm by Sample Splitting

To deal with selection bias and perform statistical inference on the causal
effect estimate, we consider the sample splitting version of the proposed
method. Specifically, we split the sample into two halves: $H_{1}$ and
$H_{2}$. The sizes of $H_{1}$ and $H_{2}$ are denoted by $n_{1}$ and $n_{2}$
respectively, which satisfy $n_{1}+n_{2}=n$ and $n_{1}=cn_{2}$ for some
constant $c>0$. In our numerical experiments, we take $c=1.5$. Let
$\bm{\mathbb{Z}}^{(1)}\in\mathbb{R}^{n_{1}\times
p},\bm{\mathrm{D}}^{(1)}\in\mathbb{R}^{n_{1}}$ and
$\bm{\mathrm{Y}}^{(1)}\in\mathbb{R}^{n_{1}}$ denote the collection of data on
$H_{1}$, and correspondingly $\bm{\mathbb{Z}}^{(2)}\in\mathbb{R}^{n_{2}\times
p},\bm{\mathrm{D}}^{(2)}\in\mathbb{R}^{n_{2}}$ and
$\bm{\mathrm{Y}}^{(2)}\in\mathbb{R}^{n_{2}}$ represent the data on $H_{2}$. We
use $H_{1}$ to find spurious instruments. We generate $p$ pseudo variables on
$H_{1}$ by randomly permuting the rows of $\bm{\mathbb{Z}}^{(1)}$. Let
$\bar{\bm{\mathbb{Z}}}^{(1)}$ be the design matrix associated with pseudo
variables. Combine the $p$-dimensional pseudo variables with the initial
$p$-dimensional candidate instruments, and we use
$\widetilde{\bm{\mathbb{Z}}}^{(1)}\in\mathbb{R}^{n_{1}\times 2p}$ to denote
the resulting design matrix defined on $H_{1}$. Identification of valid IVs
and estimation of causal effect together with inference are performed on
$H_{2}$. The procedure with sample splitting is summarized in Algorithm 1. We
discuss this alternative method in details in Sections A.1 and A.2.

###  A.1 Estimation of valid instruments

We first apply marginal screening to the subsample $H_{1}$ to reduce the
dimension of $\widetilde{\bm{Z}}^{(1)}\in\mathbb{R}^{2p}$, which is
concatenation of the initial $p$-dimensional candidate instruments and pseudo
variables. Let $\check{\mathcal{S}}_{1}\in\\{1,\ldots,2p\\}$ be the index set
of candidates selected by marginal screening and assume that
$\check{\mathcal{S}}_{1}$ include all the relevant instruments. We use
$\widetilde{\bm{\mathbb{Z}}}_{\check{\mathcal{S}}_{1}}^{(1)}$ to represent the
design matrix formed by columns of $\widetilde{\bm{\mathbb{Z}}}^{(1)}$ indexed
by $\check{\mathcal{S}}_{1}$. Then we fit the reduced form models (2.3) and
(2.4) on the first subsample $H_{1}$ with design matrix
$\widetilde{\bm{\mathbb{Z}}}_{\check{\mathcal{S}}_{1}}^{(1)}$, and obtain de-
biased lasso estimators
$\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{1}}^{(1)}$ and
$\widehat{\bm{\Gamma}}_{\check{\mathcal{S}}_{1}}^{(1)}$. We also calculate the
standard error $\text{SE}(\widehat{\gamma}_{\check{\mathcal{S}}_{1},l}^{(1)})$
for $1\leq l\leq s$. We refer to (3.1) and (3.2) in the main paper for more
details on the de-biased lasso estimators and their estimated standard errors.
To further estimate the set of relevant instruments that are conditionally
associated with the exposure, we implement the joint thresholding and let
$\check{\mathcal{S}}_{2}=\\{l:1\leq l\leq
s,|\widehat{\gamma}_{\check{\mathcal{S}}_{1},l}^{(1)}|\geq\delta_{n}\times\text{SE}(\widehat{\gamma}_{\check{\mathcal{S}}_{1},l}^{(1)})\\}$
denote the set of candidates passing the joint thresholding, where
$\delta_{n}=\sqrt{\omega\log(\max\\{n_{1},s\\})}$. Next, locate spurious
instruments by estimating the range of causal effect estimates from pseudo
variables in $\check{\mathcal{S}}_{2}$ and remove spurious IVs from the
candidate set together with pseudo variables. The remaining set is represented
by $\check{\mathcal{S}}_{3}$ with cardinality $q$ and let
$\bm{Z}_{\check{\mathcal{S}}_{3}}\in\mathbb{R}^{q}$ denote the associated
candidate instruments. Further assume the plurality rule holds in
$\check{\mathcal{S}}_{3}$, i.e. the largest group of invalid or remaining
spurious instruments with the same Wald ratio is smaller than the group of
valid instruments. The procedure implemented on $H_{1}$ is outlined as Steps
00–3 in Algorithm 1.

Henceforth, we focus on the second subsample $H_{2}$. We have shown in Lemma
D.1 that relevant variables can pass joint thresholding with high probability.
By assuming that $\check{\mathcal{S}}_{3}$ includes all relevant variables,
the models (2.1) and (2.2) hold with $\bm{Z}$ replaced by
$\bm{Z}_{\check{\mathcal{S}}_{3}}$, and parameters
$\bm{\gamma}^{*},\bm{\pi}^{*}$ substituted by
$\bm{\gamma}_{\check{\mathcal{S}}_{3}}^{*},\bm{\pi}_{\check{\mathcal{S}}_{3}}^{*}\in\mathbb{R}^{q}$.
Similarly, the reduced form parameter $\bm{\Gamma}^{*}$ would become
$\bm{\Gamma}_{\check{\mathcal{S}}_{3}}^{*}=\bm{\pi}_{\check{\mathcal{S}}_{3}}^{*}+\beta^{*}\bm{\gamma}_{\check{\mathcal{S}}_{3}}^{*}$.
We use
$\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}\in\mathbb{R}^{n_{2}\times q}$
to denote the design matrix associated with $\check{\mathcal{S}}_{3}$ on
$H_{2}$. As $q$ is much smaller than $n_{2}$, we apply OLS to
$\bm{\mathrm{D}}^{(2)},\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}$ and
$\bm{\mathrm{Y}}^{(2)},\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}$
respectively in order to obtain estimates
$\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3}}^{(2)},\widehat{\bm{\Gamma}}_{\check{\mathcal{S}}_{3}}^{(2)}\in\mathbb{R}^{q}$,
and calculate the standard error of
$\widehat{\gamma}_{\check{\mathcal{S}}_{3},l}^{(2)}$ for $1\leq l\leq q$,

$\text{SE}\big{(}\widehat{\gamma}_{\check{\mathcal{S}}_{3},l}^{(2)}\big{)}=\sqrt{\frac{\Big{[}\Big{(}\bm{\mathbb{Z}}^{(2)\mathrm{\scriptstyle
T}}_{\check{\mathcal{S}}_{3}}\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}\Big{)}^{-1}\Big{]}_{ll}\big{\|}\bm{\mathrm{D}}^{(2)}-\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3}}^{(2)}\big{\|}^{2}}{n_{2}}}.$

Then we apply another joint thresholding by comparing
$|\widehat{\gamma}_{\check{\mathcal{S}}_{3},l}^{(2)}|$ with
$\delta_{n}\times\text{SE}(\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},l}^{(2)})$,
presented as Step 4 in Algorithm 1. Here $\delta_{n}=\sqrt{\omega\log n_{2}}$.
Let $\widetilde{\mathcal{S}}_{3}\subset\\{1,\ldots,q\\}$ be the index set of
candidates passing this joint thresholding on the second portion of data.

In Step 5 of Algorithm 1, we estimate valid IVs via mode finding. For
$j\in\widetilde{\mathcal{S}}_{3}$, calculate causal estimate
$\widehat{\beta}^{(2,j)}=\widehat{\Gamma}_{\check{\mathcal{S}}_{3},j}^{(2)}/\widehat{\gamma}_{\check{\mathcal{S}}_{3},j}^{(2)},$
Take $\widehat{\beta}^{(2,j)}$ as the underlying truth and compute the bias
for estimate $\widehat{\beta}^{(2,l)}$ for $l\in\widetilde{\mathcal{S}}_{3}$,
denoted by
$\widehat{b}^{(2,j,l)}=\widehat{\beta}^{(2,l)}-\widehat{\beta}^{(2,j)}$, and
the standard error $\text{SE}(\widehat{b}^{(2,j,l)})$. We discuss details of
$\text{SE}(\widehat{b}^{(2,j,l)})$ in Section E. We say two estimates
$\widehat{\beta}^{(2,l)}$ and $\widehat{\beta}^{(2,j)}$ are close if
$\big{|}\widehat{b}^{(2,j,l)}\big{|}\leq\text{SE}\big{(}\widehat{b}^{(2,j,l)}\big{)}\sqrt{\omega^{2}\log
n_{2}}.$ Here we use $\omega^{2}$ instead $\omega$ since we are conducting
pairwise tests. Note that $\widehat{b}^{(2,j,l)}$ remains the same if we swap
$j,l$, which implies that if $|\widehat{b}^{(2,j,l)}|$ pass the threshold,
then the candidate $j$ would be also in favor of $l$ by symmetry. The mode is
then estimated to be the effect estimate $\widehat{\beta}^{(2,l)}$ that is
close to the other estimates most frequently. Candidate instruments in
$\widetilde{\mathcal{S}}_{3}$ whose corresponding causal effect estimates are
close to the mode are identified as valid IVs, denoted by
$\check{\mathcal{S}}_{4}$.

###  A.2 Estimation of causal effect and statistical inference

We borrow the idea from Guo et al. (2018) for the construction of causal
effect estimate and its associated confidence interval. On the second
subsample $H_{2}$, we estimate the causal effect as follows

$\widehat{\beta}_{\text{split}}=\frac{\big{\\{}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}\big{\\}}^{\mathrm{\scriptstyle
T}}\widehat{\bm{W}}\widehat{\bm{\Gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}}{\big{\\{}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}\big{\\}}^{\mathrm{\scriptstyle
T}}\widehat{\bm{W}}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}},$
(A.1)

where
$\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}$
is a subvector of $\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3}}^{(2)}$
indexed by $\check{\mathcal{S}}_{4}$ and similarly for
$\widehat{\bm{\Gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}$.
The weighting matrix

$\widehat{\bm{W}}=[\widehat{\bm{\Sigma}}_{\check{\mathcal{S}}_{3}}^{(2)}]_{\check{\mathcal{S}}_{4},\check{\mathcal{S}}_{4}}-[\widehat{\bm{\Sigma}}_{\check{\mathcal{S}}_{3}}^{(2)}]_{\check{\mathcal{S}}_{4},\check{\mathcal{S}}_{4}^{c}}\Big{\\{}[\widehat{\bm{\Sigma}}_{\check{\mathcal{S}}_{3}}^{(2)}]_{\check{\mathcal{S}}_{4}^{c},\check{\mathcal{S}}_{4}^{c}}\Big{\\}}^{-1}[\widehat{\bm{\Sigma}}_{\check{\mathcal{S}}_{3}}^{(2)}]_{\check{\mathcal{S}}_{4}^{c},\check{\mathcal{S}}_{4}}$
(A.2)

with $\widehat{\bm{\Sigma}}_{\check{\mathcal{S}}_{3}}^{(2)}$ representing the
sample covariance matrix of candidate instruments corresponding to
$\check{\mathcal{S}}_{3}$ on the second subsample $H_{2}$, i.e.,
$\widehat{\bm{\Sigma}}_{\check{\mathcal{S}}_{3}}^{(2)}=\big{\\{}\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}\big{\\}}^{\mathrm{\scriptstyle
T}}\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}/n_{2}$, and
$[\widehat{\bm{\Sigma}}_{\check{\mathcal{S}}_{3}}^{(2)}]_{\mathcal{I}_{1},\mathcal{I}_{2}}$
is formed by rows and columns of
$\widehat{\bm{\Sigma}}_{\check{\mathcal{S}}_{3}}^{(2)}$ indexed by
$\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ respectively.

The asymptotic variance of $\widehat{\beta}_{\text{split}}$ is given by
$(\bm{\alpha}_{Y}^{*\mathrm{\scriptstyle
T}}\Sigma_{u}\bm{\alpha}_{Y}^{*}+\sigma_{Y}^{2})/v_{\text{split}}$ with

$v_{\text{split}}=\bm{\gamma}_{\check{\mathcal{S}}_{3},\mathcal{T}_{*}}^{*\mathrm{\scriptstyle
T}}\Big{(}[\bm{\Sigma}_{\check{\mathcal{S}}_{3}}]_{\mathcal{T}_{*},\mathcal{T}_{*}}-[\bm{\Sigma}_{\check{\mathcal{S}}_{3}}]_{\mathcal{T}_{*},\mathcal{T}_{*}^{c}}\Big{\\{}[\bm{\Sigma}_{\check{\mathcal{S}}_{3}}]_{\mathcal{T}_{*}^{c},\mathcal{T}_{*}^{c}}\Big{\\}}^{-1}[\bm{\Sigma}_{\check{\mathcal{S}}_{3}}]_{\mathcal{T}_{*}^{c},\mathcal{T}_{*}}\Big{)}\bm{\gamma}_{\check{\mathcal{S}}_{3},\mathcal{T}_{*}}^{*},$
(A.3)

where $\mathcal{T}_{*}=\\{l:\gamma_{\check{\mathcal{S}}_{3},l}^{*}\neq
0,\text{ for }1\leq l\leq q\\}\subset\\{1,\ldots,q\\}$ is the index set of
relevant instruments defined on $\check{\mathcal{S}}_{3}$, and
$v_{\text{split}}$ is estimated by

$\widehat{v}_{\text{split}}=\widehat{\mathrm{var}}(\widehat{\beta}_{\text{split}})=\frac{\widehat{\Theta}_{11}+\widehat{\beta}_{\text{split}}^{2}\widehat{\Theta}_{22}-2\widehat{\beta}_{\text{split}}\widehat{\Theta}_{12}}{\big{\\{}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}\big{\\}}^{\mathrm{\scriptstyle
T}}\widehat{\bm{W}}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}},$
(A.4)

where

$\displaystyle\widehat{\Theta}_{11}=\frac{1}{n_{2}}\Big{\|}\bm{\mathrm{Y}}^{(2)}-\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}\widehat{\bm{\Gamma}}_{\check{\mathcal{S}}_{3}}^{(2)}\Big{\|}_{2}^{2},$
$\displaystyle\widehat{\Theta}_{22}=\frac{1}{n_{2}}\Big{\|}\bm{\mathrm{D}}^{(2)}-\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3}}^{(2)}\Big{\|}_{2}^{2},$
$\displaystyle\widehat{\Theta}_{12}=\frac{1}{n_{2}}\Big{(}\bm{\mathrm{Y}}^{(2)}-\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}\widehat{\bm{\Gamma}}_{\check{\mathcal{S}}_{3}}^{(2)}\Big{)}^{\mathrm{\scriptstyle
T}}\Big{(}\bm{\mathrm{D}}^{(2)}-\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3}}^{(2)}\Big{)}.$
(A.5)

We then can further obtain the $100(1-\alpha)\%$ confidence interval for
$\beta$:

$\Big{[}\widehat{\beta}_{\text{split}}-z_{1-\alpha/2}\sqrt{\widehat{v}_{\text{split}}/n_{2}},\widehat{\beta}_{\text{split}}+z_{1-\alpha/2}\sqrt{\widehat{v}_{\text{split}}/n_{2}}\Big{]},$
(A.6)

where $z_{1-\alpha/2}$ is the $(1-\alpha/2)$-quantile of the standard normal
distribution. See Steps 6 and 7 of Algorithm 1 for an outline of estimation of
causal effect and construction of confidence interval.

INPUT: Design matrix $\bm{\mathbb{Z}}\in\mathbb{R}^{n\times p}$, observed
exposure $\bm{\mathrm{D}}\in\mathbb{R}^{n}$, observed outcome
$\bm{\mathrm{Y}}\in\mathbb{R}^{n}$.

  * 00.

Divide data into two halves:

$H_{1}$ of size $n_{1}$ and $H_{2}$ of size $n_{2}$;

Denote the associated data as $\bm{\mathbb{Z}}^{(1)}\in\mathbb{R}^{n_{1}\times
p},\bm{\mathrm{D}}^{(1)}\in\mathbb{R}^{n_{1}},\bm{\mathrm{Y}}^{(1)}\in\mathbb{R}^{n_{1}}$,
and $\bm{\mathbb{Z}}^{(2)}\in\mathbb{R}^{n_{2}\times
p},\bm{\mathrm{D}}^{(2)}\in\mathbb{R}^{n_{2}},\bm{\mathrm{Y}}^{(2)}\in\mathbb{R}^{n_{2}}$.

  * 0.

Pseudo variables generation on $H_{1}$:

Generate $p$ pseudo variables $Z_{p+1}^{(1)},\ldots,Z_{2p}^{(1)}$ by randomly
permuting the rows of $\bm{\mathbb{Z}}^{(1)}$;

Denote the resulting design matrix as ${\bar{\bm{\mathbb{Z}}}^{(1)}}$, and the
concatenation of $\bm{\mathbb{Z}}^{(1)}$ and ${\bar{\bm{\mathbb{Z}}}^{(1)}}$
as $\widetilde{\bm{\mathbb{Z}}}^{(1)}$.

  * 1.

Marginal screening on $H_{1}$:

For $j=1,\ldots,2p$, compute
$\widehat{\rho}_{j}=|\sum_{i=1}^{n_{1}}\widetilde{Z}_{ij}^{(1)}D_{i}^{(1)}|/\sum_{i=1}^{n_{1}}\Big{(}\widetilde{Z}_{ij}^{(1)}\Big{)}^{2}$;

Select the $s$ variables with the largest $\widehat{\rho}_{j}$ and denote them
as $\check{\mathcal{S}}_{1}$.

  * 2.

Joint thresholding on $H_{1}$:

Fit a linear model for
$\bm{\mathrm{D}}^{(1)}\sim\widetilde{\bm{\mathbb{Z}}}_{\check{\mathcal{S}}_{1}}^{(1)}$:
let $\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{1}}^{(1)}$ be the debiased
lasso estimate constructed similarly to (3.1) in the main paper; here
$\widetilde{\bm{\mathbb{Z}}}_{\check{\mathcal{S}}_{1}}^{(1)}$ is formed by
columns of $\widetilde{\bm{\mathbb{Z}}}^{(1)}$ indexed by
$\check{\mathcal{S}}_{1}$. Also obtain their standard error estimate
$\text{SE}(\widehat{\gamma}_{\check{\mathcal{S}}_{1},l}^{(1)})$.

Fit a linear model for
$\bm{\mathrm{Y}}^{(1)}\sim\widetilde{\bm{\mathbb{Z}}}_{\check{\mathcal{S}}_{1}}^{(1)}$
and obtain the estimate
$\widehat{\bm{\Gamma}}_{\check{\mathcal{S}}_{1}}^{(1)}$ in a similar fashion.

Estimate the relevant variables
$\check{\mathcal{S}}_{2}=\big{\\{}l\in[s]:|\widehat{\gamma}_{\check{\mathcal{S}}_{1},l}^{(1)}|\geq\delta_{n}\times\text{SE}(\widehat{\gamma}_{\check{\mathcal{S}}_{1},l}^{(1)})\big{\\}},$
where $\delta_{n}=\sqrt{\omega\log(\max\\{n_{1},s\\})}$, where $\omega$ is a
tuning parameter. Denote pseudos in $\check{\mathcal{S}}_{2}$ by
$\widehat{\mathcal{P}}$.

  * 3.

Irrelevant variables removal:

For $l\in\check{\mathcal{S}}_{2}$, calculate causal effect estimate
$\widehat{\beta}^{(1,l)}={\widehat{\Gamma}_{\check{\mathcal{S}}_{1},l}^{(1)}}/{\widehat{\gamma}_{\check{\mathcal{S}}_{1},l}^{(1)}};$

Remove variables similar to pseudo variables in their causal effect estimate
$\widehat{\beta}^{(1,l)}$ from $\check{\mathcal{S}}_{2}$:

$\check{\mathcal{S}}_{3}=\check{\mathcal{S}}_{2}\setminus\left\\{l\in\check{\mathcal{S}}_{2}:\min_{j\in\widehat{\mathcal{P}}}\widehat{\beta}^{(1,j)}\leq\widehat{\beta}^{(1,l)}\leq\max_{j\in\widehat{\mathcal{P}}}\widehat{\beta}^{(1,j)}\right\\}.$

Denote the cardinality of $\check{\mathcal{S}}_{3}$ as $q$.

Algorithm 1 Fighting noise with noise – sample splitting

  * 4.

Joint thresholding on $H_{2}$:

Let $\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}\in\mathbb{R}^{n_{2}\times
q}$ denote the associated design matrix of $\check{\mathcal{S}}_{3}$ on
$H_{2}$. Fit a linear model for
$\bm{\mathrm{D}}^{(2)}\sim\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}$:
obtain the OLS estimator
$\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3}}^{(2)}$ and standard error
estimate $\text{SE}(\widehat{\gamma}_{\check{\mathcal{S}}_{3},l}^{(2)})$ for
$1\leq l\leq q$;

Fit a linear model for
$\bm{\mathrm{Y}}^{(2)}\sim\bm{\mathbb{Z}}^{(2)}_{\check{\mathcal{S}}_{3}}$ and
obtain the estimate $\widehat{\bm{\Gamma}}_{\check{\mathcal{S}}_{3}}^{(2)}$ in
a similar fashion;

Further approximate the set of relevant instruments as
$\widetilde{\mathcal{S}}_{3}=\big{\\{}l\in[q]:|\widehat{\gamma}_{\check{\mathcal{S}}_{3},l}^{(2)}|\geq\delta_{n}\times\text{SE}(\widehat{\gamma}_{\check{\mathcal{S}}_{3},l}^{(2)})\big{\\}},$
where $\delta_{n}=\sqrt{\omega\log n_{2}}$.

  * 5.

Mode-finding on $H_{2}$:

For $l\in\widetilde{\mathcal{S}}_{3}$, calculate causal effect estimate
$\widehat{\beta}^{(2,l)}=\widehat{\Gamma}_{\check{\mathcal{S}}_{3},l}^{(2)}/\widehat{\gamma}_{\check{\mathcal{S}}_{3},l}^{(2)};$

For $j,l\in\widetilde{\mathcal{S}}_{3}$, compute
$\widehat{b}^{(2,j,l)}=\widehat{\beta}^{(2,l)}-\widehat{\beta}^{(2,j)}$ and
$\text{SE}(\widehat{b}^{(2,j,l)})$;

For $j\in\widetilde{\mathcal{S}}_{3},$ find the number of candidates leading
to a similar causal effect estimate with $j$:

$V_{j}=\Big{\lvert}\Big{\\{}l\in\widetilde{\mathcal{S}}_{3}:|\widehat{b}^{(2,j,l)}|/\text{SE}(\widehat{b}^{(2,j,l)})\leq\sqrt{\omega^{2}\log
n_{2}}\Big{\\}}\Big{\rvert};$

Find valid IVs using the plurality rule:

$\check{\mathcal{S}}_{4}=\Big{\\{}j:V_{j}=\max_{l\in\widetilde{\mathcal{S}}_{3}}V_{l}\Big{\\}}.$

  * 6.

Causal effect estimation:

$\widehat{\beta}_{\text{split}}=\frac{\big{\\{}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}\big{\\}}^{\mathrm{\scriptstyle
T}}\widehat{\bm{W}}\widehat{\bm{\Gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}}{\big{\\{}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}\big{\\}}^{\mathrm{\scriptstyle
T}}\widehat{\bm{W}}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}},$

where
$\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}$
is a subvector of $\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3}}^{(2)}$
indexed by $\check{\mathcal{S}}_{4}$ and similarly for
$\widehat{\bm{\Gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}$.
See (A.2) for the explicit form of the weighting matrix $\widehat{\bm{W}}$.

Estimated variance of $\widehat{\beta}_{\text{split}}$:

$\widehat{v}_{\text{split}}=\widehat{\mathrm{var}}(\widehat{\beta}_{\text{split}})=\frac{\widehat{\Theta}_{11}+\widehat{\beta}_{\text{split}}^{2}\widehat{\Theta}_{22}-2\widehat{\beta}_{\text{split}}\widehat{\Theta}_{12}}{\big{\\{}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}\big{\\}}^{\mathrm{\scriptstyle
T}}\widehat{\bm{W}}\widehat{\bm{\gamma}}_{\check{\mathcal{S}}_{3},\check{\mathcal{S}}_{4}}^{(2)}},$

where $\widehat{\Theta}_{11},\widehat{\Theta}_{22},\widehat{\Theta}_{12}$ are
given in (A.2).

  * 7.

$100(1-\alpha)\%$ confidence interval:
$[\widehat{\beta}_{\text{split}}-z_{1-\alpha/2}\sqrt{\widehat{v}_{\text{split}}/n_{2}},\widehat{\beta}_{\text{split}}+z_{1-\alpha/2}\sqrt{\widehat{v}_{\text{split}}/n_{2}}].$

Algorithm S1 (continue)

##  Appendix B More Details on Data Application

###  B.1 Relationship between BMI and HUI-3

We display the observed relationship between BMI and Health Utility Index Mark
3 (HUI-3) in Figure 6.

(a) BMI vs. HUI-3 (b) BMI (overweight) vs. HUI-3

Figure 6: The plot on the left panel depicts the relationship between BMI and
Health Utility Index Mark 3 (HUI-3) for the WLS data used in our analysis. On
the right panel, we zoom in on BMI above 25 and HUI-3 accordingly.

We can learn from Figure 6-(a) that HUI-3 has a positive association with BMI
when BMI is lower than 25, and is negatively associated with BMI when BMI is
higher than 25. This is consistent with the fact that both underweight and
overweight are linked with a range of adverse health conditions. Meanwhile,
the healthy weight range is between 18.5 and 24.9. We provide additional plot
in Figure 6-(b) to emphasize that as BMI increases over 25, the risk for lower
quality of life increases, which is our focus in the data application.

###  B.2 WLS Data Pre-processing

For genomic data, we use imputed genetic variants extracted from whole-genome
sequencing on 22 chromosomes, in which unobserved genotypes in the study
sample are estimated from the haplotype or genotype reference panel (Browning,
2008; Li et al., 2009). This allows us to examine the evidence for association
at biomarkers that are not directly genotyped (Scheet & Stephens, 2006;
Marchini et al., 2007; Browning, 2008). The input of missing data imputation
are observed variants that meet criterion of quality control: for example,
missing call rate $<$ 2%, minor allele frequency (MAF) $\geq$ 0.01, Hardy-
Weinberg Equilibrium (HWE) p-value $\geq 10^{-4}$, with known and unique
positions, etc. See the genotype quality control report attached to Herd et
al. (2014), which follows the quality filters described in Laurie et al.
(2010). Genotype imputation were performed using IMPUTE2 software (Howie et
al., 2009, 2011). To account for uncertainty of the imputation output, we
exclude variants with BEAGLE allelic $r^{2}$ less than 0.3, which is an
assessment of genotype imputation performance (de Bakker et al., 2008; Li et
al., 2010).

For the remaining missing values in genotype data after matching with
available phenotypes, we replace them by random sampling from the empirical
allele distribution. Namely, the missing values at a genetic variant $j$ are
sampled from the associated marginal allele empirical distribution. In
addition, there might be identical columns in our current SNP matrix, in which
rows are indexed by participants considered in our analysis and columns are
indexed by genotypes. To proceed, we keep unique ones in our analysis and
recheck at the end whether any excluded variants have exactly the same (0,1,2)
genotype representation as our estimated valid IVs. If there are any, they
will also be identified as valid instruments.

##  Appendix C Notations

For any vector $\bm{v}\in\mathbb{R}^{p}$, $v_{j}$ represents the $j$th element
of $\bm{v}$. $\|\bm{v}\|_{1}$, $\|\bm{v}\|_{2}$ and $\|\bm{v}\|_{\infty}$
represent the usual 1-, 2- and $\infty$\- norms respectively. The $\ell_{0}$
norm $\|\bm{v}\|_{0}$ refers to the number of nonzero components in $\bm{v}$.
For a set $S$, $|S|$ denotes its cardinality and $S^{c}$ is the complement of
$S$. If $S\subset\\{1,\ldots,p\\}$ is an index set, then
$\bm{v}_{S}\in\mathbb{R}^{|S|}$ is a subvector of $\bm{v}\in\mathbb{R}^{p}$
indexed by $S$ and $v_{S,l}$ denotes the $l$th element of $\bm{v}_{S}$ for
$1\leq l\leq|S|$. For an $n\times p$ matrix
$\bm{\mathrm{A}}\in\mathbb{R}^{n\times p}$, we denote the $(i,j)$th element of
$\bm{\mathrm{A}}$ by $\bm{\mathrm{A}}_{ij}$, the $i$th row as
$\bm{\mathrm{A}}_{i.}$ and the $j$th column as $\bm{\mathrm{A}}_{.j}$.
$\bm{\mathrm{A}}^{\mathrm{\scriptstyle T}}$ denotes the transpose of
$\bm{\mathrm{A}}$. $\|\bm{A}\|_{\infty}=\max_{1\leq i,j\leq p}|A_{ij}|$
represents the maximum of the absolute element of matrix $\bm{A}$ and
$\|\bm{A}\|_{1}=\max_{1\leq j\leq n}\sum_{i=1}^{p}|A_{ij}|$ is the maximum
absolute column sum of $\bm{\mathrm{A}}$. For a symmetric matrix $\bm{A}$, let
$\lambda_{\min}(\bm{A})$ and $\lambda_{\max}(\bm{A})$ denote the smallest and
the largest eigenvalues respectively. The spectral norm
$\|\bm{A}\|_{2}=\lambda_{\max}(\bm{A})$. For functions $f(n)$ and $h(n)$, we
write $f(n)=O(h(n))$ if $\lim_{n\rightarrow\infty}f(n)/h(n)=c$, where $c$ is
some constant. Particularly, if $h(n)=c$, the notation becomes $f(n)=O(1)$.
Similarly, if $\lim_{n\rightarrow\infty}f(n)/h(n)=0$, we write $f(n)=o(h(n))$.
Specifically, if $h(n)=c$, we write $f(n)=o(1)$. The notation $a_{n}\lesssim
c_{n}$ represents that there exists a constant $C$ such that $a_{n}\leq
Cb_{n}$ and $a_{n}\gtrsim b_{n}$ means that there exists a constant $c$ such
that $a_{n}\geq cb_{n}$. And $a_{n}\asymp d_{n}$ represents
$a_{n}=\widetilde{c}b_{n}$ for some constant $\widetilde{c}$.
$\text{cor}(X,Y)$ represents correlation between random variables $X$ and $Y$,
and let $\mathrm{cov}(X,Y)$ denote the associated covariance. If
$\bm{X}\in\mathbb{R}^{p}$ and $\bm{Y}\in\mathbb{R}^{q}$ are random vectors,
then $\mathrm{cov}(\bm{X},\bm{Y})\in\mathbb{R}^{p\times q}$ is the cross-
covariance matrix.

##  Appendix D Technical Details

###  D.1 Setup for Technical Details

We first revisit some definitions and notations. The index set of candidate
instruments passing marginal screening is
$\widehat{\mathcal{S}}_{1}\subset\\{1,\ldots,p\\}$ that has size $s$. Let
$\bm{Z}_{\widehat{\mathcal{S}}_{1}}$ denote the subvector of
$\bm{Z}\in\mathbb{R}^{p}$ indexed by $\widehat{\mathcal{S}}_{1}$. Consider the
models (2.1) and (2.2) with $\bm{Z}\in\mathbb{R}^{p}$ replaced by
$\bm{Z}_{\widehat{\mathcal{S}}_{1}}\in\mathbb{R}^{s}$. The parameters
$\bm{\gamma}^{*},\bm{\Gamma}^{*}\in\mathbb{R}^{p}$ change accordingly to
$\bm{\gamma}_{\widehat{\mathcal{S}}_{1}}^{*},\bm{\Gamma}_{\widehat{\mathcal{S}}_{1}}^{*}\in\mathbb{R}^{s}$.
The debiased-lasso estimators of
$\bm{\gamma}_{\widehat{\mathcal{S}}_{1}}^{*},\bm{\Gamma}_{\widehat{\mathcal{S}}_{1}}^{*}$
under models (2.1) and (2.2) are denoted by
$\widehat{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}},\widehat{\bm{\Gamma}}_{\widehat{\mathcal{S}}_{1}}\in\mathbb{R}^{s}$
respectively, and
$\widetilde{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}},\widetilde{\bm{\Gamma}}_{\widehat{\mathcal{S}}_{1}}\in\mathbb{R}^{s}$
are the lasso estimators defined after (3.1). The matrix
$\bm{M}\in\mathbb{R}^{s\times s}$ is given in (G.3) of the main paper and
$\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}\in\mathbb{R}^{n\times s}$ is
formed by columns of $\bm{\mathbb{Z}}\in\mathbb{R}^{n\times p}$ indexed by
$\widehat{\mathcal{S}}_{1}$.

It follows from (4) and (5) of van de Geer et al. (2014) that the debiased-
lasso estimators satisfy the following decomposition

$\displaystyle\sqrt{n}(\widehat{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}}-\bm{\gamma}_{\widehat{\mathcal{S}}_{1}}^{*})=\underbrace{\sqrt{n}(\bm{M}\widehat{\bm{\Sigma}}_{\widehat{\mathcal{S}}_{1}}-\bm{\mathrm{I}})(\bm{\gamma}_{\widehat{\mathcal{S}}_{1}}^{*}-\widetilde{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}})}_{\Delta_{1}}+\frac{1}{\sqrt{n}}\bm{M}\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{U}}\bm{\alpha}_{D}^{*}+\bm{\epsilon}_{D})$ (D.1)
$\displaystyle\sqrt{n}(\widehat{\bm{\Gamma}}_{\widehat{\mathcal{S}}_{1}}-\bm{\Gamma}_{\widehat{\mathcal{S}}_{1}}^{*})=\underbrace{\sqrt{n}(\bm{M}\widehat{\bm{\Sigma}}_{\widehat{\mathcal{S}}_{1}}-\bm{\mathrm{I}})(\bm{\Gamma}_{\widehat{\mathcal{S}}_{1}}^{*}-\widetilde{\bm{\Gamma}}_{\widehat{\mathcal{S}}_{1}})}_{\Delta_{2}}+\frac{1}{\sqrt{n}}\bm{M}\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}^{\mathrm{\scriptstyle
T}}\\{\bm{\mathrm{U}}(\bm{\alpha}_{Y}^{*}+\beta^{*}\bm{\alpha}_{D}^{*})+(\beta^{*}\bm{\epsilon}_{D}+\bm{\epsilon}_{Y})\\},$
(D.2)

where $\bm{\mathrm{U}}\in\mathbb{R}^{n\times g}$,
$\widehat{\bm{\Sigma}}_{\widehat{\mathcal{S}}_{1}}=\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}^{\mathrm{\scriptstyle
T}}\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}/n$, and
$\bm{\mathrm{I}}\in\mathbb{R}^{s\times s}$ is an identity matrix. Define
$\bm{v}_{l}(\widehat{\mathcal{S}}_{1})=\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}\bm{M}^{\mathrm{\scriptstyle
T}}e_{l}/n\in\mathbb{R}^{n}$ for $1\leq l\leq s$, where
$e_{l}\in\mathbb{R}^{s}$ is a unit vector with the $l$th component equal to
one and zero otherwise. Hereafter, we use $\bm{v}_{l}$ to refer to
$\bm{v}_{l}(\widehat{\mathcal{S}}_{1})$ for ease of notation. The $\ell_{2}$
norm of $\bm{v}_{l}$ is
$\|\bm{v}_{l}\|^{2}=[\bm{M}\widehat{\bm{\Sigma}}_{\widehat{\mathcal{S}}_{1}}\bm{M}^{\mathrm{\scriptstyle
T}}]_{ll}/n$. Here $[\cdot]_{ll}$ denotes the $l$th diagonal element of the
matrix.

Let $\eta_{max}=1\vee
B\sup_{\widehat{\mathcal{S}}_{1}:|\widehat{\mathcal{S}}_{1}|=s}\max_{1\leq
j\leq s}\|[\bm{\Sigma}_{\widehat{\mathcal{S}}_{1}}^{-1}]_{-j,j}\|_{1}$. Next,
we define the events, on which we establish our theoretical results.

$\displaystyle\mathcal{A}=\big{\\{}\mathcal{S}_{*}\subset\widehat{\mathcal{S}}_{1}\big{\\}}$
(D.3)
$\displaystyle\mathcal{B}_{1}=\bigg{\\{}\sup_{\widehat{\mathcal{S}}_{1}:|\widehat{\mathcal{S}}_{1}|=s}\|\widetilde{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}}-\bm{\gamma}_{\widehat{\mathcal{S}}_{1}}^{*}\|_{1}\leq\frac{16s_{1}^{*}\lambda_{\gamma}}{C_{L}}\bigg{\\}}$
(D.4)
$\displaystyle\mathcal{B}_{2}=\bigg{\\{}\sup_{\widehat{\mathcal{S}}_{1}:|\widehat{\mathcal{S}}_{1}|=s}\|\widetilde{\bm{\Gamma}}_{\widehat{\mathcal{S}}_{1}}-\bm{\Gamma}_{\widehat{\mathcal{S}}_{1}}^{*}\|_{1}\leq\frac{16s_{2}^{*}\lambda_{\Gamma}}{C_{L}}\bigg{\\}}$
(D.5)
$\displaystyle\mathcal{B}_{3}=\bigg{\\{}\sup_{\widehat{\mathcal{S}}_{1}:|\widehat{\mathcal{S}}_{1}|=s}\frac{\|\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}(\bm{\gamma}_{\widehat{\mathcal{S}}_{1}}^{*}-\widetilde{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}})\|^{2}}{n}\leq\frac{16s_{1}^{*}\lambda_{\gamma}^{2}}{C_{L}}\bigg{\\}}$
(D.6)
$\displaystyle\widetilde{\mathcal{B}}=\bigg{\\{}\sup_{\widehat{\mathcal{S}}_{1}:|\widehat{\mathcal{S}}_{1}|=s}\bigg{|}\frac{1}{n}\|\bm{\mathrm{D}}-\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}\widetilde{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}}\|^{2}-(\bm{\alpha}_{D}^{*\mathrm{\scriptstyle
T}}\Sigma_{U}\bm{\alpha}_{D}^{*}+\sigma_{D}^{2})\bigg{|}\leq
a_{4}\sqrt{\frac{\log
n}{n}}+\frac{16\lambda_{\gamma}^{2}s_{1}^{*}}{C_{L}}\bigg{\\}}$ (D.7)
$\displaystyle\mathcal{C}_{1}=\bigg{\\{}\sup_{\widehat{\mathcal{S}}_{1}:|\widehat{\mathcal{S}}_{1}|=s}\|\Delta_{1}\|_{\infty}\leq\frac{70a_{2}s_{1}^{*}\sqrt{\log
p}\lambda_{\gamma}}{C_{L}^{2}}\eta_{max}\bigg{\\}}$ (D.8)
$\displaystyle\mathcal{C}_{2}=\bigg{\\{}\sup_{\widehat{\mathcal{S}}_{1}:|\widehat{\mathcal{S}}_{1}|=s}\|\Delta_{2}\|_{\infty}\leq\frac{70a_{2}s_{2}^{*}\sqrt{\log
p}\lambda_{\Gamma}}{C_{L}^{2}}\eta_{max}\bigg{\\}}$ (D.9)
$\displaystyle\mathcal{C}_{3}=\bigg{\\{}\frac{1}{n}\max_{1\leq j\leq
p}|\bm{Z}_{j}^{\mathrm{\scriptstyle T}}\bm{\epsilon}|\leq
a_{5}\sqrt{\frac{\log p}{n}}\bigg{\\}}$ (D.10)
$\displaystyle\mathcal{D}=\bigg{\\{}\frac{1}{C_{R}}-\mu\leq
n\|\bm{v}_{l}\|^{2}\leq\frac{1}{C_{L}}+\mu,\forall\widehat{\mathcal{S}}_{1}\text{
with }|\widehat{\mathcal{S}}_{1}|=s\text{ and }\forall 1\leq l\leq
s\bigg{\\}}$ (D.11)
$\displaystyle\mathcal{E}_{1}=\bigg{\\{}\sup_{\widehat{\mathcal{S}}_{1}:|\widehat{\mathcal{S}}_{1}|=s}\|\bm{M}\widehat{\bm{\Sigma}}_{\widehat{\mathcal{S}}_{1}}-\bm{\mathrm{I}}_{s}\|_{\infty}\leq\frac{\lambda_{\text{NL}}}{0.99C_{L}}\bigg{\\}}$
(D.12)
$\displaystyle\mathcal{E}_{2}=\bigg{\\{}\sup_{\widehat{\mathcal{S}}_{1}:|\widehat{\mathcal{S}}_{1}|=s}\|\bm{M}\|_{1}\leq\frac{10\eta_{max}}{C_{L}}\bigg{\\}},$
(D.13)
$\displaystyle\mathcal{E}_{3}=\bigg{\\{}\sup_{\widehat{\mathcal{S}}_{1}:|\widehat{\mathcal{S}}_{1}|=s}\max_{1\leq
j\leq
s}|[\bm{M}\widehat{\bm{\Sigma}}_{\widehat{\mathcal{S}}_{1}}\bm{M}^{\mathrm{\scriptstyle
T}}-\bm{\Sigma}_{\widehat{\mathcal{S}}_{1}}^{-1}]_{jj}|\leq\mu\bigg{\\}},$
(D.14)

where $\mathcal{S}_{*}$ denotes the index set of relevant variables to the
exposure that is defined on $\\{1,\ldots,p\\}$; $s_{1}^{*},s_{2}^{*}$
represent the number of relevant variables to the exposure and the number of
invalid instruments, respectively; $C_{L},C_{R}$ lower and upper bound the
eigenvalues of $\bm{\Sigma}\in\mathbb{R}^{p\times p}$ and $B$ is an upper
bound placed on $[\bm{\Sigma}]_{jj}$, as shown in assumption A2;
$\lambda_{\gamma}=C_{\gamma}\sqrt{(\log p)/n}$ and
$\lambda_{\Gamma}=C_{\Gamma}\sqrt{(\log p)/n}$ are lasso tuning parameters,
where constants $C_{\gamma}\geq
5.7\sqrt{B(\bm{\alpha}_{D}^{*\mathrm{\scriptstyle
T}}\Sigma_{U}\bm{\alpha}_{D}^{*}+\sigma_{D}^{2})}$, and $C_{\Gamma}\geq
5.7\sqrt{B\\{(\beta^{*}\bm{\alpha}_{D}^{*}+\bm{\alpha}_{Y}^{*})^{\mathrm{\scriptstyle
T}}\Sigma_{U}(\beta^{*}\bm{\alpha}_{D}^{*}+\bm{\alpha}_{Y}^{*})+\beta^{*2}\sigma_{D}^{2}+\sigma_{Y}^{2}\\}}$;
$\lambda_{\text{NL}}=4a_{2}\eta_{\max}\sqrt{(\log p)/n}$;
$\bm{\Sigma}_{\widehat{\mathcal{S}}_{1}}^{-1}$ denotes the inverse of the
submatrix $\bm{\Sigma}_{\widehat{\mathcal{S}}_{1}}$ formed by deleting rows
and columns of $\bm{\Sigma}$ indexed by $\widehat{\mathcal{S}}_{1}$ and
$[]_{-j,j}$ refers to the $j$th column without the $j$th component of the
matrix; $a_{2},a_{4},a_{5}$ are positive constants and particularly,
$a_{2}>4\sqrt{3}eB\kappa^{2}$ where $\kappa$ is the sub-Gaussian norm of
$\bm{\Sigma}^{-1/2}\bm{Z}_{1}$ and is a finite, positive constant, and
$a_{5}>\sqrt{1.01B\mathrm{var}(\epsilon)}$; $\mu$ is defined as

$\mu=\zeta_{1}\bigg{(}\frac{\log p}{n}\bigg{)}^{1/4}\eta_{\max}^{2}$ (D.15)

for some constant $\zeta_{1}>0$.

On the event $\mathcal{A}$, the set of candidates passing marginal screening,
$\widehat{\mathcal{S}}_{1}$, includes all relevant variables to the exposure,
and the models (2.1) and (2.2) hold with $\bm{Z}\in\mathbb{R}^{p}$ replaced by
$\bm{Z}_{\widehat{\mathcal{S}}_{1}}\in\mathbb{R}^{s}$. The events
$\mathcal{B}_{1},\mathcal{B}_{2},\mathcal{B}_{3}$ provide convergence rates
for lasso post marginal screening estimators as well as the prediction error.
Control of the variance estimate for $\bm{U}^{\mathrm{\scriptstyle
T}}\bm{\alpha}_{D}^{*}+\epsilon_{D}$ using lasso post marginal screening
estimators is introduced in the event $\widetilde{\mathcal{B}}$. Events
$\mathcal{C}_{1},\mathcal{C}_{2}$ give the convergence rates for the error
terms $\Delta_{1},\Delta_{2}$ as shown in (D.1) and (D.2) uniformly over all
possible $\widehat{\mathcal{S}}_{1}$ of size $s$. The event $\mathcal{C}_{3}$
places an upper bound for the maximal correlation between candidate $Z_{j}$
and a random error term $\epsilon$ that is independent of $Z_{j}$ for all
$j=1,\ldots,p$, which is the key ingredient in controlling convergence rates
of the lasso type estimators. The error term $\epsilon$ will vary place to
place to accommodate with different bounds of interest. The bounded support
for $n\|\bm{v}_{l}\|^{2}$ is provided in the event $\mathcal{D}$, which will
be used to explore the order of the joint threshold.

Let $\mathcal{E}$ be the intersection of all events defined in (D.3)–(D.14).
Later we will show that there exists a constant $C_{1}>0$ such that when $n,p$
are sufficiently large, we have

$P(\mathcal{E})\geq 1-e^{-C_{1}\log(n\wedge p)}.$ (D.16)

###  D.2 Proof for Theorem 3.1 and Theorem 3.2

In this section, we give proof for the two main theorems. Before proceeding,
we first derive the bias of the ratio estimates associated with the candidate
instruments passing marginal screening. Under assumptions A3, on the event
$\mathcal{A}$, the difference between the ratio estimate for each candidate
instrument in $\widehat{\mathcal{S}}_{1}$ and the true causal effect
$\beta^{*}$, i.e.,
$\widehat{\Gamma}_{\widehat{\mathcal{S}}_{1},l}/\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}-\beta^{*}$
takes the following form

$\displaystyle\frac{\widehat{\Gamma}_{\widehat{\mathcal{S}}_{1},l}}{\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}}-\beta^{*}$
$\displaystyle=\frac{e_{l}^{\mathrm{\scriptstyle
T}}\sqrt{n}(\widehat{\bm{\Gamma}}_{\widehat{\mathcal{S}}_{1}}-\bm{\Gamma}_{\widehat{\mathcal{S}}_{1}}^{*})+\sqrt{n}\Gamma_{\widehat{\mathcal{S}}_{1},l}^{*}}{e_{l}^{\mathrm{\scriptstyle
T}}\sqrt{n}(\widehat{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}}-\bm{\gamma}_{\widehat{\mathcal{S}}_{1}}^{*})+\sqrt{n}\gamma_{\widehat{\mathcal{S}}_{1},l}^{*}}-\beta^{*}$
$\displaystyle=\frac{e_{l}^{\mathrm{\scriptstyle
T}}\Delta_{2}+e_{l}^{\mathrm{\scriptstyle
T}}\frac{1}{\sqrt{n}}\bm{M}\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}^{\mathrm{\scriptstyle
T}}\\{\beta^{*}(\bm{\mathrm{U}}\bm{\alpha}_{D}^{*}+\bm{\epsilon}_{D})+\bm{\mathrm{U}}\bm{\alpha}_{Y}^{*}+\epsilon_{Y}\\}+\sqrt{n}\Gamma_{\widehat{\mathcal{S}}_{1},l}^{*}}{e_{l}^{\mathrm{\scriptstyle
T}}\Delta_{1}+e_{l}^{\mathrm{\scriptstyle
T}}\frac{1}{\sqrt{n}}\bm{M}\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{U}}\bm{\alpha}_{D}^{*}+\bm{\epsilon}_{D})+\sqrt{n}\gamma_{\widehat{\mathcal{S}}_{1},l}^{*}}-\beta^{*}$
$\displaystyle=\frac{e_{l}^{\mathrm{\scriptstyle
T}}\Delta_{2}+e_{l}^{\mathrm{\scriptstyle
T}}\frac{1}{\sqrt{n}}\bm{M}\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}^{\mathrm{\scriptstyle
T}}\\{\beta^{*}(\bm{\mathrm{U}}\bm{\alpha}_{D}^{*}+\bm{\epsilon}_{D})+\bm{\mathrm{U}}\bm{\alpha}_{Y}^{*}+\epsilon_{Y}\\}+\sqrt{n}(\pi_{\widehat{\mathcal{S}}_{1},l}^{*}+\beta^{*}\gamma_{\widehat{\mathcal{S}}_{1},l}^{*})}{e_{l}^{\mathrm{\scriptstyle
T}}\Delta_{1}+e_{l}^{\mathrm{\scriptstyle
T}}\frac{1}{\sqrt{n}}\bm{M}\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{U}}\bm{\alpha}_{D}^{*}+\bm{\epsilon}_{D})+\sqrt{n}\gamma_{\widehat{\mathcal{S}}_{1},l}^{*}}-\beta^{*}$
$\displaystyle=\beta^{*}+\frac{e_{l}^{\mathrm{\scriptstyle
T}}(\Delta_{2}-\beta^{*}\Delta_{1})+e_{l}^{\mathrm{\scriptstyle
T}}\frac{1}{\sqrt{n}}\bm{M}\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{U}}\bm{\alpha}_{Y}^{*}+\bm{\epsilon}_{Y})+\sqrt{n}\pi_{\widehat{\mathcal{S}}_{1},l}^{*}}{e_{l}^{\mathrm{\scriptstyle
T}}\Delta_{1}+e_{l}^{\mathrm{\scriptstyle
T}}\frac{1}{\sqrt{n}}\bm{M}\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{U}}\bm{\alpha}_{D}^{*}+\bm{\epsilon}_{D})+\sqrt{n}\gamma_{\widehat{\mathcal{S}}_{1},l}^{*}}-\beta^{*}$
$\displaystyle=\frac{e_{l}^{\mathrm{\scriptstyle
T}}(\Delta_{2}-\beta^{*}\Delta_{1})+\sqrt{n}\bm{v}_{l}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{U}}\bm{\alpha}_{Y}^{*}+\bm{\epsilon}_{Y})+\sqrt{n}\pi_{\widehat{\mathcal{S}}_{1},l}^{*}}{\sqrt{n}\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}}$
(D.17)

for $1\leq l\leq s$, where the second and the last equality come from the
decomposition (D.1) and (D.2); the third equality follows from
$\bm{\Gamma}_{\widehat{\mathcal{S}}_{1}}^{*}=\bm{\pi}_{\widehat{\mathcal{S}}_{1}}^{*}+\beta^{*}\bm{\gamma}_{\widehat{\mathcal{S}}_{1}}^{*}$.

Recall the definition of the set of candidates passing joint thresholding,

$\widehat{\mathcal{S}}_{2}=\bigg{\\{}l:1\leq l\leq
s,|\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}|\geq\delta_{n}\|\bm{v}_{l}\|\frac{\|\bm{\mathrm{D}}-\bm{\mathbb{Z}}_{\widehat{\mathcal{S}}_{1}}\widetilde{\bm{\gamma}}_{\widehat{\mathcal{S}}_{1}}\|}{\sqrt{n}}\bigg{\\}}$

where $\delta_{n}=\sqrt{\omega\log s}$ and $\omega>0$ is a constant. In the
coming two sections, we explore the behavior of causal estimates for spurious
and relevant instruments in $\widehat{\mathcal{S}}_{2}$.

####  D.2.1 Limiting Behavior of Causal Effect Estimates from Spurious
Instruments

Now we give the proof of Theorem 3.1, which claims that bias of causal effect
estimates from spurious instruments that pass both marginal screening and
joint thresholding concentrate in a region with the center
$C_{*}=\mathrm{cov}(\bm{\alpha}_{D}^{*\mathrm{\scriptstyle
T}}\bm{U},\bm{\alpha}_{Y}^{*\mathrm{\scriptstyle
T}}\bm{U})/\mathrm{var}(\bm{\alpha}_{D}^{*\mathrm{\scriptstyle
T}}\bm{U}+\epsilon_{D})$. Without loss of generality, assume $\beta^{*}>0$.

###### Proof of Theorem 3.1.

Recall the definition of the set of irrelevant variables,
$\check{\mathcal{I}}=\\{l:1\leq l\leq
s,\gamma_{\widehat{\mathcal{S}}_{1},l}^{*}=0\\}$. To prove that there exists a
constant $d>0$ such that the causal estimates
$\widehat{\Gamma}_{\widehat{\mathcal{S}}_{1},l}/\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}$
for $l\in\check{\mathcal{I}}\cap\widehat{\mathcal{S}}_{2}$ concentrate in the
interval $[\beta^{*}+C_{*}-d,\beta^{*}+C_{*}+d]$ simultaneously, we need to
control the probability below

$P\bigg{(}\cup_{l\in\check{\mathcal{I}}\cap\widehat{\mathcal{S}}_{2}}\bigg{\\{}\bigg{|}\frac{\widehat{\Gamma}_{\widehat{\mathcal{S}}_{1},l}}{\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}}-\beta^{*}-C_{*}\bigg{|}>d\bigg{\\}}\bigg{)},$
(D.18)

where
$\widehat{\Gamma}_{\widehat{\mathcal{S}}_{1},l},\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}$
denote the debiased-lasso estimators.

It follows from (D.2) and $\pi_{\widehat{\mathcal{S}}_{1},l}^{*}=0$ for
irrelevant variables that

$\displaystyle\quad P\bigg{(}\bigg{\\{}\exists
l\in\check{\mathcal{I}}\cap\widehat{\mathcal{S}}_{2}:\bigg{|}\frac{\widehat{\Gamma}_{\widehat{\mathcal{S}}_{1},l}}{\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}}-\beta^{*}-C_{*}\bigg{|}>d\bigg{\\}}\cap\mathcal{E}\bigg{)}$
$\displaystyle\leq P\bigg{(}\bigg{\\{}\exists
l\in\check{\mathcal{I}}\cap\widehat{\mathcal{S}}_{2}:\bigg{|}\frac{e_{l}^{\mathrm{\scriptstyle
T}}(\Delta_{2}-\beta^{*}\Delta_{1})+\sqrt{n}\bm{v}_{l}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{U}}\bm{\alpha}_{Y}^{*}+\bm{\epsilon}_{Y})}{\sqrt{n}\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}}-C_{*}\bigg{|}>d\bigg{\\}}\cap\mathcal{E}\bigg{)}$
$\displaystyle=P\bigg{(}\bigg{\\{}\exists
l\in\check{\mathcal{I}}\cap\widehat{\mathcal{S}}_{2}:\bigg{|}\frac{e_{l}^{\mathrm{\scriptstyle
T}}(\Delta_{2}-\beta^{*}\Delta_{1})}{\sqrt{n}}+\bm{v}_{l}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{U}}\bm{\alpha}_{Y}^{*}+\bm{\epsilon}_{Y})-C_{*}\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}\bigg{|}>d|\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}|\bigg{\\}}\cap\mathcal{E}\bigg{)}$
$\displaystyle=P\bigg{(}\bigg{\\{}\exists
l\in\check{\mathcal{I}}\cap\widehat{\mathcal{S}}_{2}:\bigg{|}\frac{e_{l}^{\mathrm{\scriptstyle
T}}\\{\Delta_{2}-(\beta^{*}+C_{*})\Delta_{1}\\}}{\sqrt{n}}+\bm{v}_{l}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{U}}\bm{\alpha}_{Y}^{*}+\bm{\epsilon}_{Y})$
$\displaystyle\qquad\qquad\qquad-C_{*}\bm{v}_{l}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{U}}\bm{\alpha}_{D}^{*}+\bm{\epsilon}_{D})\bigg{|}>d|\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}|\bigg{\\}}\cap\mathcal{E}\bigg{)}$
$\displaystyle\leq\underbrace{P\bigg{(}\bigg{\\{}\exists
l\in\check{\mathcal{I}}\cap\widehat{\mathcal{S}}_{2}:\bigg{|}\frac{e_{l}^{\mathrm{\scriptstyle
T}}\\{\Delta_{2}-(\beta^{*}+C_{*})\Delta_{1}\\}}{\sqrt{n}}\bigg{|}>\frac{d|\widehat{\gamma}_{\widehat{\mathcal{S}}_{1},l}|}{2}\bigg{\\}}\cap\mathcal{E}\bigg{)}}_{\eqref{HDirrConc2}_{a}}$
$\displaystyle\quad+\underbrace{P\bigg{(}\bigg{\\{}\exists
l\in\check{\mathcal{I}}\cap\widehat{\mathcal{S}}_{2}:|\bm{v}_{l}^{\mathrm{\scriptstyle
T}}(\bm{\mathrm{U}}\bm{\alpha}_{Y}^{*}+\bm{\epsilon}_{Y})-C_{*}\bm{v}_{l}^{\mathrm{\scriptstyle