Robustness Against Weak or Invalid Instruments: Exploring Nonlinear Treatment Models with Machine Learning

We aim to devise a robust IV framework that leads to causal identification even if the IVs proposed by domain experts may violate assumptions (A1) to (A3). Our framework provides a robustness guarantee even when all proposed IVs are invalid, including the most common regime with a single IV that is possibly invalid. It is well known that the treatment effect is not identifiable when there is no constraint on how the IVs violate assumptions (A2) and (A3). Our key identification assumption is that the violations of (A2) and (A3) arise from simpler forms than the association between the treatment and the IVs; that is, we exclude “special coincidences” such as the IVs violate (A2) and (A3) by linear forms and the conditional mean model of the treatment given IVs is also linear. This identification condition can be evaluated with the generalized IV strength measure introduced in Section 3.3 for a user-specific violation form.

We propose a novel two-stage curvature identification (TSCI) method to infer the treatment effect. An important operational step is to fit the treatment model with a machine learning (ML) algorithm, e.g., random forests, boosting, or deep neural networks. This first-stage ML model enhances the IV’s strength by learning a general conditional mean model of treatment given the IVs, instead of restricting to a linear association model. In the second stage, we leverage the ML prediction and adjust for different IV violation forms. A main novelty of our proposal is to estimate the bias resulting from high complexity of the first-stage ML and implement a corresponding bias correction step; see more details in Section 3.2.

We show that the TSCI estimator is asymptotically Gaussian when the generalized IV strength measure is sufficiently large. We further devise a data-dependent way of choosing the most suitable violation form among a collection of IV violation forms. By including the valid IV setting into the selection, our proposed general TSCI methodology does not exclude the possibility of the IVs being valid but provides a robust IV method against different invalidity forms.

To sum up, the contribution of the current paper is two-fold:

1. 1.

   TSCI explores a nonlinear treatment model with ML and leads to more reliable causal conclusions than existing methods assuming valid IVs.

2. 2.

   TSCI provides an efficient way of integrating the first-stage ML. A methodological novelty of TSCI is its bias correction in its second stage, which addresses the issue of high complexity and potential overfitting of ML.

There is relevant literature on causal inference when assumptions (A2) and (A3) are violated. When all IVs are invalid, Lewbel (2012); Tchetgen et al. (2021) proposed an identification strategy by assuming that the treatment model has heteroscedastic error but the covariance between the treatment and outcome model errors is homoscedastic. Our proposed TSCI is based on a different idea by exploring the nonlinear structure between the treatment and the IVs, not requiring anything on homo- and heteroscedasticity. Bowden et al. (2015) and Kolesár et al. (2015) assumed that the direct effect of the IVs on the outcome is nearly orthogonal to the effect of the IVs on the treatment. Han (2008); Bowden et al. (2016); Kang et al. (2016); Windmeijer et al. (2019); Guo et al. (2018); Windmeijer et al. (2021) considered the setting that a proportion of IVs are valid and conducted causal inference by selecting valid IVs in a data-dependent way. Along this line, Guo (2023) proposed a uniform inference procedure that is robust to the IV selection error. More recently, Sun et al. (2023) leveraged the number of valid IVs and used the interaction of genetic markers to identify the causal effect. In contrast to assuming the orthogonality and existence of valid IVs, our proposed TSCI is effective even if all IVs are invalid and the effect orthogonality condition does not hold. Such a setting is especially useful in handling econometric applications with a single IV.

The nonlinear structure of the treatment model has been explored in the econometric and the more recent ML literature. In the IV literature, Kelejian (1971); Amemiya (1974); Newey (1990); Hausman et al. (2012) considered constructing a non-parametric treatment model and Belloni et al. (2012) proposed constructing the optimal IVs with a Lasso-type first-stage estimator. More recently, Chen et al. (2023); Liu et al. (2020) proposed constructing the IV as the first-stage prediction given by the ML algorithm. All of the above works are focused on the valid IV setting. In contrast, our current paper applies to the more robust regime where the IVs are allowed to violate assumptions (A2) or (A3) in a range of forms. Even in the context of valid IVs, our proposed TSCI may lead to a more accurate estimator than directly using the predicted value for the treatment model as the new IV (Chen et al., 2023; Liu et al., 2020); see Section 3.5 for details.

ML algorithms integrated into the IV analysis to better accommodate the complicated outcome model. In Sections 3.5 and 5.1, we provide detailed comparisons to Double Machine Learning (DML) proposed in Chernozhukov et al. (2018). Athey et al. (2019) proposed the generalized random forests to infer heterogeneous treatment effects. Both works are based on valid IVs, while our current paper assumes a homogeneous treatment effect and focuses on the different robust IV framework with possibly invalid IVs.

Finally, our proposed TSCI methodology can be used to check the IV validity for the just-identification case. This is distinct from the Sargan test, J test, or specification test (Hansen, 1982; Sargan, 1958; Woutersen and Hausman, 2019) which are used to test IV validity in the over-identification case.

To explain the main idea, we start with the special case with no covariates,

$$Y_{i}=D_{i}\beta+Z_{i}^{\intercal}\pi+\epsilon_{i},\quad\text{and}\quad D_{i}=f(Z_{i})+\delta_{i},\quad\text{for}\quad 1\leq i\leq n,$$

