#  RobustIV and controlfunctionIV: Causal Inference for Linear and Nonlinear
Models with Invalid Instrumental Variables

Taehyeon Koo<EMAIL_ADDRESS>
Department of Statistics
Rutgers University
Piscataway, NJ 08854 Youjin Lee<EMAIL_ADDRESS>
Department of Biostatistics
Brown University
Providence, RI 02912 Dylan S. Small<EMAIL_ADDRESS>
Department of Statistics
The Wharton School, University of Pennsylvania
Philadelphia, PA 19104 Zijian Guo<EMAIL_ADDRESS>
Department of Statistics
Rutgers University
Piscataway, NJ 08854

###### Abstract

We present R software packages RobustIV and controlfunctionIV for causal
inference with possibly invalid instrumental variables. RobustIV focuses on
the linear outcome model. It implements the two-stage hard thresholding method
to select valid instrumental variables from a set of candidate instrumental
variables and make inferences for the causal effect in both low- and high-
dimensional settings. Furthermore, RobustIV implements the high-dimensional
endogeneity test and the searching and sampling method, a uniformly valid
inference method robust to errors in instrumental variable selection.
controlfunctionIV considers the nonlinear outcome model and makes inferences
about the causal effect based on the control function method. Our packages are
demonstrated using two publicly available economic data sets together with
applications to the Framingham Heart Study.

Keywords: Instrumental variable selection, confidence interval, nonlinear
outcome model, control function, maximum clique

##  1 Introduction

A common problem in making causal inferences from observational studies is
that there may be unmeasured confounders. The instrumental variable (IV)
method is one of the most useful methods to estimate the causal effect when
there might exist unmeasured confounding. The validity of IV methods relies on
that the constructed IVs satisfy the following three assumptions
simultaneously (Wooldridge, 2010, e.g.): conditioning the measured covariates,
(A1) the IVs are associated with the treatment;
(A2) the IVs are independent with the unmeasured confounders;
(A3) the IVs have no direct effect on the outcome.

The main challenge of applying IV-based methods in practice is that the
proposed IVs might not satisfy the above assumptions (A1)-(A3). For example,
in studying the causal effect of education on earning, the proximity of school
(Angrist and Krueger, 1991; Card, 1999) has been used as an instrumental
variable. However, this instrument might be related to other factors, such as
socioeconomic status, which could affect one’s earnings. Also, there might be
other advantages due to the proximity; for instance, people living close to
college could be more likely to be exposed to vocational programs linked to
colleges. So, the instrument could have a direct effect on earnings. In
addition, the problem of IVs not satisfying assumptions (A1) to (A3) is a
fundamental problem in Mendelian Randomization (MR), whose goal is to estimate
the causal effect of exposure on the disease by using genetic variants as
instruments. These genetic variants might violate assumptions (A2) and (A3)
due to pleiotropic effects (Bowden, Davey Smith, and Burgess, 2015, 2016;
Kang, Zhang, Cai, and Small, 2016).

This paper presents the R packages RobustIV and controlfunctionIV,
implementing robust causal inference approaches proposed in Guo, Kang, Cai,
and Small (2018a, b); Guo (2021); Guo and Small (2016); Li and Guo (2020). The
implemented inference methods choose the valid IVs among a set of candidate
IVs that may violate the assumptions (A2) and (A3). The proposed methods
target both linear and nonlinear causal effects. We also include the algorithm
implementation for settings with high-dimensional covariates and IVs.

In the package RobustIV, we implement robust and high-dimensional IV
algorithms for models assuming a constant and linear treatment effect. We
implement the Two Stage Hard Thresholding (TSHT) proposed in Guo et al.
(2018b), which selects valid IVs based on a voting method. The selected IVs
are then used to infer the linear treatment effect. Additionally, RobustIV
implements uniformly valid confidence intervals proposed in Guo (2021), which
guarantees valid coverage even if there are errors in selecting valid IVs.
RobustIV also contains the high-dimensional endogeneity test proposed in Guo
et al. (2018a), generalizing the Durbin-Wu-Hausman test (Durbin, 1954; Wu,
1973; Hausman, 1978).

We implement several control function methods in the package controlfunctionIV
to infer causal effects under nonlinear outcome models. We implement the
control function for the continuous outcome variable by showing it as the two-
stage least squares (TSLS) estimator with an augmented set of IVs (Guo and
Small, 2016). We further follow Guo and Small (2016) to test the validity of
the augmented set of IVs and construct the pretest estimator by comparing the
control function estimator and the TSLS estimator. We implement the probit
control function method for the binary outcome and make inferences for the
conditional average treatment effect (CATE) with possibly invalid IVs.
Moreover, the controlfunctionIV package implements the SpotIV method proposed
in Li and Guo (2020) for the semi-parametric outcome model with possibly
invalid IVs.

In R, there are well-developed IV methods when all IVs are assumed to satisfy
the assumptions (A1) to (A3), such as AER by Kleiber and Zeileis (2008) and
ivmodel by Kang, Jiang, Zhao, and Small (2020). The main difference in our
packages RobustIV and controlfunctionIV is that we allow for invalid IVs and
leverage the multiple IVs to learn the validity of the candidate IVs. We shall
mention other R packages implementing causal inference approaches with
possibly invalid IVs: sisvive implemented the method proposed in Kang et al.
(2016) to estimate the treatment effect under the majority rule; CIIV by
Windmeijer, Liang, Hartwig, and Bowden (2021) considered the causal inference
approaches with possibly invalid IVs for the low-dimensional linear outcome
model. In contrast, our packages RobustIV and controlfunctionIV are designed
under a broader framework by allowing for linear and nonlinear outcome models
with low- and high-dimensional IVs and covariates. Moreover, our package
RobustIV provides a uniformly valid confidence interval robust to the errors
in separating valid and invalid IVs.

The GitHub repository at https://github.com/bluosun/MR-GENIUS implemented the
MR Genius method (Tchetgen, Sun, and Walter, 2021), generalizing the method in
Lewbel (2012) and leveraging the heteroscedastic regression errors in the
treatment model to identify the causal parameter. The R package TSCI
implemented the two-stage curvature identification method proposed in Guo and
Bühlmann (2022), which leveraged the machine learning methods to capture the
nonlinearity in the treatment and identify the treatment effect with possibly
invalid instruments. In contrast, our packages use the different
identification conditions from Tchetgen et al. (2021); Lewbel (2012); Guo and
Bühlmann (2022).

The paper is organized as follows. In Section 2, we review the methods
implemented in RobustIV under the linear outcome model; in Section 3, we
discuss the inference approaches for the nonlinear outcome models implemented
in controlfunctionIV. In Section 4, we demonstrate the usage of RobustIV and
controlfunctionIV by analyzing economics data sets from Angrist and Krueger
(1991) and Mroz data. In Section 5, we demonstrate our packages with an MR
application to analyze the data from Framingham Heart Study (FHS).

#### Notation.

Let $\mathbb{R}^{p}$ be the set of real numbers with dimension $p$. For any
vector $v\in\mathbb{R}^{p}$, $v_{j}$ denotes its $j$th element, $v_{-j}$
denotes whole $v$ except for $j$-th index, and $\|v\|_{0}$ denotes the number
of non-zero elements in $v$. For any $n\times p$ matrix $M$, denote the
$(i,j)$ entry by $M_{ij}$, the $i$th row by $M_{i\cdot}$ , the $j$th column by
$M_{\cdot j}$, and the transpose of $M$ by $M^{\mkern-1.5mu\mathsf{T}}$; also,
$M_{IJ}$ denotes the submatrix of $M$ consisting of rows specified by the set
$I\subset\\{1,...,n\\}$ and columns specified by the set
$J\subset\\{1,...,p\\}$, $M_{I\cdot}$ denotes the submatrix of $M$ consisting
of rows indexed by the set I and all columns, and $M_{\cdot J}$ denotes the
submatrix of M consisting of columns specified by the set $J$ and all rows.
${\rm I}_{p}$ denotes $p\times p$ identity matrix. $\mathbbm{1}$ denotes the
indicator function. $\Phi$ denotes the CDF of the standard normal
distribution. For a sequence of random variable $X_{n}$, we use
$X_{n}\overset{d}{\rightarrow}X$ to denote that $X_{n}$ converges to X in
distribution.

##  2 Linear outcome models

Throughout the paper, we consider $n$ i.i.d. observations. For $1\leq i\leq
n,$ let $Y_{i}\in\mathbb{R}$, $D_{i}\in\mathbb{R},$
$Z_{i\cdot}\in\mathbb{R}^{p_{z}},$ and $X_{i\cdot}\in\mathbb{R}^{p_{x}}$
denote the outcome, the treatment, the instruments, and the baseline
covariates, respectively. This section reviews the robust instrumental
variable approaches in Guo et al. (2018a, b); Guo (2021), which are
implemented in the RobustIV package. We demonstrate the usage of RobustIV in
Sections 4.1 and 4.2.

###  2.1 Model assumption

We assume the following outcome model with possibly invalid IVs (Small, 2007;
Kang et al., 2016; Guo et al., 2018a; Windmeijer et al., 2021)

$\displaystyle
Y_{i}=D_{i}\beta+Z_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\pi+X_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\phi+\epsilon_{i},\quad$
$\displaystyle\mathbb{E}[\epsilon_{i}Z_{i\cdot}]=0,\
\mathbb{E}[\epsilon_{i}X_{i\cdot}]=0.$ (1)

This is the linear structural model in econometrics (Wooldridge, 2010). Here,
we aim to estimate the constant causal effect $\beta\in\mathbb{R}$. If $D_{i}$
is correlated with $\epsilon_{i}$ in the model (1), we say it is an endogenous
variable, and we cannot use popular estimators such as the OLS estimator. We
also assume the linear association model for the treatment

$\displaystyle
D_{i}=Z_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\gamma+X_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\psi+\delta_{i},\quad$
$\displaystyle\mathbb{E}[\delta_{i}Z_{i\cdot}]=0,\
\mathbb{E}[\delta_{i}X_{i\cdot}]=0.$ (2)