(1)

where the errors $\epsilon_{i}$ and $\delta_{i}$ satisfy ${\mathbf{E}}(\epsilon_{i}\mid Z_{i})=0$ and ${\mathbf{E}}(\delta_{i}\mid Z_{i})=0.$ These errors $\epsilon_{i}$ and $\delta_{i}$ are typically correlated due to unobserved confounding. The outcome model in (1) is commonly used in the invalid IV literature (Small, 2007; Kang et al., 2016; Guo et al., 2018; Windmeijer et al., 2019), with $\pi\neq 0$ standing for the violation of IV assumptions (A2) and (A3). The treatment model in (1) is not a causal generating model between $D_{i}$ and $Z_{i}$ but only stands for a probability model between $D_{i}$ and $Z_{i}$. Specifically, the treatment assignment might be directly influenced by unmeasured confounders, but the treatment model in (1) does not explicitly state this generating process. Instead, for the joint distribution of $D_{i}$ and $Z_{i}$, we define the conditional expectation $f(Z_{i})={\mathbf{E}}(D_{i}\mid Z_{i}),$ which leads to the treatment model in (1). We discuss the structural equation model interpretation of the model (1) in Section 2.3, which explicitly characterizes how unmeasured confounders may affect the treatment assignment.

We introduce some projection notations to facilitate the discussion. For a random variable $U_{i}$, define its best linear approximation by $Z_{i}$ as $\gamma^{*}=\operatorname*{arg\,min}_{\gamma\in\mathbb{R}^{p_{z}+1}}{\mathbf{E}}(U_{i}-\widetilde{Z}_{i}^{\intercal}\gamma)^{2}$ with $\widetilde{Z}_{i}=(1,Z_{i}^{\intercal})^{\intercal}$. For $1\leq i\leq n$, write $\mathcal{P}_{Z_{i}}U_{i}=\widetilde{Z}_{i}^{\intercal}\gamma^{*}$ and $\mathcal{P}^{\perp}_{Z_{i}}U_{i}=U_{i}-\widetilde{Z}_{i}^{\intercal}\gamma^{*}.$ We apply $\mathcal{P}^{\perp}_{Z_{i}}$ to the model (1) and obtain $\mathcal{P}^{\perp}_{Z_{i}}Y_{i}=\mathcal{P}^{\perp}_{Z_{i}}D_{i}\beta+\mathcal{P}^{\perp}_{Z_{i}}\epsilon_{i}$, where the linear violation form $Z_{i}^{\intercal}\pi$ in (1) is removed after applying the transformation $\mathcal{P}^{\perp}_{Z_{i}}$. If $f$ is nonlinear, then ${\rm Var}\left(\mathcal{P}^{\perp}_{Z_{i}}f(Z_{i})\right)>0$ and the adjusted variable $\mathcal{P}^{\perp}_{Z_{i}}f(Z_{i})$ can be used as a valid IV to identify the effect $\beta$ via the following estimating equation

$${\mathbf{E}}\left[\mathcal{P}^{\perp}_{Z_{i}}f(Z_{i})(Y_{i}-D_{i}\beta)\right]={\mathbf{E}}\left[\mathcal{P}^{\perp}_{Z_{i}}f(Z_{i})(Z_{i}^{\intercal}\pi+\epsilon_{i})\right]=0.$$

(2)

Note that the last equality follows from the population orthogonality between $\mathcal{P}^{\perp}_{Z_{i}}f(Z_{i})$ and $Z_{i}$, together with ${\mathbf{E}}(\epsilon_{i}\mid Z_{i})=0.$

The estimation equation in (2) reveals that the treatment effect $\beta$ can be identified by exploring the nonlinear structure of $f$ under the model (1). We propose to learn $f$ using a ML prediction model and measure the variability of $\mathcal{P}^{\perp}_{Z_{i}}f(Z_{i})$ by defining a notation of generalized IV strength in Section 3.3. We shall emphasize that using ML algorithms to learn the possibly nonlinear function $f$ is critical in the current framework, allowing for invalid IVs. The ML algorithms help capture complicated structures in $f$ and typically retain enough variability contained in the adjusted variable $\mathcal{P}^{\perp}_{Z_{i}}f(Z_{i})$. In Section A.1 in the supplement, we present another example for causal identification in the context of the randomized experiments with non-compliance and invalid IVs.

We now introduce the general outcome model, allowing for more complicated violation forms than the linear violation in (1),

$$Y_{i}=D_{i}\beta+g(Z_{i},X_{i})+\epsilon_{i},\quad\text{with}\quad{\mathbf{E}}(\epsilon_{i}\mid X_{i},Z_{i})=0,\quad\text{for}\quad 1\leq i\leq n$$

(3)

where $\beta\in\mathbb{R}$ is the homogeneous treatment effect and the function $g:\mathbb{R}^{p_{x}+p_{z}}\rightarrow\mathbb{R}$ encodes how the IVs violate the assumptions (A2) and (A3). The classical valid IV setting requires that the function $g(z,x)$ does not change with different assignment of $z$. In the invalid IV literature, the commonly used outcome model with a linear violation form can be viewed as a special case of (3) when taking $g(Z_{i},X_{i})=Z_{i}^{\intercal}\pi_{z}+X_{i}^{\intercal}\pi_{x}.$

For the treatment, we define its conditional mean function $f(Z_{i},X_{i})\coloneqq{\mathbf{E}}(D_{i}\mid Z_{i},X_{i})$ for $1\leq i\leq n,$ leading to the following treatment model:

$$D_{i}=f(Z_{i},X_{i})+\delta_{i}\quad\text{with}\quad{\mathbf{E}}(\delta_{i}\mid Z_{i},X_{i})=0.$$

(4)

The model (4) is flexible as $f$ might be any unknown function of $Z_{i}$ and $X_{i}$, and the treatment variable can be continuous, binary, or a count variable. Similarly to (1), the model (4) is only a probability relation instead of a causal generating model and the treatment assignment is allowed to be influenced by unmeasured confounders.

It is well known that, for the outcome model (3), the identification of $\beta$ is generally impossible without additional identifiability conditions on the function $g$. In this paper, we allow the existence of invalid IVs but constrain the functional form of violating (A2) and (A3) to a pre-specified class. An operational step is to specify a set of basis functions for $g(\cdot)$ in the outcome model (3).

Define $\phi(X_{i})={\mathbf{E}}[g(Z_{i},X_{i})\mid X_{i}]$ and decompose $g$ as $g(Z_{i},X_{i})=h(Z_{i},X_{i})+\phi(X_{i})$ with $h(Z_{i},X_{i})=g(Z_{i},X_{i})-\phi(X_{i}).$ When the IVs $Z_{i}$ are valid, $g(Z_{i},X_{i})$ does not directly depend on the assignment of $Z_{i}$, which implies $\phi(X_{i})=g(Z_{i},X_{i})$ and $h(\cdot)=0.$ In this case, we just need to approximate $g(z,x)=\phi(x)$ by a set of basis functions $\overrightarrow{\bf w}(x)\in\mathbb{R}^{p_{w}}$. With $\overrightarrow{\bf w}(x)$, we approximate $g(Z_{i},X_{i})=\phi(X_{i})$ by a linear combination of $W_{i}=\overrightarrow{\bf w}(X_{i})\in\mathbb{R}^{p_{w}}$ for $1\leq i\leq n$ and $\{\phi(X_{i})\}_{1\leq i\leq n}$ by the matrix $W=(W^{\intercal}_{1},\ldots,W^{\intercal}_{n})^{\intercal}$. For example, if $\phi(x)$ is linear, we set $\overrightarrow{\bf w}(x)=\{1,x_{1},x_{2},\cdots,x_{p_{x}}\}$. Furthermore, we consider the additive model $\phi(x)=\sum_{j=1}^{p_{x}}\phi_{j}(x_{j})$ with smooth functions $\{\phi_{j}(\cdot)\}_{1\leq j\leq p_{x}}.$ For $1\leq j\leq p_{x},$ we construct a set of basis functions $\overrightarrow{{\mathbf{b}}_{j}}=\{b_{j,l}(\cdot)\}_{1\leq l\leq M_{j}}$ for $\phi_{j}(x_{j})$ with $M_{j}\geq 1$ denoting the number of basis functions. Examples of the basis functions include the polynomial or B spline basis. We then approximate $\phi(x)$ with $\overrightarrow{\bf w}(x)=\{1,\overrightarrow{{\mathbf{b}}_{1}}(x_{1}),\cdots,\overrightarrow{{\mathbf{b}}_{p_{x}}}(x_{p_{x}})\}.$

When the IV is invalid with $h\not\equiv 0$, in addition to approximating $\phi(x)$ by $\overrightarrow{\bf w}(x)$, we further approximate the function $h(z,x)$ by a set of pre-specified basis functions $v_{1}(z,x),\cdots,v_{L}(z,x)$ and then approximate $g(z,x)$ by

$$\mathcal{V}\coloneqq\{v_{1}(z,x),\cdots,v_{L}(z,x),\overrightarrow{\bf w}(x)\}\quad\text{with}\quad v_{l}:\mathbb{R}^{p_{x}+p_{z}}\rightarrow\mathbb{R}\quad\text{for}\quad 1\leq l\leq L.$$

(5)

We now present different choices of $g(\cdot)$ in (3), or equivalently, the basis set $\mathcal{V}$ in (5). With them, the model (3) can accommodate a wide range of valid or invalid IV forms.

1. (1)

   Linear violation. We take $g(z,x)=\sum_{l=1}^{p_{z}}\pi_{l}z_{l}+\phi(x)$ and $\mathcal{V}=\{z_{1},\cdots,z_{p_{z}},\overrightarrow{\bf w}(x)\}.$

2. (2)

   Polynomial violation (for a single IV). We consider $g(z,x)=h(z)+\phi(x)$, where the IV and baseline covariates affect the outcome in an additive manner. We set $\mathcal{V}=\{v_{1}(z),\cdots,v_{L}(z),\overrightarrow{\bf w}(x)\}$ and may take $v_{1}(z),\cdots,v_{L}(z)$ as (piecewise) polynomials of various orders.