As a remark, the errors in (1) and (2) are allowed to be heteroscedastic. In
(1) and (2), $\pi_{j}=0$ if $j$-th IV satisfies the exclusion restriction
conditions (A2) and (A3), and $\gamma_{j}\neq 0$ if it satisfies the strong IV
assumption (A1).

We discuss the causal interpretation of the above model (1) using the
potential outcome framework (Small, 2007; Kang et al., 2016). Let
$Y_{i}^{(d,z)}$ be the potential outcome if individual $i$ were to receive the
treatment $d$ and the instruments $z$. For two possible values of the
treatment $d^{\prime},d$ and instruments $z^{\prime},z$, if we assume the
following potential outcomes model

$\displaystyle
Y_{i}^{(d^{\prime},z^{\prime})}-Y_{i}^{(d,z)}=(d^{\prime}-d)\beta+(z^{\prime}-z)^{\mkern-1.5mu\mathsf{T}}\kappa,\quad\mathbb{E}[Y_{i}^{(0,0)}|Z_{i\cdot},X_{i\cdot}]=X_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\phi+Z_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\eta,$
(3)

and define $\pi=\kappa+\eta$, and
$\epsilon_{i}=Y_{i}^{(0,0)}-\mathbb{E}[Y_{i}^{(0,0)}|Z_{i\cdot},X_{i\cdot}]$,
we obtain the model (1).

By combining (1) and (2), we obtain the reduced form models of $Y$ and $D$ as

$\displaystyle Y_{i}$
$\displaystyle=Z_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\Gamma+X_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\Psi+\xi_{i},\quad\mathbb{E}[\xi_{i}Z_{i\cdot}]=0,\mathbb{E}[\xi_{i}X_{i\cdot}]=0,$
(4) $\displaystyle D_{i}$
$\displaystyle=Z_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\gamma+X_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\psi+\delta_{i},\quad\mathbb{E}[\delta_{i}Z_{i\cdot}]=0,\mathbb{E}[\delta_{i}X_{i\cdot}]=0.$
(5)

Here, $\Gamma=\beta\gamma+\pi$, $\Psi=\beta\psi+\phi$ are reduced form
parameters and $\xi_{i}=\beta\delta_{i}+\epsilon_{i}$ is the reduced form
error term.

We introduce identifiability conditions for models (4) and (5). Let
$\mathcal{S}$ be the set of relevant IVs, i.e., $\mathcal{S}=\\{1\leq j\leq
p_{z}:\gamma_{j}\neq 0\\}$ and $\mathcal{V}$ be the set of relevant and valid
IVs, i.e., $\mathcal{V}=\\{j\in\mathcal{S}:\pi_{j}=0\\}$. The set
$\mathcal{S}$ contains all candidate IVs that are strongly associated with the
treatment. The set $\mathcal{V}$ is a subset of $\mathcal{S}$, which contains
all candidate IVs satisfying all classical IV assumptions. The main challenge
is that the set $\mathcal{V}$ is not known a priori in the data analysis.
Additional identifiability conditions are needed for identifying the causal
effect without any prior knowledge of $\mathcal{V}$. The majority rule is
introduced to identify causal effects with invalid IVs (Bowden et al., 2016;
Kang et al., 2016).

######  Condition 1 (Majority Rule).

More than half of the relevant IVs are valid: $|\mathcal{V}|>|\mathcal{S}|/2$.

The following plurality rule is a weaker identification condition than the
majority rule (Hartwig, Davey Smith, and Bowden, 2017; Guo, Kang, Cai, and
Small, 2018b).

######  Condition 2 (Plurality Rule).

The valid instruments form a plurality compared to the invalid instruments:
$|\mathcal{V}|>\max_{c\neq 0}|\\{j\in\mathcal{S}:\pi_{j}/\gamma_{j}=c\\}|.$

We present two inference methods for $\beta$ utilizing the majority and
plurality: two stage hard thresholding in Section 2.2 and searching and
sampling in Section 2.3. To present the methods, we consider the reduced form
estimators
$(\widehat{\Gamma}^{\intercal},\widehat{\gamma}^{\intercal})^{\intercal}$
satisfying

$\displaystyle\sqrt{n}\left(\begin{pmatrix}\widehat{\Gamma}\\\
\widehat{\gamma}\end{pmatrix}-\begin{pmatrix}\Gamma\\\
\gamma\end{pmatrix}\right)\overset{d}{\rightarrow}N_{2p_{z}}\left(0_{2p_{z}},\begin{pmatrix}\mathbf{V}^{\Gamma}&\mathbf{C}\\\
\mathbf{C}^{\mkern-1.5mu\mathsf{T}}&\mathbf{V}^{\gamma}\end{pmatrix}\right).$
(6)

We use $\widehat{\mathbf{V}}^{\Gamma},\widehat{\mathbf{C}},$ and
$\widehat{\mathbf{V}}^{\gamma}$ to denote consistent estimators of asymptotic
covariance matrix terms. In low dimensions, we estimate the reduced form
$({\Gamma}^{\intercal},{\gamma}^{\intercal})^{\intercal}$ by the OLS estimator
$(\widehat{\Gamma}^{\intercal},\widehat{\gamma}^{\intercal})^{\intercal}$ and
estimate the variance covariance matrices by sandwich estimators; see the
detailed construction in Section 2 of Guo (2021). In high-dimensional
settings, we can construct
$(\widehat{\Gamma}^{\intercal},\widehat{\gamma}^{\intercal})^{\intercal}$ as
the debiased Lasso estimator (Belloni, Chernozhukov, and Wang, 2011; Javanmard
and Montanari, 2014; Guo, Kang, Cai, and Small, 2018b); see more details in
Section 4.1 of Guo et al. (2018b).

###  2.2 Two stage hard thresholding (TSHT)

The TSHT consists of two steps: the first step is to screen out the weak IVs,
and the second step is to screen out invalid IVs. Specifically, the first step
of TSHT is to estimate the set $\mathcal{S}$ of relevant IVs by
$\widehat{\mathcal{S}}=\left\\{1\leq j\leq
p_{z}:|\widehat{\gamma}_{j}|\geq\lambda_{1}\sqrt{\widehat{\mathbf{V}}_{jj}^{\gamma}/n}\right\\},$
where $\lambda_{1}>0$ is a tuning parameter adjusting the testing
multiplicity.

The second thresholding step estimates the set $\mathcal{V}$ of valid
instruments. Our main strategy is to assume that one IV is valid and evaluate
whether the other IVs are valid from the point of view of that IV.
Particularly, for $j\in\widehat{\mathcal{S}}$, we assume the $j$-th IV to be
valid (i.e., $\pi_{j}=0$) and construct an estimator $\widehat{\pi}_{-j}$ of
$\pi_{-j}$ using the equation $\pi_{-j}=\Gamma_{-j}-\beta^{[j]}\gamma_{-j}$
with $\beta^{[j]}=\Gamma_{j}/\gamma_{j}$, and get the standard error of the
estimator. We test whether $\pi_{-j}=0$ by comparing $\widehat{\pi}_{-j}$ to a
threshold, calculated as multiplying the standard error of
$\widehat{\pi}_{-j}$ by a tuning parameter $\lambda_{2}$, which is a
Bonferroni correction adjusting for testing multiplicity. Using the above test
procedures, we construct a voting matrix
$\tilde{\Pi}\in\mathbb{R}^{|\widehat{\mathcal{S}}|\times|\widehat{\mathcal{S}}|}$
where $\tilde{\Pi}_{j,k}=1$ indicates that the $k$-th and $j$-th IVs agree
with each other to be valid. Finally, we get a symmetric voting matrix
$\widehat{\Pi}$ by setting
$\widehat{\Pi}_{j,k}=\min\\{\tilde{\Pi}_{j,k},\tilde{\Pi}_{k,j}\\}$.

Once we get $\widehat{\Pi}$, we estimate $\mathcal{V}$ by two options. Let
$\text{VM}_{k}$ denote the number of votes that the $k^{\rm th}$ IV, with
$k\in\widehat{\mathcal{S}}$, received from other candidates of IVs. First, we
define $\widehat{\mathcal{V}}$ by the set of IVs that receive a majority and a
plurality of votes (Guo et al., 2018b)

$\displaystyle\widehat{\mathcal{V}}^{\texttt{MP}}:=\\{k\in\widehat{\mathcal{S}}:\text{VM}_{k}>|\widehat{\mathcal{S}}|/2\\}\cup\\{k\in\widehat{\mathcal{S}}:\text{VM}_{k}=\max_{l\in\widehat{\mathcal{S}}}\text{VM}_{l}\\}.$
(7)

The next method is to estimate $\mathcal{V}$ by the maximum clique method. We
can generate a graph $\mathcal{G}$ with indexes belonging to
$\widehat{\mathcal{S}}$ and the adjacency matrix as $\widehat{\Pi}$. That is,
the indexes $j,k\in\widehat{\mathcal{S}}$ are connected if and only if
$\widehat{\Pi}_{j,k}=1$. Then as suggested in Windmeijer et al. (2021), we can
estimate $\widehat{\mathcal{V}}^{\texttt{MC}}$ as the maximum clique of the
graph $\mathcal{G}$, which is the largest fully connected sub-graph of
$\mathcal{G}$ (Csardi and Nepusz, 2006). Note that there might be several
maximum cliques. In this case, each maximum clique forms an estimator of
${\mathcal{V}}$ and our proposal reports several causal effect estimators
based on each maximum clique.