3. (3)

   Interaction violation (for a single IV). We consider $g(z,x)=\sum_{l=1}^{p_{x}}\alpha_{l}\cdot z\cdot x_{l}+\phi(x)$, allowing the IV to affect the outcome interactively with the base covariates. For such a violation form, we set $\mathcal{V}=\{z\cdot x_{1},\cdots,z\cdot x_{p_{x}},\overrightarrow{\bf w}(x)\}.$

We focus on the single IV setting for the polynomial and interaction violation to simplify the notations. It can be extended to the setting with multiple IVs.

Two important remarks are in order for the choice of $g$ and $\mathcal{V}$. Firstly, the choice of $g$ and the corresponding basis set $\mathcal{V}$ is part of the outcome model specification (3), reflecting the users’ belief about how the IVs can potentially violate the assumptions (A2) and (A3). In applications, empirical researchers might apply domain knowledge and construct the pre-specified set of basis functions in (5). For example, Angrist and Lavy (1999) applied the Maimonides’ rule to construct the IV as the transformation of a variable “enrollment” and then adjust with the possible violation form generated by a linear, quadratic, or piecewise linear transformation of “enrollment”. Importantly, we do not have to assume the knowledge of $\mathcal{V}$ and will present in Section 3.4 a data-dependent way to choose the best $\mathcal{V}$ from a collection of candidate basis sets. Secondly, the linear violation form with $g(z,x)=\sum_{l=1}^{p_{z}}\pi_{l}z_{l}+\phi(x)$ is the most commonly used violation form in the context of the multiple IV framework (Small, 2007; Kang et al., 2016; Guo et al., 2018; Windmeijer et al., 2019). Most existing methods require at least a proportion of $Z_{i}$ to be valid, i.e., the corresponding $\pi_{l}$ being zero. In contrast, our framework allows all IVs to be invalid (i.e., all $\pi_{l}$ to be nonzero) by taking $\mathcal{V}=\{z_{1},\cdots,z_{p_{z}},\overrightarrow{\bf w}(x)\}.$

With the set ${\cal V}$ of basis functions in (5), we define the matrix $V$ which evaluates the functions in ${\cal V}$ at the observed variables:

$$V=\begin{pmatrix}V_{1}&\cdots&V_{n}\end{pmatrix}^{\intercal}\in\mathbb{R}^{n\times(L+p_{w})}\quad\text{with}\quad V_{i}=\left(v_{1}(Z_{i},X_{i}),\cdots,v_{L}(Z_{i},X_{i}),W_{i}^{\intercal}\right)^{\intercal},$$

(6)

for $1\leq i\leq n$. Note that we always use an intercept $W_{i,1}=1$. Instead of assuming $V_{i}$ perfectly generating $g(Z_{i},X_{i})$, we go with a broader framework by allowing an approximation error of $g(X_{i},Z_{i})$ by the best linear combination of $V_{i}$, defined as

$$R_{i}(V)\coloneqq g(Z_{i},X_{i})-V^{\intercal}_{i}\pi\quad\text{with}\quad\pi\coloneqq\operatorname*{arg\,min}_{b}{\mathbf{E}}\left[g(Z_{i},X_{i})-V_{i}^{\intercal}b\right]^{2}.$$

(7)

Define $R(V)=(R_{1}(V),\cdots,R_{n}(V))\in\mathbb{R}^{n}$. We shall omit the dependence on $V$ when there is no confusion. With (7), we express the outcome model (3) as

$$Y_{i}=D_{i}\beta+V_{i}^{\intercal}\pi+R_{i}(V)+\epsilon_{i},\quad\text{for}\quad 1\leq i\leq n.$$

(8)

If $g(Z_{i},X_{i})$ is well approximated by a linear combination of $V_{i},$ $\|R(V)\|_{2}/\sqrt{n}$ is close to zero or even $\|R(V)\|_{2}=0$.

We first provide the structural equation model (SEM) interpretation of the models (3) and (4). We consider the following SEM for $Y_{i}$ and $D_{i}$ and $X_{i}$,

$$Y_{i}\leftarrow a_{0}+D_{i}\beta+g_{1}(Z_{i},X_{i})+\nu_{1}(H_{i})+\epsilon^{0}_{i},\quad\text{and}\quad D_{i}\leftarrow f_{1}(Z_{i},X_{i})+\nu_{2}(H_{i})+\delta_{i}^{0},$$

(9)

where $H_{i}$ denotes some unmeasured confounders, $\epsilon^{0}_{i}$ and $\delta_{i}^{0}$ are random errors independent of $D_{i},Z_{i},X_{i},H_{i},$ and $a_{0}$ is the intercept such that ${\mathbf{E}}g_{1}(Z_{i},X_{i})={\mathbf{E}}\nu_{1}(H_{i})=0$. Define $g_{2}(Z_{i},X_{i})={\mathbf{E}}\left(\nu_{1}(H_{i})\mid Z_{i},X_{i}\right)$ and $f_{2}(Z_{i},X_{i})={\mathbf{E}}\left(\nu_{2}(H_{i})\mid Z_{i},X_{i}\right).$ In (9), the unmeasured confounders $H_{i}$ might affect both the outcome and treatment, $g_{1}$ encodes a direct effect of the IVs on the outcome, and $g_{2}$ captures the association between the IVs and the unmeasured confounders. The SEM (9) together with the definitions of $f_{2}$ and $g_{2}$ imply the outcome model (3) with $\epsilon_{i}=\epsilon^{0}_{i}+\nu_{1}(H_{i})-g_{2}(Z_{i},X_{i}).$ The treatment model $D_{i}=f(Z_{i},X_{i})+\delta_{i}$ arises with $f(Z_{i},X_{i})=f_{1}(Z_{i},X_{i})+f_{2}(Z_{i},X_{i})$ and $\delta_{i}=\delta^{0}_{i}+\nu_{2}(H_{i})-f_{2}(Z_{i},X_{i}).$