We further illustrate the definitions of $\widehat{\mathcal{V}}^{\texttt{MP}}$
and $\widehat{\mathcal{V}}^{\texttt{MC}}$ using the following example.
Consider $p_{z}=8$ with $\\{z_{1},z_{2},z_{3},z_{4}\\}$ being valid and
$\\{z_{5},z_{6},z_{7}\\}$ being invalid with the same invalidity level, and
$z_{8}$ being invalid IV with a different invalidity level. The left side of
Table 1 corresponds to an ideal setting where the valid and invalid IVs are
well separated and the valid IVs $\\{z_{1},z_{2},z_{3},z_{4}\\}$ only vote for
each other. In this case,
$\widehat{\mathcal{V}}^{\texttt{MC}}=\widehat{\mathcal{V}}^{\texttt{MP}}=\\{z_{1},z_{2},z_{3},z_{4}\\}$.
On the right side of Table 1, we consider the setting that the invalidity
level of $z_{5}$ might be mild and the IV $z_{5}$ receives the votes from
three valid IVs $\\{z_{2},z_{3},z_{4}\\}$. In this case,
$\widehat{\mathcal{V}}^{\texttt{MP}}=\\{z_{2},z_{3},z_{4},z_{5}\\}$. In
contrast, there are two maximum cliques $\\{z_{1},z_{2},z_{3},z_{4}\\}$ and
$\\{z_{2},z_{3},z_{4},z_{5}\\}$ and $\widehat{\mathcal{V}}^{\texttt{MC}}$ can
be either of these two.

| $z_{1}$ | $z_{2}$ | $z_{3}$ | $z_{4}$ | $z_{5}$ | $z_{6}$ | $z_{7}$ | $z_{8}$
---|---|---|---|---|---|---|---|---
$z_{1}$ | ✓ | ✓ | ✓ | ✓ | X | X | X | X
$z_{2}$ | ✓ | ✓ | ✓ | ✓ | X | X | X | X
$z_{3}$ | ✓ | ✓ | ✓ | ✓ | X | X | X | X
$z_{4}$ | ✓ | ✓ | ✓ | ✓ | X | X | X | X
$z_{5}$ | X | X | X | X | ✓ | ✓ | ✓ | X
$z_{6}$ | X | X | X | X | ✓ | ✓ | ✓ | X
$z_{7}$ | X | X | X | X | ✓ | ✓ | ✓ | X
$z_{8}$ | X | X | X | X | X | X | X | ✓
Votes | 4 | 4 | 4 | 4 | 3 | 3 | 3 | 1

| $z_{1}$ | $z_{2}$ | $z_{3}$ | $z_{4}$ | $z_{5}$ | $z_{6}$ | $z_{7}$ | $z_{8}$
---|---|---|---|---|---|---|---|---
$z_{1}$ | ✓ | ✓ | ✓ | ✓ | X | X | X | X
$z_{2}$ | ✓ | ✓ | ✓ | ✓ | ✓ | X | X | X
$z_{3}$ | ✓ | ✓ | ✓ | ✓ | ✓ | X | X | X
$z_{4}$ | ✓ | ✓ | ✓ | ✓ | ✓ | X | X | X
$z_{5}$ | X | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | X
$z_{6}$ | X | X | X | X | ✓ | ✓ | ✓ | X
$z_{7}$ | X | X | X | X | ✓ | ✓ | ✓ | X
$z_{8}$ | X | X | X | X | X | X | X | ✓
Votes | 4 | 5 | 5 | 5 | 6 | 3 | 3 | 1

Table 1: The left voting matrix $\widehat{\Pi}$ denotes that all valid IVs
$\\{z_{1},z_{2},z_{3},z_{4}\\}$ vote each other but not any other invalid IV.
The right voting matrix $\widehat{\Pi}$ denotes that the locally invalid IV
$z_{5}$ receives votes from valid IVs $\\{z_{2},z_{3},z_{4}\\}$ and invalid
IVs $\\{z_{6},z_{7}\\}$.

Once we have $\widehat{\mathcal{V}}$, we can construct an efficient point
estimator $\widehat{\beta}$ for $\beta$ in a low-dimensional setting via one-
step iteration as follows. First, we construct an initial estimator
$\tilde{\beta}=\frac{\widehat{\gamma}_{\widehat{\mathcal{V}}}^{\mkern-1.5mu\mathsf{T}}\tilde{A}\widehat{\Gamma}_{\widehat{\mathcal{V}}}}{\widehat{\gamma}_{\widehat{\mathcal{V}}}^{\mkern-1.5mu\mathsf{T}}\tilde{A}\widehat{\gamma}_{\widehat{\mathcal{V}}}},$
where
$\tilde{A}=\widehat{\Sigma}_{\widehat{\mathcal{V}},\widehat{\mathcal{V}}}-\widehat{\Sigma}_{\widehat{\mathcal{V}},\widehat{\mathcal{V}}^{c}}\widehat{\Sigma}_{\widehat{\mathcal{V}}^{c},\widehat{\mathcal{V}}^{c}}^{-1}\widehat{\Sigma}_{\widehat{\mathcal{V}}^{c},\widehat{\mathcal{V}}}$,
$\widehat{\Sigma}=\frac{1}{n}\sum_{i=1}^{n}W_{i\cdot}W_{i\cdot}^{\mkern-1.5mu\mathsf{T}}$,
and
$W_{i\cdot}=(Z_{i\cdot}^{\mkern-1.5mu\mathsf{T}},X_{i\cdot}^{\mkern-1.5mu\mathsf{T}})^{\mkern-1.5mu\mathsf{T}}$.

Next, we get a point estimator $\widehat{\beta}$ by one-step iteration
(Holland and Welsch, 1977)

$\displaystyle\widehat{\beta}=\frac{\widehat{\gamma}_{\widehat{\mathcal{V}}}^{\mkern-1.5mu\mathsf{T}}\widehat{A}\widehat{\Gamma}_{\widehat{\mathcal{V}}}}{\widehat{\gamma}_{\widehat{\mathcal{V}}}^{\mkern-1.5mu\mathsf{T}}\widehat{A}\widehat{\gamma}_{\widehat{\mathcal{V}}}},\quad\text{where}\quad\widehat{A}=[(\widehat{\mathbf{V}}^{\Gamma}-2\tilde{\beta}\widehat{\mathbf{C}}+\tilde{\beta}^{2}\widehat{\mathbf{V}}^{\gamma})_{\widehat{\mathcal{V}},\widehat{\mathcal{V}}}]^{-1}.$
(8)

Finally, the $1-\alpha$ confidence interval for $\beta$ is

$\displaystyle(\widehat{\beta}-z_{1-\alpha/2}\widehat{\rm
SE},\widehat{\beta}+z_{1-\alpha/2}\widehat{\rm
SE})\quad\text{where}\quad\widehat{\rm
SE}=\sqrt{\frac{\widehat{\gamma}_{\widehat{\mathcal{V}}}^{\mkern-1.5mu\mathsf{T}}\widehat{A}(\widehat{\mathbf{V}}^{\Gamma}-2\widehat{\beta}\widehat{\mathbf{C}}+\widehat{\beta}^{2}\widehat{\mathbf{V}}^{\gamma})_{\widehat{\mathcal{V}},\widehat{\mathcal{V}}}\widehat{A}\widehat{\gamma}_{\widehat{\mathcal{V}}}}{n(\widehat{\gamma}_{\widehat{\mathcal{V}}}^{\mkern-1.5mu\mathsf{T}}\widehat{A}\widehat{\gamma}_{\widehat{\mathcal{V}}})^{2}}}.$
(9)

As a remark, $\widehat{\mathbf{V}}^{\Gamma},\widehat{\mathbf{V}}^{\gamma},$
and $\widehat{\mathbf{C}}$ are heteroscedasticity-robust covariance estimators
and hence (9) is also robust to heteroscedastic errors in a low-dimensional
setting. In a high-dimensional setting, we set $\widehat{A}={\rm I}$ in (8),
and $\widehat{\mathbf{V}}^{\Gamma},\widehat{\mathbf{V}}^{\gamma},$ and
$\widehat{\mathbf{C}}$ are constructed under the homoscedastic error
assumptions; see more details in Guo et al. (2018b).

###  2.3 Searching and Sampling

We now review the searching and sampling method proposed in Guo (2021), which
provides uniformly valid conference intervals even if there are errors in
separating valid and invalid IVs. The right-hand side of Table 1 illustrates
an example of the invalid IVs not being separated from valid IVs in finite
samples. In the following, we review the idea of searching and sampling under
the majority rule and the more general method with the plurality rule can be
found in Guo (2021).

Let $\alpha\in(0,1)$ denote the pre-specified significance level. Given
$\beta\in\mathbb{R}$ and the reduced form estimator $\widehat{\Gamma}$ and
$\widehat{\gamma}$, we estimate $\pi_{j}$ with $j\in\widehat{\mathcal{S}}$ by

$\widehat{\pi}_{j}(\beta)=(\widehat{\Gamma}_{j}-\beta\widehat{\gamma}_{j})\mathbbm{1}(|\widehat{\Gamma}_{j}-\beta\widehat{\gamma}_{j}|\geq\widehat{\rho}_{j}(\beta,\alpha)),$
(10)

where
$\widehat{\rho}_{j}(\beta,\alpha)=\Phi^{-1}\left(1-\frac{\alpha}{2|\widehat{\mathcal{S}}|}\right)\widehat{\text{SE}}(\widehat{\Gamma}_{j}-\beta\widehat{\gamma}_{j})$
with $\widehat{\text{SE}}(\widehat{\Gamma}_{j}-\beta\widehat{\gamma}_{j})$
denoting a consistent estimator of the standard error of
$\widehat{\Gamma}_{j}-\beta\widehat{\gamma}_{j}$. We search for the value of
$\beta$ leading to enough valid IVs and construct the searching confidence
interval as

$\displaystyle\text{CI}^{\text{search}}=\left\\{\beta\in\mathbb{R}:\left\lVert\widehat{\pi}_{\widehat{\mathcal{S}}}(\beta)\right\rVert_{0}<|\widehat{\mathcal{S}}|/2\right\\},$
(11)

which collects all $\beta$ values such that more than half of IVs in
$\widehat{\mathcal{S}}$ are selected as valid.

Based on the searching method, Guo (2021) proposed a sampling confidence
interval, which retains the uniform coverage property and improves the
precision of the confidence interval. In particular, we sample