We now present how to interpret the invalid IV model with the potential outcome perspective. For the $i$-th subject with the baseline covariates $X_{i}$, we use $Y^{(z,d)}(X_{i})$ to denote the potential outcome with the IVs and the treatment assigned to $z\in\mathbb{R}^{p_{z}}$ and $d\in\mathbb{R}$, respectively. For $1\leq i\leq n,$ we consider the potential outcome model (Splawa-Neyman et al., 1990; Rubin, 1974),

$$Y_{i}^{(z,d)}(X_{i})=Y_{i}^{(0,0)}(X_{i})+d\beta+g_{1}(z,X_{i})\quad\text{and}\quad{\mathbf{E}}(Y_{i}^{(0,0)}(X_{i})\mid Z_{i},X_{i})=g_{2}(Z_{i},X_{i}),$$

(10)

where $\beta\in\mathbb{R}$ denotes the treatment effect, $g_{1}:\mathbb{R}^{p_{\rm z}+p_{\rm x}}\rightarrow\mathbb{R}$, and $g_{2}:\mathbb{R}^{p_{\rm z}+p_{\rm x}}\rightarrow\mathbb{R}.$ If $g_{1}(z,x)$ changes with $z$, the IVs directly affect the outcome, violating the assumption (A3). If $g_{2}(z,x)$ changes with $z$, the IVs are associated with unmeasured confounders, violating the assumption (A2). The model (10) extends the Additive LInear, Constant Effects (ALICE) model of Holland (1988) by allowing for a general class of invalid IVs. By the consistency assumption $Y_{i}=Y^{(Z_{i},D_{i})}_{i}(X_{i}),$ (10) implies (3) with $g(\cdot)=g_{1}(\cdot)+g_{2}(\cdot)$ and $\epsilon_{i}=Y_{i}^{(0,0)}(X_{i})-g_{2}(Z_{i},X_{i}).$ We can easily generalize (10) by considering $Y_{i}^{(z,d)}(X_{i})=Y_{i}^{(0,0)}(X_{i})+d\beta_{i}+g_{1}(z,X_{i}),$ where $\beta_{i}$ denotes the individual effect for the $i$-th individual. If we consider the random effect $\beta_{i}$ with ${\mathbf{E}}\beta_{i}=\beta$ and $\beta_{i}-\beta$ being independent of $(Z_{i},X_{i},D_{i}),$ then we obtain the model (3) with $\epsilon_{i}=Y_{i}^{(0,0)}(X_{i})-g_{2}(Z_{i},X_{i})+(\beta_{i}-\beta)D_{i}$ and our proposed method can be applied to make inference for $\beta={\mathbf{E}}\beta_{i}.$

We propose in Section 3 a novel methodology called two-stage curvature identification (TSCI) to identify $\beta$ under models (3) and (4). We first assume knowledge of the basis set $\mathcal{V}$ in (5) that spans the function $g(\cdot).$ The proposal consists of two stages: in the first stage, we employ ML algorithms to learn the conditional mean $f(\cdot)$ in the treatment model (4); in the second stage, we adjust for the IV invalidity form encoded by $\mathcal{V}$ and construct the confidence interval for $\beta$ by leveraging the first-stage ML fits.

The success of TSCI relies on the following identification condition.

The condition states that the basis set $\mathcal{V}$, used for generating $g(\cdot)$, does not fully span the conditional mean function $f(\cdot)$ of the treatment model. For any user-specified $\mathcal{V}$ in (5), Condition 1 can be partially examined by evaluating a generalized IV strength introduced in Section 3.3. When the generalized IV strength is sufficiently large, our proposed TSCI methodology leads to an estimator $\widehat{\beta}(\mathcal{V})$ in (19) such that

$$(\widehat{\beta}(\mathcal{V})-\beta)/\widehat{\rm SE}(\mathcal{V})\overset{d}{\to}N(0,1),$$

(11)

with the estimated standard error $\widehat{\rm SE}(\mathcal{V})$ depending on the generalized IV strength. Importantly, the basis set $\mathcal{V}$ used for generating $g(\cdot)$ does not need to be fixed in advance: we discuss in Section 3.4 a data-driven strategy for choosing a “good” basis set $\mathcal{V}_{\widehat{q}}$ among a collection of candidate sets $\mathcal{V}_{0}\subset\mathcal{V}_{1}\subset\cdots\subset\mathcal{V}_{Q}$ for a positive integer $Q$, where $\mathcal{V}_{0}$ corresponds to the valid IV setting. Our proposed TSCI compares estimators assuming valid IVs and adjusting for violation forms generated from the family $\{\mathcal{V}_{q}\}_{1\leq q\leq Q}.$ Thus, we encompass the setting of valid IVs and through this comparison, the TSCI methodology provides more robust causal inference tools than directly assuming IVs’ validity.