$\displaystyle\begin{pmatrix}\widehat{\Gamma}^{[m]}\\\
\widehat{\gamma}^{[m]}\end{pmatrix}\overset{iid}{\sim}N\left[\begin{pmatrix}\widehat{\Gamma}\\\
\widehat{\gamma}\end{pmatrix},\frac{1}{n}\begin{pmatrix}\widehat{\mathbf{V}}^{\Gamma}&\widehat{\mathbf{C}}\\\
\widehat{\mathbf{C}}^{\mkern-1.5mu\mathsf{T}}&\widehat{\mathbf{V}}^{\gamma}\end{pmatrix}\right],\quad\text{for}\quad
1\leq m\leq M.$

For $1\leq m\leq M$ and $j\in\widehat{\mathcal{S}}$, we modify (10) and define

$\widehat{\pi}_{j}^{[m]}(\beta,\lambda)=(\widehat{\Gamma}_{j}^{[m]}-\beta\widehat{\gamma}_{j}^{[m]})\mathbbm{1}(|\widehat{\Gamma}_{j}^{[m]}-\beta\widehat{\gamma}_{j}^{[m]}|\geq\lambda\cdot\widehat{\rho}_{j}(\beta,\alpha))$

with the shrinkage parameter $\lambda\asymp\left({\log
n}/{M}\right)^{\frac{1}{2|\widehat{\mathcal{S}}|}}.$ A data-dependent way of
choosing $\lambda$ can be found in Remark 3 of Guo (2021). For each $1\leq
m\leq M$, we construct a searching interval
$(\beta_{\min}^{[m]},\beta_{\max}^{[m]})$ where

$\beta_{\min}^{[m]}=\min_{\beta\in\mathcal{B}_{\lambda}^{[m]}}\beta\quad\text{and}\quad\beta_{\max}^{[m]}=\max_{\beta\in\mathcal{B}_{\lambda}^{[m]}}\beta$

with
$\mathcal{B}_{\lambda}^{[m]}=\left\\{\beta\in\mathbb{R}:\left\lVert\widehat{\pi}_{\widehat{\mathcal{S}}}^{[m]}(\beta,\lambda)\right\rVert_{0}<|\widehat{\mathcal{S}}|/2\right\\}$.
Then the sampling CI is defined as

$\text{CI}^{\text{sample}}=\left(\min_{m\in\mathcal{M}}\beta_{\min}^{[m]},\max_{m\in\mathcal{M}}\beta_{\max}^{[m]}\right).$
(12)

with $\mathcal{M}=\\{1\leq m\leq
M:(\beta_{\min}^{[m]},\beta_{\max}^{[m]})\neq\varnothing\\}.$ The sampling
confidence intervals in general improve the precision of the searching
confidence intervals. But both intervals can provide uniformly valid coverage
robust to the errors in separating valid and invalid IVs.

###  2.4 Endogeneity test in high dimensions

We review the high-dimensional endogeneity test proposed in Guo et al.
(2018a). We focus on the homoscedastic error setting by writing
$\Theta_{11}=\text{Var}[\xi_{i}|Z_{i\cdot},X_{i\cdot}]$,
$\Theta_{22}=\text{Var}[\delta_{i}|Z_{i\cdot},X_{i\cdot}]$, and
$\Theta_{12}=\text{Cov}[\xi_{i},\delta_{i}|Z_{i\cdot},X_{i\cdot}]$ for the
reduced form models (4) and (5). With the same estimators from TSHT in the
high-dimensional settings in Section 2.2, we can estimate the covariance
$\sigma_{12}$ by
$\widehat{\sigma}_{12}=\widehat{\Theta}_{12}-\widehat{\beta}\widehat{\Theta}_{22}$
where
$\widehat{\beta}=\frac{\sum_{j\in\widehat{\mathcal{V}}}\widehat{\gamma}_{j}\widehat{\Gamma}_{j}}{\sum_{j\in\widehat{\mathcal{V}}}\widehat{\gamma}_{j}^{2}}$
and $\widehat{\Theta}_{12}$ and $\widehat{\Theta}_{22}$ are consistent
estimators of the reduced form covariance ${\Theta}_{22}$ and $\Theta_{12}.$
We establish the asymptotic normality of $\widehat{\sigma}_{12}-{\sigma}_{12}$
in Guo et al. (2018a) and propose a testing procedure for
$H_{0}:{\sigma}_{12}=0.$

##  3 Nonlinear outcome models

This section reviews the control function IV methods (Guo and Small, 2016; Li
and Guo, 2020) implemented in the controlfunctionIV package, whose usage is
demonstrated in Sections 4.3 and 4.4.

###  3.1 Control function and pretest estimators

We consider the following nonlinear outcome and treatment models:

$\displaystyle Y_{i}$
$\displaystyle=\mathbf{G}(D_{i})^{\mkern-1.5mu\mathsf{T}}\bm{\beta}+X_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\phi+u_{i},\
\mathbb{E}[u_{i}Z_{i\cdot}]=\mathbb{E}[u_{i}X_{i\cdot}]=0,$ (13)
$\displaystyle D_{i}$
$\displaystyle=\mathbf{H}(Z_{i\cdot})^{\mkern-1.5mu\mathsf{T}}\bm{\gamma}+X_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\mathbf{\psi}+v_{i},\
\mathbb{E}[v_{i}Z_{i\cdot}]=\mathbb{E}[v_{i}X_{i\cdot}]=0,$ (14)

where
$\mathbf{G}(D_{i})=(D_{i},g_{2}(D_{i}),...,g_{k}(D_{i}))^{\mkern-1.5mu\mathsf{T}}$,
$\mathbf{H}(Z_{i\cdot})=(Z_{i\cdot},h_{2}(Z_{i\cdot}),...,h_{k}(Z_{i\cdot}))^{\mkern-1.5mu\mathsf{T}}$
with $\\{g_{j}(\cdot)\\}_{2\leq j\leq k}$ and $\\{h_{j}(\cdot)\\}_{2\leq j\leq
k}$ denoting the known nonlinear transformations. Under the models (13) and
(14), the IVs are assumed to be valid and the causal effect of increasing the
value of $D$ from $d_{2}$ to $d_{1}$ is defined as
$\mathbf{G}(d_{1})^{\mkern-1.5mu\mathsf{T}}\bm{\beta}-\mathbf{G}(d_{2})^{\mkern-1.5mu\mathsf{T}}\bm{\beta}.$

The control function (CF) method is a two-stage procedure. In the first stage,
regress $D$ on $\mathbf{H}(Z)$ and $X$, and obtain the predicted value
$\widehat{D}$ and its associated residual $\widehat{v}=D-\widehat{D}$. In the
second stage, we use $\widehat{v}$ as the proxies for the unmeasured
confounders and regress $Y$ on $\mathbf{G}(D)$, $X$, and $\widehat{v}$. We use
$\widehat{\beta}^{\text{CF}}$ to denote the estimated regression coefficient
corresponding to $D$. Guo and Small (2016) showed that
$\widehat{\beta}^{\text{CF}}$ is equivalent to the TSLS estimator with the
augmented set of IVs. Even if all IVs satisfy the classical assumptions
(A1)-(A3), there is no guarantee of the validity of the augmented IVs
generated by the CF estimator. Guo and Small (2016) applied the Hausman test
to assess the validity of the augmented set of IVs generated by the CF
estimator. The test statistic is defined as

$\displaystyle
H(\widehat{\beta}^{\text{CF}},\widehat{\beta}^{\text{TSLS}})=(\widehat{\beta}^{\text{CF}}-\widehat{\beta}^{\text{TSLS}})^{\mkern-1.5mu\mathsf{T}}[\text{Cov}(\widehat{\beta}^{\text{TSLS}})-\text{Cov}(\widehat{\beta}^{\text{CF}})]^{-}(\widehat{\beta}^{\text{CF}}-\widehat{\beta}^{\text{TSLS}}),$
(15)

where $\widehat{\beta}^{\text{TSLS}}$ is the two stage least square estimator,
$\text{Cov}(\widehat{\beta}^{\text{TSLS}})$ and
$\text{Cov}(\widehat{\beta}^{\text{CF}})$ are the covariance matrices of
$\widehat{\beta}^{\text{TSLS}}$ and $\widehat{\beta}^{\text{CF}}$, and $A^{-}$
denote the Moore-Penrose pseudoinverse. If the $p$-value
$P\left(\chi_{1}^{2}\geq
H(\widehat{\beta}^{\text{CF}},\widehat{\beta}^{\text{TSLS}})\right)$ is less
than $\alpha=0.05$, then we define the level $\alpha$ pretest estimator
$\widehat{\beta}^{\text{Pretest}}$ as $\widehat{\beta}^{\text{TSLS}}$;
otherwise, $\widehat{\beta}^{\text{CF}}$ defined above (Guo and Small, 2016).

###  3.2 Probit CF and SpotIV

We now consider the binary outcome model and continuous treatment model,

$\displaystyle\mathbb{E}\left[Y_{i}|D_{i}=d,W_{i\cdot}=w,u_{i}=u\right]=\mathbbm{1}(d\beta+w^{\mkern-1.5mu\mathsf{T}}\kappa+u>0),\quad\text{and}\quad
D_{i}=W_{i\cdot}^{\mkern-1.5mu\mathsf{T}}\gamma+v_{i},$ (16)