Intuitively, Condition 1 requires $f(\cdot)$ to be more non-linear than $g(\cdot)$, which cannot be fully tested since $g(\cdot)$ is unknown. We shall provide two important remarks on the plausibility of Condition 1. Firstly, when the proposed IVs are valid with $g(z,x)$ being independent of $z$, Condition 1 is satisfied as long as $f(z,x)$ depends on $z$. More interestingly, when $g(z,x)$ has a relatively weak dependence on $z$ (i.e., the IVs violate (A2) and (A3) locally), Condition 1 holds if the machine learning algorithm captures a strong dependence of $f(z,x)$ on $z$. In practice, domain scientists identify IVs based on their knowledge, and these IVs, even when violating assumptions (A2) and (A3), may have a relatively weak direct effect on the outcome or are weakly associated with the unmeasured confounders. In such a situation, TSCI provides robust causal inference by allowing the proposed IV to be invalid. However, we are not suggesting the users take any covariate as an invalid IV to implement our proposal. Secondly, even if Condition 1 does not hold, our proposed TSCI estimator is reduced to an estimator assuming valid IVs which then has a similar performance to TSLS and DML. In Section 5.3, we explore the performance of TSCI when Condition 1 does not hold.

We estimate the conditional mean $f(\cdot)$ in the treatment model (4) by a general ML fit. We randomly split the data into two disjoint subsets $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$. Throughout the paper we use $\mathcal{A}_{1}=\{1,2,\cdots,n_{1}\}$ with $n_{1}=\lfloor 2n/3\rfloor$ and set $\mathcal{A}_{2}=\{n_{1}+1,\cdots,n\}.$ Our results can be extended to any other splitting with $|\mathcal{A}_{1}|\asymp|\mathcal{A}_{2}|.$ Sample splitting is essential for removing the endogeneity of the ML predicted values. Without sample splitting, the ML predicted value for the treatment can be close to the treatment (due to overfitting) and hence highly correlated with unmeasured confounders, leading to a TSCI estimator suffering from a significant bias.

The main step is to construct the ML prediction model with data belonging to $\mathcal{A}_{2}$ and express the ML estimator of $f_{\mathcal{A}_{1}}=({f}(Z_{1},X_{1}),\cdots,{f}(Z_{n_{1}},X_{n_{1}}))^{\intercal}$ as

$$\widehat{f}_{\mathcal{A}_{1}}=\Omega D_{\mathcal{A}_{1}}\quad\text{for some matrix}\quad\Omega\in\mathbb{R}^{n_{1}\times n_{1}},$$

(12)

where the data belonging to $\mathcal{A}_{2}$ is used to train the ML algorithm and construct the transformation matrix $\Omega$. Our proposed TSCI method mainly relies on expressing the first-stage prediction as the linear estimator in (12). A wide range of first-stage ML algorithms are shown to have the expression in (12). In the following, we follow the results in Lin and Jeon (2006); Meinshausen (2006); Wager and Athey (2018) and express the sample split random forests in the form (12). In Sections A.6, A.7, and A.8 in the supplement, we present the explicit definitions of $\Omega$ for boosting, deep neural network (DNN), and approximation with B-splines, respectively.

Importantly, the transformation matrix $\Omega\in\mathbb{R}^{n_{1}\times n_{1}}$ has a related meaning as the hat matrix in linear regression procedures: however, $\Omega$ is not a projection matrix for nonlinear ML algorithms, which creates additional challenges in adopting the first-stage ML fit; see Section A.2.

The construction of honest trees and removal of self-prediction help remove the endogeneity of the ML predicted values $\{\widehat{f}(Z_{i},X_{i})\}_{i\in\mathcal{A}_{1}}$. Firstly, due to the sample-splitting, the construction of $\Omega$ does not directly depend on $\{D_{i}\}_{i\in\mathcal{A}_{1}},$ which removes the endogeneity contained in the ML predicted values. Secondly, consider an extreme setting with $Z_{i}$ and $X_{i}$ not being predictive for $D_{i}$. In such a scenario, without removing the self-prediction, the estimator $\widehat{f}(Z_{i},X_{i})$ will be dominated by its own treatment observation $D_{i}$ and $\widehat{f}(Z_{i},X_{i})$ is still highly correlated with the unmeasured confounders, leading to a biased TSCI estimator. The removal of self-prediction can be viewed as an ML version of the “leave-one-out” estimator for the treatment model, which was proposed in Angrist et al. (1999) to reduce the bias of TSLS with many IVs. In the numerical studies, we have observed that the exclusion of self-prediction leads to a more accurate estimator than the corresponding TSCI estimator with self-prediction, regardless of the IV strength; see Table S1 for the detailed comparison.

In the second stage, we leverage the first-stage ML fits in the form of (12) and adjust for the possible IV invalidity forms. In the remaining of the construction, we consider the outcome model in the form of (8) and assume that $V_{i}^{\intercal}\pi$ provides an accurate approximation of $g(Z_{i},X_{i})$, resulting in zero or sufficiently small $R_{i}(V)$. Hence $R_{i}(V)$ does not enter in the following construction of the estimator for $\beta$. We quantify the effect of the error term $R_{i}(V)$ in the theoretical analysis in Section 4 and propose in Section 3.4 a data-dependent way of choosing $V$ such that $R_{i}(V)$ is sufficiently small.

Applying the transformation $\Omega$ in (12) to the model (8) with data in $\mathcal{A}_{1}$, we obtain

$$\widehat{Y}_{\mathcal{A}_{1}}=\widehat{f}_{\mathcal{A}_{1}}\beta+\widehat{V}_{\mathcal{A}_{1}}\pi+\widehat{R}_{\mathcal{A}_{1}}+\widehat{\epsilon}_{\mathcal{A}_{1}},$$

(15)

where $\widehat{Y}_{\mathcal{A}_{1}}=\Omega Y_{\mathcal{A}_{1}},\widehat{f}_{\mathcal{A}_{1}}=\Omega D_{\mathcal{A}_{1}},\widehat{V}_{\mathcal{A}_{1}}=\Omega V_{\mathcal{A}_{1}},\widehat{R}_{\mathcal{A}_{1}}=\Omega R_{\mathcal{A}_{1}},$ and $\widehat{\epsilon}_{\mathcal{A}_{1}}=\Omega\epsilon_{\mathcal{A}_{1}}.$ For a matrix $V$, we use $P_{V}$ and $P_{V}^{\perp}$ to denote the projection to the column spaces of $V$ and its orthogonal complement, respectively. Based on (15), we project out the columns of $\widehat{V}_{\mathcal{A}_{1}}$ and estimate the effect $\beta$ by

$$\widehat{\beta}_{\rm init}(V)\coloneqq\frac{\widehat{Y}_{\mathcal{A}_{1}}^{\intercal}P^{\perp}_{\widehat{V}_{\mathcal{A}_{1}}}\widehat{f}_{\mathcal{A}_{1}}}{\widehat{f}_{\mathcal{A}_{1}}^{\intercal}P^{\perp}_{\widehat{V}_{\mathcal{A}_{1}}}\widehat{f}_{\mathcal{A}_{1}}}=\frac{{Y}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}}{{D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}}\quad\text{with}\quad\mathbf{M}(V)=\Omega^{\intercal}P^{\perp}_{\widehat{V}_{\mathcal{A}_{1}}}\Omega.$$

(16)

When $R_{i}(V)=0$, then we decompose the error of the above estimator as

$$\widehat{\beta}_{\rm init}(V)-\beta=\frac{{{\epsilon}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){\delta}_{\mathcal{A}_{1}}}}{{{D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}}}+\frac{{\epsilon}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){f}_{\mathcal{A}_{1}}}{{D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}}.$$

(17)

In the above decomposition, the first term on the right-hand side is a bias component appearing due to the use of the first-stage ML. To explain this, we consider the homoscedastic correlation ${\rm Cov}(\epsilon_{i},\delta_{i}\mid Z_{i},X_{i})={\rm Cov}(\epsilon_{i},\delta_{i})$ and obtain the explicit bias expression,

$$\frac{{\epsilon}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){\delta}_{\mathcal{A}_{1}}}{{D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}}\approx\frac{{\rm Cov}(\delta_{i},\epsilon_{i})\cdot{\rm Tr}[\mathbf{M}(V)]}{{D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}},$$

(18)

where ${\rm Tr}[\mathbf{M}(V)]$ denotes the trace of $\mathbf{M}(V)$ defined in (16). Note that ${\rm Tr}[\mathbf{M}(V)]$ can be viewed as degrees of freedom or complexity measure for the ML algorithm. It may become particularly large in the overfitting regime. The bias in (18) also becomes large when the denominator ${D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}$ is small, which means a weak generalized IV strength in the following (22). Thus, both the numerator and denominator in (18) can lead to a large bias of $\widehat{\beta}_{\rm init}(V)$. In Figure 1, we illustrate the bias of $\widehat{\beta}_{\rm init}(V)$ in settings with different values of the IV strength scaled to ${D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}$.

To address this finite-sample bias of $\widehat{\beta}_{\rm init}(V)$, we propose the following bias-corrected estimator,

$$\widehat{\beta}(V)=\frac{{Y}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}}{{D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}}-\frac{\sum_{i=1}^{n_{1}}[\mathbf{M}(V)]_{ii}\widehat{\delta}_{i}[\widehat{\epsilon}(V)]_{i}}{{{D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}}}.$$

(19)

with $\widehat{\delta}_{i}=D_{i}-\widehat{f}_{i}$ for $1\leq i\leq n_{1}$ and $\widehat{\epsilon}(V)=P^{\perp}_{V}[Y-D\widehat{\beta}_{\rm init}(V)]$. We show in Figure 1 that our proposal effectively corrects the bias of $\widehat{\beta}_{\rm init}(V)$. Importantly, our proposed bias correction is effective for both homoscedastic and heteroscedastic correlations. We present a simplified bias correction assuming homoscedastic correlation in Section A.10 in the supplement.

Centering at $\widehat{\beta}_{\rm RF}(V)$ defined in (19), we construct the confidence interval

$${\rm CI}(V)=\left(\widehat{\beta}(V)-z_{\alpha/2}\widehat{\rm SE}(V),\widehat{\beta}(V)+z_{\alpha/2}\widehat{\rm SE}(V)\right),$$