where
$W_{i\cdot}=(Z_{i\cdot}^{\mkern-1.5mu\mathsf{T}},X_{i\cdot}^{\mkern-1.5mu\mathsf{T}})^{\mkern-1.5mu\mathsf{T}}$,
the errors $(u_{i},v_{i})^{\intercal}$ are bivariate normal random variables
with zero means and independent of $W_{i\cdot}$,
$\kappa=(\kappa_{z}^{\mkern-1.5mu\mathsf{T}},\kappa_{x}^{\mkern-1.5mu\mathsf{T}})^{\mkern-1.5mu\mathsf{T}}$
is the coefficient vector of the IVs and measured covariates, and
$\gamma=(\gamma_{z}^{\mkern-1.5mu\mathsf{T}},\gamma_{x}^{\mkern-1.5mu\mathsf{T}})^{\mkern-1.5mu\mathsf{T}}$
is a parameter representing the association between $D_{i}$ and $W_{i\cdot}$.
When $\kappa_{z}\neq 0,$ the instruments are invalid. Since $u_{i}$ and
$v_{i}$ are bivariate normal, we write $u_{i}=\rho v_{i}+e_{i}.$ The model
(16) implies
$\mathbb{E}\left[Y_{i}|W_{i\cdot},v_{i}\right]=\Phi(D_{i}\beta^{*}+W_{i\cdot}^{\intercal}\Gamma^{*}+\rho^{*}v_{i})$
where $\beta^{*}=\beta/\sigma_{e}$,
$\Gamma^{*}=\kappa/\sigma_{e}+\beta^{*}\cdot\gamma$, and
$\rho^{*}=\rho/\sigma_{e}+\beta^{*}$ with $\sigma_{e}$ denoting the standard
error of $e_{i}=u_{i}-\rho v_{i}.$ That is, the conditional outcome model of
$Y_{i}$ given $W_{i,\cdot}$ and $v_{i}$ is a probit regression model.

Our goal is to estimate the conditional average treatment effect (CATE) from
$d_{2}$ to $d_{1}$

$\displaystyle\text{CATE}(d_{1},d_{2}|w)$
$\displaystyle:=\mathbb{E}[Y_{i}|D_{i}=d_{1},W_{i\cdot}=w]-\mathbb{E}[Y_{i}|D_{i}=d_{2},W_{i\cdot}=w].$
(17)

We first construct the OLS estimator $\widehat{\gamma}$ of $\gamma$. We
compute its residual $\widehat{v}=D-W\widehat{\gamma}$ and define
$\widehat{\Sigma}=\frac{1}{n}\sum_{i=1}^{n}W_{i\cdot}W_{i\cdot}^{\mkern-1.5mu\mathsf{T}}$.
We estimate $\mathcal{S}=\\{1\leq j\leq p_{z}:(\gamma_{z})_{j}\neq 0\\}$ by

$\displaystyle\widehat{\mathcal{S}}=\left\\{1\leq j\leq
p_{z}:|\widehat{\gamma}_{j}|\geq\widehat{\sigma}_{v}\sqrt{2\\{\widehat{\Sigma}^{-1}\\}_{j,j}\log
n/n}\right\\}$ (18)

with $\widehat{\sigma}_{v}^{2}=\sum_{i=1}^{n}\widehat{v}_{i}^{2}/n$. Next, as
CF in Section 3.1, we use $\widehat{v}$ as the proxy for unmeasured
confounders and implement the probit regression $Y$ on $W$ and $\widehat{v}$.
We use $\widehat{\Gamma}$ and $\widehat{\rho}$ to denote the probit regression
coefficients of $W$ of $\widehat{v}$ respectively. We apply the majority rule
and compute $\widehat{\beta}$ as the median of
$(\widehat{\Gamma}_{j}/\widehat{\gamma}_{j})_{j\in\widehat{\mathcal{S}}}.$ We
then estimate
$\widehat{\kappa}=\widehat{\Gamma}-\widehat{\gamma}\widehat{\beta}$. Finally,
we estimate CATE defined in (17) by the partial mean method (Newey, 1994;
Mammen, Rothe, and Schienle, 2012),

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[\Phi(d_{1}\widehat{\beta}+w^{\mkern-1.5mu\mathsf{T}}\widehat{\kappa}+\widehat{v}_{i}\widehat{\rho})\right]-\frac{1}{n}\sum_{i=1}^{n}\left[\Phi(d_{2}\widehat{\beta}+w^{\mkern-1.5mu\mathsf{T}}\widehat{\kappa}+\widehat{v}_{i}\widehat{\rho})\right]$

and construct the confidence interval by bootstrap (Li and Guo, 2020).

Li and Guo (2020) has proposed a more general methodology, named SpotIV, to
conduct robust causal inference with possibly invalid IVs. The model
considered in Li and Guo (2020) includes the probit outcome model in (16) as a
special case. In particular, Li and Guo (2020) replaced the known probit
transformation in (16) with the more general non-parametric function, which is
possibly unknown. Moreover, Li and Guo (2020) allows some instruments to be
correlated with the unmeasured confounders $u_{i}$ in the outcome model.

##  4 RobustIV and controlfunctonIV Usage

In this section, we illustrate the basic usage of RobustIV with the data set
from Angrist and Krueger (1991) and simulated high-dimensional data. Also, we
use the Mroz data set from Wooldridge (2010) to demonstrate the usage of
controlfunctionIV.

###  4.1 TSHT and SearchingSampling

In the following, we introduce usages of the R functions TSHT and
SearchingSampling with the data used in Angrist and Krueger (1991). Angrist
and Krueger (1991) studied the causal effect of the years of education (EDUC)
on the log weekly earnings (LWKLYWGE). Following Angrist and Krueger (1991),
we take 30 interactions (QTR120-QTR129, QTR220-QTR229, QTR320-QTR329) between
three quarter-of-birth dummies (QTR1-QTR3) and ten year-of-birth dummies
(YR20-YR29) as the instruments $Z$. For example, QTR120 is element-wise
product of QTR1 and YR20. Here, the quarter-of-birth dummies are the
indicators of whether the observed person was born in the first, second, and
third quarter of the year respectively, and the year-of-birth dummies are
indicators of which year the subject was born from 1940 to 1949 respectively.
We also include the following baseline covariates $X$: 9 year-of-birth dummies
(YR20-YR28), a race dummy (RACE), a marital status dummy (MARRIED), a dummy
for residence in an SMSA (SMSA), and eight region-of-residence dummies
(NEWENG, MIDATL, ENOCENT, WNOCENT, SOATL, ESOCENT, WSOCENT, MT). We first
apply the function TSHT.

⬇

R> Y <- as.vector(LWKLYWGE); D <- as.vector(EDUC)

R> Z <- sapply(paste0(”QTR”, c(seq(120,129), seq(220,229), seq(320,329))),

function(x){get(x)})

R> X <- cbind(sapply(paste0(”YR”,seq(20,28)),function(x){get(x)}),RACE,
MARRIED,

SMSA, NEWENG, MIDATL, ENOCENT, WNOCENT, SOATL, ESOCENT, WSOCENT, MT)

R> pz <- ncol(Z)

R> out.TSHT <- TSHT(Y=Y,D=D,Z=Z,X=X,

tuning.1st = sqrt(2.01*log(pz)), tuning.2nd = sqrt(2.01*log(pz)))

R> summary(out.TSHT)

betaHat Std.Error CI(2.50.0874 0.019 0.0502 0.1247 QTR120 QTR121 QTR122 QTR220
QTR222 QTR227 QTR322

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

Detected invalid IVs: QTR126 QTR226

Here, tuning.1st and tuning.2nd are tuning parameters $\lambda_{1}$ and
$\lambda_{2}$ used for the thresholds to get $\widehat{\mathcal{S}}$ and
$\widehat{\mathcal{V}}$ in Section 2.2 respectively. The default values for
these parameters are $\sqrt{\log n}$ in the low-dimensional setting. However,
in theory, any value above $\sqrt{2\log p}$ and diverging to infinity would
suffice. Since the data has 486926 observations, we choose the tuning
parameters as $\sqrt{2.01\log p_{z}}$ to avoid too conservative threshold
levels due to the huge sample. Once TSHT is implemented, we can call summary
to see the outputs of TSHT including the point estimator, its standard error,
confidence interval, and valid IVs as we discussed in Section 2.2.

The above result shows that TSHT selected QTR120, QTR121, QTR122, QTR220,
QTR222, QTR227, and QTR322 as valid IVs. Thus, valid IVs are interactions with
the first quarter of birth and dummies representing births in 1940, 1941, and
1942, interactions with the second quarter of birth and dummies representing
births in 1940, 1942, and 1947, and finally the interaction between the third
quarter-of-birth and dummy representing births in 1942. On the other hand, it
is reported that QTR126 and QTR226 are invalid IVs. That is, interactions with
the first and the second quarter of birth and dummy representing births in
1946 are relevant but invalid IVs. The remaining IVs have been screened out of
the first-stage selection as individually weak IVs.

The detection of invalid IVs implies that using whole $Z$ as valid IVs can
cause the estimate to be biased. In Angrist and Krueger (1991), the TSLS
estimate by using whole $Z$ as valid IVs is 0.0393. In contrast, our procedure
is more robust to the existence of possibly invalid IVs, giving the causal
estimate as 0.0874. Our $95\%$ confidence interval is above zero, indicating a
positive effect of education on earning.

In addition to the above output, the class object TSHT has other values that
are not reported by summary, for example, whether the majority rule is
satisfied or not, and the voting matrix to construct $\widehat{\mathcal{V}}$
in Section 2.2. These can be checked by directly calling TSHT.

As discussed in Section 2.2, there are different voting options to get
$\widehat{\mathcal{V}}$, where the default option voting = ’MaxClique’ stands
for $\widehat{\mathcal{V}}^{\texttt{MC}}$ and voting = ’MP’ stands for
$\widehat{\mathcal{V}}^{\texttt{MP}}$ in (7). If there are several maximum
cliques, summary returns results corresponding to each maximum clique.
Furthermore, since the default argument for which estimator to use is method =
’OLS’, one can choose other estimators by method = ’DeLasso’ for the debiased
Lasso estimator with SIHR R package (Rakshit, Cai, and Guo, 2021) and method =
’Fast.DeLasso’ for the fast computation of the debiased Lasso estimator
(Javanmard and Montanari, 2014). The above methods are useful in a high-
dimensional setting.

Next, we implement the uniformly valid confidence intervals by calling the
function SearchingSampling. We start with the searching CI defined in (11)
with the argument Sampling = FALSE.

⬇

R> out1 = SearchingSampling(Y=Y, D=D, Z=Z, X=X, Sampling=FALSE,

tuning.1st = sqrt(2.01*log(pz)), tuning.2nd = sqrt(2.01*log(pz)))

R> summary(out1)

Confidence Interval for Causal Effect: [-0.0964,0.2274]

With the default argument Sampling = TRUE, one can use the following code to
implement the more efficient sampling CI in (12).

⬇

R> set.seed(1)

R> out.SS = SearchingSampling(Y=Y, D=D, Z=Z, X=X,

tuning.1st = sqrt(2.01*log(pz)), tuning.2nd = sqrt(2.01*log(pz)))

R> summary(out.SS)

Confidence Interval for Causal Effect: [0.0135,0.1775]

The SearchingSampling confidence intervals are generally wider than that of
the TSHT since they are robust to the IV selection errors. The function
summary displays confidence interval for $\beta$, which are discussed in
Section 2.3. As in TSHT, one can use the argument method to employ the high-
dimensional debiased estimators instead of OLS.

###  4.2 endo.test

In the following, we show the usage of endo.test, a function for the
endogeneity test in high dimension with a simulated example. The corresponding
model and method are presented in Section 2.4. We consider the models (1) and
(2) and set $p_{z}=600$ with only the first 10 IVs being relevant. Among these
10 IVs, the first 3 IVs are invalid but the remaining IVs are valid. Moreover,
we set $\text{Corr}\left(\epsilon_{i},\delta_{i}\right)=0.8$, which indicates
a level of endogeneity. The function endo.test generates a class object with
same arguments in TSHT. The class object from endo.test can be used by calling
summary function, which enable us to see a brief result of ento.test.

⬇

R> set.seed(5)

R> n = 500; L = 600; s = 3; k = 10; px = 10; epsilonSigma =
matrix(c(1,0.8,0.8,1),2,2)

R> beta = 1; gamma = c(rep(1,k),rep(0,L-k))

R> phi = (1/px)*seq(1,px)+0.5; psi = (1/px)*seq(1,px)+1

R> Z = matrix(rnorm(n*L),n,L); X = matrix(rnorm(n*px),n,px);

R> epsilon = MASS::mvrnorm(n,rep(0,2),epsilonSigma)

R> D = 0.5 + Z R> Y = -0.5 + Z R> endo.test.model <- endo.test(Y,D,Z,X,
invalid = TRUE)

R> summary(endo.test.model)

P-value Test Valid IVs

0 H0 rejected Z4 Z5 Z6 Z7 Z8 Z9 Z10

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

Detected invalid IVs: Z1 Z2 Z3

When we call summary function, $p$-value, it reports the test result with
significance level $\alpha$ (default alpha = 0.05), the valid IVs, and
detected invalid IVs. H0 rejected means that the treatment is endogenous,
otherwise not. Since we set invalid = TRUE, ento.test allows some of IVs to be
invalid and conducts the endogeneity test with the selected
$\widehat{\mathcal{V}}$ defined in Section 2.2. With invalid = FALSE, the
function assumes that all IVs are valid. As in the above sections, one can use
method argument to employ other estimators both in the low and high
dimensions.

###  4.3 cf and pretest

In this section, we introduce usages of cf and pretest in the package
controlfunctionIV. The Mroz data was introduced in Mroz (1987) and then used
in various works of literature including Wooldridge (2010), which has $n=428$
individuals after removing the data with NA. Following Wooldridge (2010), we
estimate the causal effect of education on the log earnings of married working
women. The data is available in the Wooldridge package.

Here, the outcome $Y$ is log earnings (lwage), and the exposure $D$ is years
of schooling (educ). Moreover, there are other variables such as the father’s
education (fatheduc), the mother’s education (motheduc), the husband’s
education (huseduc), actual labor market experience (exper), its square
(expersq), and the women’s age (age).

Following Example 5.3 in Wooldridge (2010), we assume motheduc, fatheduc, and
huseduc to be valid IVs, denoted as
$Z_{i}=(Z_{i1},Z_{i2},Z_{i3})^{\mkern-1.5mu\mathsf{T}}$; we use and exper,
expersq, and age as baseline covariates, denoted as
$X_{i}=(X_{i1},X_{i2},X_{i3})^{\mkern-1.5mu\mathsf{T}}.$ Also assume that the
outcome and treatment models are (13) and (14) respectively with
$\mathbf{G}(D_{i})=(D_{i},D_{i}^{2})^{\mkern-1.5mu\mathsf{T}}$ and
$\mathbf{H}(Z_{i\cdot})=(Z_{i1},Z_{i2},Z_{i3},Z_{i1}^{2},Z_{i2}^{2},Z_{i3}^{2})^{\mkern-1.5mu\mathsf{T}}.$

Then we can implement the cf function by inputting a formula object, which has
the same form as that of ivreg in AER package. The function summary gives us
information on coefficients of the control function estimators, including the
point estimator, its standard error, $t$ value, and $p$ value.

⬇

R> library(wooldridge); library(controlfunctionIV); data(mroz); mroz <-
na.exclude(mroz)

R> Y <- mroz[,”lwage”]; D <- mroz[,”educ”]

R> Z <- as.matrix(mroz[,c(”motheduc”,”fatheduc”,”huseduc”)])

R> X <- as.matrix(mroz[,c(”exper”,”expersq”,”age”)])

R> cf.model <- cf(Y~D+I(D^2)+X|Z+I(Z^2)+X)

R> summary(cf.model)

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

Coefficients of the control function estimators:

\parEstimate Std.Error t value Pr(>|t|)

(Intercept) 1.2573907 0.7871438 1.597 0.055457 .

D -0.1434395 0.1102058 1.302 0.096884 .

I(D^2) 0.0086426 0.0041004 2.108 0.017817 *

Xexper 0.0438690 0.0131574 3.334 0.000465 ***

Xexpersq -0.0008713 0.0003984 2.187 0.014631 *

Xage -0.0011636 0.0048634 0.239 0.405511

—

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The following code infers the causal effect
$\mathbf{G}(d_{1})^{\mkern-1.5mu\mathsf{T}}\bm{\beta}-\mathbf{G}(d_{2})^{\mkern-1.5mu\mathsf{T}}\bm{\beta}$
by changing the treatment level from $d_{2}$ to $d_{1}=d_{2}+1$. Since the
second and third coefficients are related to $D$, we use the second and third
index to get the causal effect and its standard error.

⬇

R> d2 = median(D); d1 = median(D)+1;

R> D.diff <- c(d1,d1^2)-c(d2,d2^2); CE <- (D.diff)R> CE.sd <-sqrt(D.diffR>
CE.ci <- c(CE-qnorm(0.975)*CE.sd,CE+qnorm(0.975)*CE.sd)

R> cmat <- cbind(CE,CE.sd,CE.ci[1],CE.ci[2])

R> colnames(cmat)<-c(”Estimate”,”Std.Error”,”CI(2.5R> print(cmat, digits = 4)

Estimate Std.Error CI(2.5CE 0.07263 0.02171 0.03007 0.1152

The function pretest can be used to choose between the TSLS or the control
function method. If we run pretest with the same argument above and call
summary, it will output the following result:

⬇

R> pretest.model <- pretest(Y~D+I(D^2)+X|Z+I(Z^2)+X)

R> summary(pretest.model)

\parLevel 0.05 pretest estimator is control function estimator.

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

Coefficients of the pretest estimators:

\parEstimate Std.Error t value Pr(>|t|)

(Intercept) 1.2573907 0.7871438 1.597 0.055457 .

D -0.1434395 0.1102058 1.302 0.096884 .

I(D^2) 0.0086426 0.0041004 2.108 0.017817 *

Xexper 0.0438690 0.0131574 3.334 0.000465 ***

Xexpersq -0.0008713 0.0003984 2.187 0.014631 *

Xage -0.0011636 0.0048634 0.239 0.405511

—

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The first section of the output of summary reports which estimator is chosen
after the pretesting step. The second section lists brief information on
coefficients of pretest estimators including the point estimator, its standard
error, $t$ value, and $p$-value, similar to cf. Since the pretest estimator is
the control function estimator, the second section of summary is the same as
that of summary(cf.model).

###  4.4 Probit.cf

Finally, we conclude the usage part by looking at the usage of Probit.cf,
which is designed for the binary outcome with unmeasured confounders and
possibly invalid IVs. For illustration, we use the Mroz data and define the
binary outcome variable $Y_{0}$ to take the value 1 if the continuous outcome
$Y$ is greater than the median of $Y$ and 0 otherwise. We use the same
treatment variable $D$ as in the cf example. Contrary to the cf example, we
set the candidates of IVs $Z$ as motheduc, fatheduc, huseduc, exper, and
expersq, and assume that we have covariates $X$ as age.

We implement the Probit.cf function to estimate the CATE by increasing the
treatment value from the median of $D$ to the median plus one. We can call
summary to see the result of Probit.cf. The function summary provides
information on the valid IVs $\widehat{\mathcal{V}}$, the point estimator,
standard error, and $95\%$ confidence interval for $\beta$ in (16), and the
point estimator, the standard error, and $95\%$ confidence interval of CATE.

⬇

R> Z <- as.matrix(mroz[,c(”motheduc”,”fatheduc”,”huseduc”,”exper”,”expersq”)])

R> Y0 <- as.numeric((Y>median(Y)))

R> d2 = median(D); d1 = d2+1; w0 = apply(cbind(Z,X)[which(D == d2),], 2, mean)

R> Probit.model <- Probit.cf(Y0,D,Z,X,d1 = d1,d2 = d2,w0 = w0)

R> summary(Probit.model)

Estimate Std.Error CI(2.5Beta 0.2119 0.092 0.0316 0.3922 motheduc fatheduc
huseduc

CATE 0.0844 0.033 0.0198 0.1489 motheduc fatheduc huseduc

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

No invalid IV is detected

With the option invalid = TRUE, we allow invalid IVs and choose the valid IVs
among all provided IVs. If one wants to assume all IVs are valid, one can set
invalid = FALSE.

##  5 Application to Framingham Heart Study

We analyze the Framingham Heart Study (FHS) data and illustrate our package
using genetic variants as IVs. The FHS is an ongoing cohort study of
participants from the town of Framingham, Massachusetts, that has grown over
the years to include five cohorts with a total sample of over 15,000. The FHS,
initiated in 1948, is among the most critical sources of data on
cardiovascular epidemiology (Sytkowski, Kannel, and D’Agostino, 1990; Kannel,
2000; Mahmood, Levy, Vasan, and Wang, 2014). Since the late 1980s, researchers
across human health-related fields have used genetic factors underlying
cardiovascular diseases and other disorders. Over the last two decades, DNA
has been collected from blood samples and immortalized cell lines from members
of the Original Cohort, the Offspring Cohort, and the Third Generation Cohort
(Govindaraju et al., 2008). Several large-scale genotyping projects and
genome-wide linkage analysis have been conducted, and several other recent
collaborative projects have completed thousands of SNP genotypes for candidate
gene regions in subsets of FHS subjects with available DNA. The FHS has
recently been used for Mendelian Randomization to determine causal
relationships even in the presence of unmeasured confounding thanks to the
availability of genotype and phenotype data (Holmes et al., 2014; Dalbeth et
al., 2015; Hughes et al., 2014). As candidate IVs, we will use genotype data
from the FHS associated with the phenotype of interest and apply the proposed
methods described above.

We apply the RobustIV package to investigate the effect of low-density
lipoprotein (LDL-C) on globulin levels among individuals in the Framingham
Heart Study (FHS) Offspring Cohort, as was studied in Kang et al. (2020). We
use eight SNP genotypes (rs646776, rs693, rs2228671, rs2075650, rs4299376,
rs3764261, rs12916, rs2000999) that are known to be significantly associated
with LDL-C measured in mg/dL as candidate IVs (Kathiresan et al., 2007; Ma et
al., 2010; Smith et al., 2014). See Table 2 for details. The outcome of
interest $Y_{i}$ is a continuous globulin level ($g/L$) and the exposure
variable $D_{i}$ is the LDL-C level. Globulin is known to play a crucial role
in liver function, clotting, and the immune system. We also use the age and
sex of the subjects as covariates ${X}_{i\cdot}$. The study includes $n=1445$
subjects, with an average globulin level of $27.27$ (SD: $3.74$) and an
average LDL-C of $1.55$ (SD: 0.50). An average age is $35.58$ (SD: $9.74$) and
$54.95\%$ are males.

$Z_{j}$ | SNP | Position | lm(D $\sim$ Z) | lm(Y $\sim$ Z)
---|---|---|---|---
Estimate (Std. Error) |  $t$-statistic ($p$-value) | Estimate (Std. Error) |  $t$-statistic ($p$-value)
$Z_{1}$ | rs646776 | chr1:109275908 | -5.160 (1.610) | -3.205 (0.001) | -0.001 (0.170) | -0.007 (0.994)
$Z_{2}$ | rs693 | chr2:21009323 | -3.600 (1.286) | -2.799 (0.005) | 0.318 (0.135) | 2.349 (0.019)
$Z_{3}$ | rs2228671 | chr19:11100236 | 7.138 (2.029) | 3.518 ($<$0.001) | 0.529 (0.214) | 2.474 (0.014)
$Z_{4}$ | rs2075650 | chr19:44892362 | 8.451 (2.021) | 4.183 ($<$0.001) | 0.471 (0.213) | 2.208 (0.027)
$Z_{5}$ | rs4299376 | chr2:43845437 | 3.847 (1.387) | 2.773 (0.006) | 0.110 (0.146) | 0.752(0.452)
$Z_{6}$ | rs3764261 | chr16:56959412 | 3.651 (1.429) | 2.555 (0.011) | 0.275 (0.151) | 1.829 (0.067)
$Z_{7}$ | rs12916 | chr5:75360714 | 3.363 (1.365) | 2.463 (0.014) | -0.195 (0.144) | -1.357 (0.175)
$Z_{8}$ | rs2000999 | chr16:72074194 | -2.961 (1.629) | -1.818 (0.069) | -0.119 (0.172) | -0.691 (0.489)

Table 2:  Summary of the relationship between the single nucleotide
polymorphisms (SNPs) and low-density lipoprotein. The point estimator, its
standard error, $t$ value, and $p$-value are summary statistics from running a
marginal regression model specified in the column title. Position refers to
the position of the SNP in the chromosome, denoted as chr.

By applying endo.test, we detect one invalid IV and observe the evidence for
the existence of unmeasured confounders since the null hypothesis
$H_{0}:\sigma_{12}=0$ is rejected.

⬇

R> pz <- ncol(Z)

R> globulin.endo2 <- endo.test(Y.D,Z,X, invalid = TRUE,

tuning.1st = sqrt(2.01*log(pz)), tuning.2nd = sqrt(2.01*log(pz)))

R> summary(globulin.endo2)

P-value Test Valid IVs

0.0091 H0 rejected Z.1 Z.3 Z.4 Z.5 Z.6 Z.8

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

Detected invalid IVs: Z.2

Next, we implement TSHT with the default method of "OLS" under the low-
dimensional setting. Again, the same invalid IV is detected and the confidence
interval is above zero, indicating a positive effect of LDL on the glucose
level.

⬇

R> pz <- ncol(Z)

R> TSHT2 <- TSHT(Y, D, Z, X,

tuning.1st = sqrt(2.01*log(pz)), tuning.2nd = sqrt(2.01*log(pz)))

R> summary(TSHT2)

betaHat Std.Error CI(2.50.0529 0.0146 0.0243 0.0814 Z.1 Z.3 Z.4 Z.5 Z.6 Z.8

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

Detected invalid IVs: Z.2

We also constructed the confidence interval using the searching method, which
provides robustness to the IV selection errors.

⬇

R> SS1 <- SearchingSampling(Y, D, Z, X, tuning.1st = sqrt(2.01*log(pz)),

tuning.2nd = sqrt(2.01*log(pz)), Sampling = FALSE)

R> summary(SS1)

Confidence Interval for Causal Effect: [-0.2427,0.1894]

We further implement the sampling method, which leads to a shorter uniformly
valid CI than the searching method.

⬇

R> SS2 <- SearchingSampling(Y, D, Z, X, tuning.1st = sqrt(2.01*log(pz)),

tuning.2nd = sqrt(2.01*log(pz)), Sampling = TRUE)

R> summary(SS2)

Confidence Interval for Causal Effect: [-0.0521,0.1259]

In the following, we study nonlinear causal relationships using the
controlfunctionIV package. Burgess, Davies, and Thompson (2014) investigated a
nonlinear causal relationship between BMI and diverse cardiovascular risk
factors. Here we examine BMI’s possibly nonlinear causal effect on the insulin
level. Among $n=3733$ subjects, we excluded 618 subjects with missing
information on insulin level, and 50 subjects whose insulin level is greater
than $300\text{pmol/L}$ and whose BMI is greater than $45\text{kg/m}^{2}$. We
use log-transformed insulin as the outcome of interest $Y_{i}$ measured at
Exam 2. The exposure $D_{i}$ denotes the BMI measures at Exam 1. The
covariates $X_{i\cdot}$ that we adjusted for are age and sex. As valid IVs
$Z_{i\cdot}$, we propose using four SNP genotypes known to be significantly
associated with obesity. In our analysis, we include I(D^2) and I(X^2) to
account for quadratic effects of BMI, age, and sex on the outcome. We also
include I(Z^2) to account for possible quadratic effects of SNPs on the
exposure. The result from the pretest estimator is as follows:

⬇

R> insulin.pretest = pretest( Y ~ D + I(D^2) + X + I(X^2) | Z + I(Z^2) + X + I(X^2))

R> summary(insulin.pretest)

\parLevel 0.05 pretest estimator is control function estimator.

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

\parCoefficients of Pretest Estimators:

\parEstimate Std.Err t value Pr(>|t|)

(Intercept) 2.674e+00 6.006e-01 4.453 4.39e-06 ***

D 8.295e-02 2.828e-02 2.933 0.001690 **

I(D^2) -7.784e-04 2.742e-04 2.839 0.002276 **

X1 -1.780e-02 7.816e-03 2.277 0.011427 *

I(X1^2) 2.954e-04 8.852e-05 3.337 0.000428 ***

X2 -1.361e-01 5.654e-02 2.406 0.008087 **

—

Signif. codes:

0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The pretest estimator chooses the control function over the standard TSLS. The
results also show that BMI has a positive linear effect on the outcome but a
negative quadratic effect on the outcome.

## Acknowledgement

The research of T. Koo was supported in part by NIH grants R01GM140463 and
R01LM013614. The research of D. Small was supported in part by NIH grant
5R01AG065276-02.The research of Z. Guo was partly supported by the NSF grants
DMS 1811857 and 2015373 and NIH grants R01GM140463 and R01LM013614. Z. Guo is
grateful to Dr. Frank Windmeijer for bringing up the maximum clique method.

The Framingham Heart Study is conducted and supported by the National Heart,
Lung, and Blood Institute (NHLBI) in collaboration with Boston University
(Contract No. N01-HC-25195, HHSN268201500001I, and 75N92019D00031). This
manuscript was not prepared in collaboration with investigators of the
Framingham Heart Study and does not necessarily reflect the opinions or views
of the Framingham Heart Study, Boston University, or NHLBI. Funding for SHARe
Affymetrix genotyping was provided by NHLBI Contract N02-HL64278. SHARe
Illumina genotyping was provided under an agreement between Illumina and
Boston University. Funding for Affymetrix genotyping of the FHS Omni cohorts
was provided by Intramural NHLBI funds from Andrew D. Johnson and Christopher
J. O’Donnell.

## References

  * Angrist and Krueger (1991) Joshua D. Angrist and Alan B. Krueger.  Does compulsory school attendance affect schooling and earnings?  _The Quarterly Journal of Economics_ , 106(4):979–1014, 1991.
  * Belloni et al. (2011) A. Belloni, V. Chernozhukov, and L. Wang.  Square-root lasso: pivotal recovery of sparse signals via conic programming.  _Biometrika_ , 98(4):791–806, 2011.
  * Bowden et al. (2015) Jack Bowden, George Davey Smith, and Stephen Burgess.  Mendelian randomization with invalid instruments: Effect estimation and bias detection through egger regression.  _International Journal of Epidemiology_ , 44, 2015.
  * Bowden et al. (2016) Jack Bowden, George Davey Smith, Philip C. Haycock, and Stephen Burgess.  Consistent estimation in mendelian randomization with some invalid instruments using a weighted median estimator.  _Genetic Epidemiology_ , 40(4):304–314, 2016\.
  * Burgess et al. (2014) Stephen Burgess, Neil M Davies, and Simon G Thompson.  Instrumental variable analysis with a nonlinear exposure–outcome relationship.  _Epidemiology (Cambridge, Mass.)_ , 25(6):877, 2014.
  * Card (1999) David Card.  The causal effect of education on earnings.  In O. Ashenfelter and D. Card, editors, _Handbook of Labor Economics_ , volume 3, Part A, chapter 30, pages 1801–1863. Elsevier, 1 edition, 1999.
  * Csardi and Nepusz (2006) Gabor Csardi and Tamas Nepusz.  The igraph software package for complex network research.  _InterJournal_ , Complex Systems:1695, 2006.
  * Dalbeth et al. (2015) Nicola Dalbeth, Ruth Topless, Tanya Flynn, Murray Cadzow, Mark J Bolland, and Tony R Merriman.  Mendelian randomization analysis to examine for a causal effect of urate on bone mineral density.  _Journal of Bone and Mineral Research_ , 30(6):985–991, 2015.
  * Durbin (1954) J. Durbin.  Errors in variables.  _Revue de l’Institut International de Statistique / Review of the International Statistical Institute_ , 22(1/3):23–32, 1954\.
  * Govindaraju et al. (2008) Diddahally R Govindaraju, L Adrienne Cupples, William B Kannel, Christopher J O’Donnell, Larry D Atwood, Ralph B D’Agostino Sr, Caroline S Fox, Marty Larson, Daniel Levy, Joanne Murabito, et al.  Genetics of the framingham heart study population.  _Advances in genetics_ , 62:33–65, 2008.
  * Guo (2021) Zijian Guo.  Causal Inference with Invalid Instruments: Post-selection Problems and A Solution Using Searching and Sampling.  _arXiv e-prints_ , art. arXiv:2104.06911, 2021.
  * Guo and Bühlmann (2022) Zijian Guo and Peter Bühlmann.  Two stage curvature identification with machine learning: Causal inference with possibly invalid instrumental variables.  _arXiv preprint arXiv:2203.12808_ , 2022.
  * Guo and Small (2016) Zijian Guo and Dylan S. Small.  Control function instrumental variable estimation of nonlinear causal effect models.  _J. Mach. Learn. Res._ , 17(1):3448–3482, 2016\.
  * Guo et al. (2018a) Zijian Guo, Hyunseung Kang, T Tony Cai, and Dylan S Small.  Testing endogeneity with high dimensional covariates.  _Journal of Econometrics_ , 207(1):175–187, 2018a.
  * Guo et al. (2018b) Zijian Guo, Hyunseung Kang, Tony Cai, and Dylan S. Small.  Confidence intervals for causal effects with invalid instruments by using two-stage hard thresholding with voting.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 80(4):793–815, 2018b.
  * Hartwig et al. (2017) Fernando Pires Hartwig, George Davey Smith, and Jack Bowden.  Robust inference in summary data Mendelian randomization via the zero modal pleiotropy assumption.  _International Journal of Epidemiology_ , 46(6):1985–1998, 2017.
  * Hausman (1978) J. A. Hausman.  Specification tests in econometrics.  _Econometrica_ , 46(6):1251–1271, 1978.
  * Holland and Welsch (1977) Paul W Holland and Roy E Welsch.  Robust regression using iteratively reweighted least-squares.  _Communications in Statistics-theory and Methods_ , 6(9):813–827, 1977.
  * Holmes et al. (2014) Michael V Holmes, Leslie A Lange, Tom Palmer, Matthew B Lanktree, Kari E North, Berta Almoguera, Sarah Buxbaum, Hareesh R Chandrupatla, Clara C Elbers, Yiran Guo, et al.  Causal effects of body mass index on cardiometabolic traits and events: a mendelian randomization analysis.  _The American Journal of Human Genetics_ , 94(2):198–208, 2014.
  * Hughes et al. (2014) Kim Hughes, Tanya Flynn, Janak De Zoysa, Nicola Dalbeth, and Tony R Merriman.  Mendelian randomization analysis associates increased serum urate, due to genetic variation in uric acid transporters, with improved renal function.  _Kidney international_ , 85(2):344–351, 2014\.
  * Javanmard and Montanari (2014) Adel Javanmard and Andrea Montanari.  Confidence intervals and hypothesis testing for high-dimensional regression.  _Journal of Machine Learning Research_ , 15(82):2869–2909, 2014.
  * Kang et al. (2016) Hyunseung Kang, Anru Zhang, T. Tony Cai, and Dylan S. Small.  Instrumental variables estimation with some invalid instruments and its application to mendelian randomization.  _Journal of the American Statistical Association_ , 111(513):132–144, 2016.
  * Kang et al. (2020) Hyunseung Kang, Yang Jiang, Qingyuan Zhao, and Dylan S. Small.  ivmodel: An R Package for Inference and Sensitivity Analysis of Instrumental Variables Models with One Endogenous Variable.  _arXiv e-prints_ , art. arXiv:2002.08457, 2020.
  * Kang et al. (2020) Hyunseung Kang, Youjin Lee, T Tony Cai, and Dylan S Small.  Two robust tools for inference about causal effects with invalid instruments.  _Biometrics_ , 2020.
  * Kannel (2000) William B Kannel.  The framingham study: Its 50-year legacy and future promise.  _Journal of atherosclerosis and thrombosis_ , 6(2):60–66, 2000.
  * Kathiresan et al. (2007) Sekar Kathiresan, Alisa K Manning, Serkalem Demissie, Ralph B D’agostino, Aarti Surti, Candace Guiducci, Lauren Gianniny, Nöel P Burtt, Olle Melander, Marju Orho-Melander, et al.  A genome-wide association study for blood lipid phenotypes in the framingham heart study.  _BMC Medical Genetics_ , 8(1):S17, 2007.
  * Kleiber and Zeileis (2008) Christian Kleiber and Achim Zeileis.  _Applied Econometrics with R_.  Springer-Verlag, New York, 2008.
  * Lewbel (2012) Arthur Lewbel.  Using heteroscedasticity to identify and estimate mismeasured and endogenous regressor models.  _Journal of Business & Economic Statistics_, 30(1):67–80, 2012.
  * Li and Guo (2020) Sai Li and Zijian Guo.  Causal Inference for Nonlinear Outcome Models with Possibly Invalid Instrumental Variables.  _arXiv e-prints_ , art. arXiv:2010.09922, 2020.
  * Ma et al. (2010) Li Ma, Jing Yang, H Birali Runesha, Toshiko Tanaka, Luigi Ferrucci, Stefania Bandinelli, and Yang Da.  Genome-wide association analysis of total cholesterol and high-density lipoprotein cholesterol levels using the framingham heart study data.  _BMC Medical Genetics_ , 11(1):55, 2010.
  * Mahmood et al. (2014) Syed S Mahmood, Daniel Levy, Ramachandran S Vasan, and Thomas J Wang.  The framingham heart study and the epidemiology of cardiovascular disease: a historical perspective.  _The lancet_ , 383(9921):999–1008, 2014.
  * Mammen et al. (2012) Enno Mammen, Christoph Rothe, and Melanie Schienle.  Nonparametric regression with nonparametrically generated covariates.  _The Annals of Statistics_ , 40(2):1132–1170, 2012.
  * Mroz (1987) Thomas A. Mroz.  The sensitivity of an empirical model of married women’s hours of work to economic and statistical assumptions.  _Econometrica_ , 55(4):765–799, 1987.
  * Newey (1994) Whitney K Newey.  Kernel estimation of partial means and a general variance estimator.  _Econometric Theory_ , 10(2):1–21, 1994.
  * Rakshit et al. (2021) Prabrisha Rakshit, T. Tony Cai, and Zijian Guo.  SIHR: An R Package for Statistical Inference in High-dimensional Linear and Logistic Regression Models.  _arXiv e-prints_ , art. arXiv:2109.03365, 2021.
  * Small (2007) Dylan S Small.  Sensitivity analysis for instrumental variables regression with overidentifying restrictions.  _Journal of the American Statistical Association_ , 102(479):1049–1058, 2007.
  * Smith et al. (2014) J Gustav Smith, Kevin Luk, Christina-Alexandra Schulz, James C Engert, Ron Do, George Hindy, Gull Rukh, Line Dufresne, Peter Almgren, David S Owens, et al.  Association of low-density lipoprotein cholesterol–related genetic variants with aortic valve calcium and incident aortic stenosis.  _Journal of the American Medical Association_ , 312(17):1764–1771, 2014.
  * Sytkowski et al. (1990) Pamela A Sytkowski, William B Kannel, and Ralph B D’Agostino.  Changes in risk factors and the decline in mortality from cardiovascular disease: the framingham heart study.  _New England Journal of Medicine_ , 322(23):1635–1641, 1990.
  * Tchetgen et al. (2021) Eric Tchetgen Tchetgen, BaoLuo Sun, and Stefan Walter.  The GENIUS Approach to Robust Mendelian Randomization Inference.  _Statistical Science_ , 36(3):443 – 464, 2021\.
  * Windmeijer et al. (2021) Frank Windmeijer, Xiaoran Liang, Fernando P. Hartwig, and Jack Bowden.  The confidence interval method for selecting valid instrumental variables.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 83(4):752–776, 2021.
  * Wooldridge (2010) Jeffrey M Wooldridge.  _Econometric Analysis of Cross Section and Panel Data_ , volume 1 of _MIT Press Books_.  The MIT Press, 2010.
  * Wu (1973) De-Min Wu.  Alternative tests of independence between stochastic regressors and disturbances.  _Econometrica_ , 41(4):733–750, 1973.