(20)

with $z_{\alpha/2}$ denoting the upper $\alpha/2$ quantile of the standard normal distribution and

$$\widehat{\rm SE}(V)=\frac{\sqrt{\sum_{i=1}^{n_{1}}[\widehat{\epsilon}(V)]^{2}_{i}[\mathbf{M}(V){D}_{\mathcal{A}_{1}}]_{i}^{2}}}{{D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}}\quad\text{with}\quad\widehat{\epsilon}(V)=P^{\perp}_{V}\left[Y-D\widehat{\beta}_{\rm init}(V)\right].$$

(21)

The IV strength is particularly important for identifying the treatment effect stably, and the weak IV is a major concern in practical applications of IV-based methods (Stock et al., 2002; Hansen et al., 2008). With a larger basis set $\mathcal{V}$, the IV strength will generally decrease as the information contained in $\mathcal{V}$ is projected out from the first-stage ML fit. It is crucial to introduce a generalized IV strength measure in consideration of invalid IV forms and the ML algorithm. Similarly to the classical setup, good performance of our proposed TSCI estimator also requires a relatively large generalized IV strength. We introduce a generalized IV strength measure as,

$$\mu(V)\coloneqq{{f}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){f}_{\mathcal{A}_{1}}}/{\left[\sum_{i\in\mathcal{A}_{1}}{\rm Var}(\delta_{i}\mid X_{i},Z_{i})/n_{1}\right]}.$$

(22)

If ${\rm Var}(\delta_{i}\mid X_{i},Z_{i})=\sigma_{\delta}^{2},$ then $\mu(V)$ is reduced to ${f}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){f}_{\mathcal{A}_{1}}/\sigma_{\delta}^{2}.$ A sufficiently large strength $\mu(V)$ will guarantee stable point and interval estimators defined in (19) and (20). Hence, we need to check whether $\mu(V)$ is sufficiently large. Since $f$ is unknown, we estimate ${f}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){f}_{\mathcal{A}_{1}}$ by its sample version ${D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}$ and estimate $\mu(V)$ by $\widehat{\mu(V)}\coloneqq{{D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}}/\left[{\|D_{\mathcal{A}_{1}}-\widehat{f}_{\mathcal{A}_{1}}\|_{2}^{2}/n_{1}}\right].$

In Section 4, we show that our point estimator in (19) is consistent when the IV strength $\mu(V)$ is much larger than ${\rm Tr}\left[\mathbf{M}(V)\right]$; see Condition (R2). We now develop a bootstrap test to provide an empirical assessment of this IV strength requirement. Since Rothenberg (1984) and Stock et al. (2002) suggested the concentration parameter being larger than $10$ as being “adequate”, we develop a bootstrap test for $\mu(V)\geq\max\{2{\rm Tr}\left[\mathbf{M}(V)\right],10\}.$ We apply the wild bootstrap method and construct a probabilistic upper bound for the estimation error ${D}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){D}_{\mathcal{A}_{1}}-{f}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){f}_{\mathcal{A}_{1}}=2{f}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){\delta}_{\mathcal{A}_{1}}+{\delta}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){\delta}_{\mathcal{A}_{1}}.$ For $1\leq i\leq n_{1}$, we define $\widehat{\delta}_{i}=D_{i}-\widehat{f}_{i}$ and compute $\widetilde{\delta}_{i}=\widehat{\delta}_{i}-\bar{\mu}_{\delta}$ with $\bar{\mu}_{\delta}=\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\widehat{\delta}_{i}.$ For $1\leq l\leq L,$ we generate $\delta^{[l]}_{i}=U^{[l]}_{i}\cdot\widetilde{\delta}_{i}$ for $1\leq i\leq n_{1},$ with $\{U^{[l]}_{i}\}_{1\leq i\leq n_{1}}$ generated as i.i.d. standard normal random variables, and compute $S^{[l]}=\left[{2{f}_{\mathcal{A}_{1}}^{\intercal}\mathbf{M}(V){\delta}^{[l]}+({\delta}^{[l]})^{\intercal}\mathbf{M}(V){\delta}^{[l]}}\right]/\left[{{{\|D_{\mathcal{A}_{1}}-\widehat{f}_{\mathcal{A}_{1}}\|_{2}^{2}/n_{1}}}}\right].$ We use $\mathcal{S}_{\alpha_{0}}(V)$ to denote the upper $\alpha_{0}$ empirical quantile of $\{|S^{[l]}|\}_{1\leq l\leq L}$. We conduct the generalized IV strength test $\widehat{\mu(V)}\geq\max\{2{\rm Tr}\left[\mathbf{M}(V)\right],10\}+{\mathcal{S}_{\alpha_{0}}(V)},$ with ${\mathcal{S}_{\alpha_{0}}(V)}$ being a high probability upper bound for $|\widehat{\mu(V)}-\mu(V)|.$ We use $\alpha_{0}=0.025$ throughout this paper. If the above generalized IV strength test is passed, the IV is claimed to be strong after adjusting for the matrix $V$ defined in (6); otherwise, the IV is claimed to be weak after adjusting for the matrix $V.$ Empirically, we observe reliable inference properties when the estimated generalized IV strength $\widehat{\mu(V)}$ is above 40; see Figure S1 for details.