Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
48,643 characters · 12 sections · 82 citation commands
Dual Instrumental Variable Regression
Inferring causal relationships under the influence of unobserved confounders remains one of the most challenging problems in economics, health care, and social sciences Angrist08:Harmless,Imbens15:CIS. A typical example in economics is the study of returns from schooling Card99:Labor, which attempts to measure the causal effect of education on labor market earnings. For each individual, the treatment variable $X$ represents the level of education and the outcome $Y$ represents how much they earn. However, one's level of education and income is likely confounded by the socioeconomic status or other unobserved confounding factors $H$ Angrist08:Harmless.
Since randomized control trials are often infeasible in most economic studies, economists have turned to instrumental variables (IVs) or instruments derived from naturally occurring random experiments to overcome unobserved confounding. Informally, instrumental variables $Z$ are defined as variables that are associated with the treatment $X$, affect the outcome $Y$ only through $X$ and do not share common causes with $Y$. For instance, the season-of-birth was used as an instrument in Angrist91:Schooling to estimate the impact of compulsory schooling on earnings. Because of the compulsory school attendance laws, an individual's season-of-birth, which is likely to be random, affects how long they actually remain in school, but not their earnings. Figure (ref) illustrates this example. Finding valid instruments for specific problems is an essential task in econometrics Angrist08:Harmless and epidemiology Burgess17:Mendelian.
Although IV analysis is widely used, the statistical tools employed for estimating causal effect are fairly rudimentary. Most applications of instrumental variables utilise a two-stage procedure Angrist08:Harmless,White82:IV,Singh19:KIV,Hartford17:DIV. For instance, the two-stage least squares (2SLS) relies on the assumption that the relationship between $X$ and $Y$ is linear peters2017elements. It first estimates the conditional mean $\mathbb{E}[X|Z=z]$ via linear regression and then regresses $Y$ on the estimate of $\mathbb{E}[X|Z=z]$ to obtain an estimate of the causal effect. Since the first-stage estimate is by construction independent from confounders, the resultant causal estimate is therefore free from hidden confounding. In the non-linear setting, however, a poorly-fitted first-stage regression may result in inaccurate second-stage estimates Angrist08:Harmless.
In this paper, we propose a novel procedure, DualIV, to directly estimate the structural causal function. Unlike previous works which extend 2SLS by employing non-linear models in place of their linear counterparts Hartford17:DIV,Singh19:KIV, we solve the dual problem which can be expressed as a convex-concave saddle-point problem. Based on this framework, we develop a consistent reproducing kernel Hilbert spaces-based (RKHS) algorithm. Our formulation was inspired by the mathematical resemblance of non-linear IV to two-stage problems in stochastic programming Shapiro14:LSP,Dai17:Dual,Hsu19.
The rest of the paper is organized in the following manner. Section (ref) introduces the IV regression problem, reviews related work and identifies current limitations. We present our formulation in Section (ref), followed by the kernelized estimation method in Section (ref). Then, Section (ref) reports empirical comparisons between DualIV and existing algorithms. Finally, we discuss the limitations of our procedure and suggest future directions in Section (ref). All proofs can be found in Appendix (ref).
Let $X$, $Y$, and $Z$ be treatment, outcome, and instrumental variable(s) taking values in $\ensuremath{\mathcal{X}}$, $\ensuremath{\mathcal{Y}}$, and $\ensuremath{\mathcal{Z}}$, respectively. In this work, we assume that $\ensuremath{\mathcal{Y}} \in \mathbb{R}$, and $\ensuremath{\mathcal{X}}$ and $\ensuremath{\mathcal{Z}}$ are Polish spaces. We also assume that $Y$ is bounded, i.e., $|Y| < M < \infty$ almost surely. Moreover, we denote unobserved confounder(s) by $H$. The underlying data generating process (DGP) is described by the causal graph in Figure (ref) equipped with the following structural causal model (SCM):
where $f$ is an unknown, potentially non-linear continuous function and $\varepsilon$ denotes the additive noise which depends on the hidden confounder(s) $H$. If $\mathbb{E}[\varepsilon \,|\, X]=0$, we can estimate $f$ consistently from observational data via the standard least-square regression. This allows us to identify $\mathbb{E}[Y \,|\, \text{do}(X=x)]$ where $\text{do}(X=x)$ represents an intervention on $X$ where its value is set to $x$ Pearl2000.
For non-expert readers, we elaborate that $\text{do}(X=x)$ here denotes a mathematical operator which simulates physical interventions by setting the value of $X$ to $x$, while keeping the rest of the model unchanged Pearl09:Causal. That is, the conditional expectation $\mathbb{E}[Y\,|\,\text{do}(X=x)]$ is computed with respect to the interventional distribution $\ensuremath{\mathbb{P}}(Y\,|\,\text{do}(X=x))$. We can estimate $\ensuremath{\mathbb{P}}(Y\,|\,\text{do}(X=x))$ if it is possible to directly manipulate $X$ and then observe the resulting outcome $Y$. In Figure (ref), for instance, one may assign different levels of education to people and then observe their subsequent levels of income in the labor market. Unfortunately, such experiment is not always possible and we only have access to an observational distribution $\ensuremath{\mathbb{P}}(Y\,|\,X=x)$, which can be different from $\ensuremath{\mathbb{P}}(Y\,|\,\text{do}(X=x))$. The discrepancy between interventional and observational distributions may result from the unobserved socioeconomic status, as illustrated in Figure (ref).
When hidden confounders exist between $X$ and $Y$, the error term $\varepsilon$ in (ref) is generally correlated with $X$. Hence, $\mathbb{E}[\varepsilon|X]\neq 0$ and it follows from (ref) that
which implies that $\mathbb{E}[Y \,|\, \text{do}(X=x)] \neq \mathbb{E}[Y\,|\, X=x]$. Thus, standard least-square regression no longer provides a valid estimate of $f$ for making a prediction about the outcome of an intervention on $X$ Newey03:NIV,Hartford17:DIV,Singh19:KIV. To handle hidden confounders, we assume access to an instrumental variable(s) $Z$ which satisfies the following assumptions:
The properties of $Z$ imply that $\mathbb{E}[\varepsilon \,|\,Z] = 0$. Taking the expectation of (ref) w.r.t. $Y$ conditioned on $Z$ yields the following integral equation
which is a Fredholm integral equation of the first kind. Recent works in nonparametric IV regression have adopted this perspective horowitz2011applied,Newey03:NIV,Hartford17:DIV,Singh19:KIV despite the fact that solving (ref) directly is an ill-posed problem as it involves inverting linear compact operators Kress89:LIE,Nashed74:Inverse,Newey03:NIV.
To illustrate the role of an instrument, we consider two special cases. When $X$ is perfectly correlated with $Z$, the treatment is uncorrelated with the hidden confounder. In other words, we recover the strong ignorability assumption Rubin74:Causal,Rubin05:PO required for causal inference. When $Z$ is independent of $X$, the instrument is useless as it has no predictive power over treatment so the structural function $f$ is unidentifiable from the data. Therefore, the most interesting cases lie between these two extremes, especially when $X$ and $Z$ are weakly correlated, see, e.g., Bound95:WeakInstrument,Staiger97:WeakAngrist08:Harmless.
Early applications of instrumental variables often assume linear relationships between $Z$ and $X$ as well as $X$ and $Y$ Angrist96:IV,Angrist08:Harmless. When there is a single endogeneous variable and instrument, the structural parameter can be estimated consistently by the instrumental variable (IV) estimator Angrist96:IV. Interestingly, we can obtain this estimate using a two-stage procedure: regress $X$ on $Z$ using ordinary least square (OLS) to calculate the predicted value of $X$ and used that as an explanatory variable in the structural equation to estimate the structural parameter using OLS. When there are multiple instruments, the two-stage least squares (2SLS) estimator is obtained by using all the instruments simultaneously in the first-stage regression. Wooldridge01:Econometric asserts that the 2SLS estimator is the most efficient IV estimator; see, e.g., Wooldridge01:Econometric,Angrist08:Harmless for a detailed exposition.
Recently, several extensions of 2SLS have been proposed to overcome the linearity constraint. The first line of work replaces linear regression by a linear projection onto a set of known basis functions Newey03:NIV,Blundell07:Semi-NIV,horowitz2011applied,Chen12:NIV. Chen18:NIV-Uniform provides a uniform convergence rate of this approach. However, there exists no principled way of choosing the appropriate set of basis functions. The second line of work replaces the first-stage regression by a conditional density estimate of $\mathbb{P}(X|Z)$ Hall05:NIV,Darolles11:NIV. Despite being more flexible, such approaches are known to suffer from the curse of dimensionality Tsybakov08:INE. Other extensions of 2SLS are DeepIV Hartford17:DIV and KernelIV Singh19:KIV algorithms. In Hartford17:DIV, (ref) is solved by first estimating $\mathbb{P}(X|Z)$ with a mixture of deep generative models on which $f$ is learned using another deep neural network. Instead of neural networks, Singh19:KIV proposes to model the first-stage regression using the conditional mean embedding of $\mathbb{P}(X|Z)$ Song10:KCOND,Song2013,Muandet17:KME which is then used in the second-stage kernel ridge regression.
\paragraph{The curse of two-stage methods.} Two-stage procedures have two fundamental issues. First, such procedures violate Vapnik's principle Vapnik98:SLT: “[...] when solving a problem of interest, do not solve a more general problem as an intermediate step [...]”. Specifically, estimating the conditional density Hartford17:DIV or the conditional mean embedding Singh19:KIV via regression in the first stage can be harder than estimating the parameter of interest in the second stage. The first stage is even referred as the “forbidden regression” in econometrics Angrist08:Harmless. On top of that, we usually only observe a single sample from each $\mathbb{P}(X|Z=z)$, which further increases the difficulty of the task. Second, although two-stage procedures are asymptotically consistent, the first-stage estimate creates a finite-sample bias in the second-stage estimate Angrist08:Harmless. This bias can be alleviated through sample splitting Angrist93:Split which is also used in Hartford17:DIV,Singh19:KIV. Thus, two-stage procedures are less sample efficient and could yield biased estimates when run on the smaller datasets common in economics and social sciences.
The generalized method of moments (GMM) framework provides another set of popular approaches for estimating $f$ Hansen82:GMM,Hall05:GMM. Unlike two-stage procedures, GMM-based algorithms find a function $f$ that satisfies the orthogonality condition $\mathbb{E}[\varepsilon | Z]=0$ directly. Specifically, if $g_1,g_2,\ldots,g_m$ are arbitrary real-valued functions, the orthogonality condition implies that $\mathbb{E}[(Y - f(X))g_j(Z)] = 0$ for $j=1,\ldots,m$. The GMM estimate of $f$ can then be obtained by minimizing the quadratic form $\frac{1}{2}\sum_{j=1}^m\psi(f,g_j)^2$ where $\psi(f,g) := \mathbb{E}[(Y - f(X))g(Z)]$. This estimator can be interpreted as a generalization of the 2SLS estimator in the linear setting White82:IV. Recently, extensions of GMM-based methods where both $f$ and $g$ are parameterized by deep neural networks have successfully been used to solve non-linear IV regression Lewis18:AGMM,Bennett19:DeepGMM. In contrast, Muandet20:KCM considers the set of RKHS functions which allow for an analytic formulation of the orthogonality condition.
In this section, we reformulate the integral equation (ref) as an empirical risk minimization problem and present DualIV algorithm.
Let $\ell: \mathbb{R}\times\mathbb{R}\to\mathbb{R}_+$ be a proper, convex, and lower semi-continuous loss function for any value in its first argument.\footnote{The function $f$ is proper if $\text{dom}\, f\neq\emptyset$ and $f(x)>-\infty, \,\forall x\in\ensuremath{\mathcal{X}}$. It is lower (upper) semi-continuous at $x_0\in\ensuremath{\mathcal{X}}$ if for $\varepsilon > 0$ there exists a neighborhood $N(x_0)$ of $x_0$ such that $\varepsilon < (>) f(x)-f(x_0)$ for all $x\in N(x_0)$ Rockafellar70:Convex.} Let $\mathcal{F}$ be an arbitrary class of continuous functions which we assume contains $f$ that fulfills the integral equation (ref). Then, we can formulate (ref) as
where $R(f)$ denotes the expected risk of $f$. To understand how (ref) and (ref) are related, let us consider the squared loss $\ell(y,y') = (y-y')^2$ and define $h(z) := \mathbb{E}_{X|z}[f(X)]$. Then, the solution to (ref) is the minimum mean square error (MMSE) estimator $h^{*}(z) := \mathbb{E}[Y|z]$, which is exactly the LHS of (ref). If there exists no $f\in\mathcal{F}$ for which $h^{*}(z)=\mathbb{E}[Y|z]$, we use $h^{*}(z)$ as the best MMSE approximation.
The key challenge is if $f$ is noncontinuous in $h(z)$, it is not guaranteed to be consistently estimated even if $h(z)$ is estimated correctly Newey03:NIV. We defer further discussion to Section (ref). In addition, it remains cumbersome to solve (ref) directly because of the inner expectation. To circumvent this, we look to similar two-stage problems in stochastic programming Shapiro14:LSP,Dai17:Dual. For example, in Dai17:Dual, the problem of learning from conditional distributions was formulated in a similar fashion to (ref). Moreover, Hsu19 proposes the deconditional mean embedding (DME) which solves the integral equation (ref) by performing a closed-form “inversion” of the conditional mean embedding of $\ensuremath{\mathbb{P}}(X|Z)$ (see Muandet17:KME,Song13:CME for a review). In contrast, we solve the equivalent dual formulation of (ref) instead of (ref).
To derive the dual of (ref), we employ two existing results, interchangeability and Fenchel duality, which we review; see, e.g., Dai17:Dual, Rock98:VA, and Shapiro14:LSP for more details.
Applying the interchangeability and Fenchel duality to (ref) yields the expected loss
where $\mathcal{U}$ is the space of continuous functions over $\mathcal{Y}\times\mathcal{Z}$. Hence, (ref) can be reformulated as
Following Dai17:Dual, we will refer to $u\in\mathcal{U}$ as the dual function. Note that this function depends on only the outcome $Y$ and the instrument $Z$, but not the treatment $X$.
The advantages of our formulation (ref) over (ref) and (ref) are twofold. First, there is no need to estimate $\mathbb{E}_{X|Z}[f(X)]$ or $\ensuremath{\mathbb{P}}(X|Z)$ explicitly. Second, the target function $f$ appears linearly in (ref) which makes it convex in $f$. Since $\ell^{\star}_y$ is also convex, (ref) is concave in the dual function $u$. Hence, (ref) is essentially a convex-concave saddle-point problem for which efficient solvers exist Dai17:Dual.
For the squared loss $\ell(y,y') = (y-y')^2$, we have $\ell^{\star}_y(w) = wy + \frac{1}{2}w^2$ (see Appendix (ref) for the derivation) and the saddle-point problem (ref) reduces to
To solve (ref), one can adopt an SGD-based algorithm developed by Dai17:Dual. Alternatively, we propose in Section (ref) a simple algorithm that can solve (ref) in closed form.
The dual function $u(y,z)$ plays an important role in our framework. To understand its role, we consider the minimization and maximization problems in (ref) separately. For any $f\in\mathcal{F}$, the maximization problem is $\max_{u\in\mathcal{U}} \mathbb{E}_{\mathit{XYZ}}[(f(X) - Y) u(Y,Z)] - \frac{1}{2}\|u(Y, Z)\|_{L^2(\ensuremath{\mathbb{P}}_{\mathit{YZ}})}^2$ where the first term can be viewed, loosely speaking, as a loss function and the second as a regularizer. Intuitively, we are seeking $u^{*}\in\mathcal{U}$ that is least orthogonal to the residual. Given $u^{*}$, the outer minimization problem $\min_{f\in\mathcal{F}} \; \mathbb{E}_{\mathit{XYZ}}\left[(f(X)-Y) u^{*}(Y,Z)\right]$ finds the function $f$ that yields the most orthogonal residual to $u^{*}$. Our procedure clearly differs from previous two-stage methods as the minimization and maximization stages are interdependent.
Examining the formulation in the context of instrumental variable regression, the residual contains the variation that cannot be explained by the current estimate of $f$ due to hidden confounding. We select $u$ that maximally reweights the residuals according to how inconsistent they are w.r.t. the unconfounded joint distribution of $Y$ and $Z$. Given $u$, we then select $f$ that minimizes the inconsistencies between the residuals and $u$. Hence, at the equilibrium, we are left with residuals uncorrelated with $Y$ and $Z$ which can be attributed to noise due to unobserved confounding.
Lastly, we draw a connection between (ref) and GMM. Let $g_1,g_2,\ldots,g_m$ be real-valued functions on $\ensuremath{\mathcal{Y}}\times\ensuremath{\mathcal{Z}}$ and $\psi(f,g) := \mathbb{E}[(Y-f(X))g(Y,Z)]$. When $\mathcal{U} = \text{span}\{g_1,\ldots,g_m\}$, it is not difficult to show that $\max_{u\in\mathcal{U}} \Psi(f,u) = \frac{1}{2}\bm{\psi}^\top\Lambda^{-1}\bm{\psi}$ where $\Lambda := \mathbb{E}_{\mathit{YZ}}[\mathbf{g}\otimes\mathbf{g}]$ with $\mathbf{g} := (g_1(Y,Z),\ldots,g_m(Y,Z))^\top$; see Appendix (ref). That is, minimizing the above over $f$ yields a formulation that strongly resembles the GMM objective, with the dual function $u(Y,Z)$ playing a role similar to that of an instrument. However, we must clarify that $u$ cannot act as an instrument since it depends on $Y$ and thereby violates the exclusion restriction assumption. Liao20:NeuralSEM, using an alternative formulation similar to (ref) and (ref), showed that one can obtain a dual function $u(Z)$ that can act as an instrument. Furthermore, we also note that AGMM Lewis18:AGMM and DeepGMM Bennett19:DeepGMM rely on minimax optimization, similar to (ref), but were formulated based on the GMM framework.
This section provides the conditions for which the true structural function $f^*$ can be identified by the optimum of the saddle-point problem (ref). We lay out the assumptions needed for the optimal dual function $u^{*}$ to be unique and continuous, show that the saddle-point formulation (ref) is equivalent to the problem (ref) under the squared loss and prove that the solution of (ref) given $u^{*}$ is indeed $f^*$.
Following Dai17:Dual, we define the optimal dual function for any pair $(y,z)\in\ensuremath{\mathcal{Y}}\times\ensuremath{\mathcal{Z}}$ as $u^*(y,z) \in \arg\max_{u\in\mathbb{R}} \{ \mathbb{E}_{X|z}[f(X)-y]u - (1/2)u^2 \}$. Since this is an unconstrained quadratic program, $u^*(y,z)$ takes the form $\mathbb{E}_{X|z}[f(X)] - y$. Given Assumption (ref) and the loss function $\ell$ is convex and continuously differentiable, it follows from Dai17:Dual that $u^*$ is unique and continuous.
Next, we shows that if $(f^*,u^*)$ is the saddle-point of (ref), $f^*$ minimizes the original objective (ref). The result follows from plugging $u^* = \mathbb{E}_{X|z}[f(X)]-y$ into the dual loss $\Psi(f,u)$ in (ref); see Appendix (ref) for the detailed proof.
By Proposition (ref) and the convexity of the loss $\ell(y,y')$, we obtain the following result.
By virtue of Theorem (ref), we can identify the true function $f^*$ under relatively weak assumptions. In contrast, previous work usually require stronger assumptions such as the completeness condition Newey03:NIV,Singh19:KIV which specifies that the first-stage conditional expectation $\mathbb{E}_{X|z}[f(X)]$ is injective, or $h(z) = \mathbb{E}_{X|z}[f(X)]$ is a smooth function of $z$ Singh19:KIV,Chen18:NIV-Uniform,Chen12:NIV. Since we do not perform first-stage regression, we only require $\ensuremath{\mathbb{P}}(X|Z)$ is continuous in $Z$ for any value of $X$. The assumption that (ref) is correctly specified, i.e., $f^{*}\in\mathcal{F}$, is standard in the literature Newey03:NIV,horowitz2011applied,Singh19:KIV.
As we can see, the optimal dual function $u^{*}(y,z) = \mathbb{E}_{X|z}[f(X)]-y$ acts as a residual function measuring the discrepancy between $y$ and $\mathbb{E}_{X|z}[f(X)]$ Dai17:Dual. Remarkably, this makes it possible to approximate $R(f)$ in (ref) without computing the expectation $\mathbb{E}_{X|z}[f(X)]$ explicitly. We later exploit this property when performing hyperparameter selection. Moreover, since $\mathbb{E}_{X|z}[f(X)]$ allows $X$ and $Z$ to have a non-linear relationship, $u$ can be non-linear even when the true structural function $f^*$ is linear. This flexibility enables $u$ to accommodate a larger class of functions that maps $Z$ to $X$. Figure (ref) illustrates this given the following generative process:
where $e \sim \mathcal{N}(0,2), Z_1\sim\mathcal{N}(0,2), \epsilon \sim \mathcal{N}(0,0.1)$, and $\eta\sim\mathcal{N}(0,0.1)$. The parameter $\rho$ controls the strength of the instrument w.r.t. hidden confounder $e$. Here, we set $n=300$, $\beta=0.7$, $\hat{\beta}=0.4$, and $\rho=0.2$ where $\hat{\beta}$ is an OLS estimate of $\beta$. Under this model, we have $u^*(y,z) \approx \hat{\beta}(1-\rho)z - y$.
To demonstrate the effectiveness of our framework, we develop a simple kernel-based algorithm using the new formulation (ref). To simplify notation, we denote by $W := (Y,Z)$ a random variable taking value in $\mathcal{W} := \mathcal{Y}\times\mathcal{Z}$. We pick $\ensuremath{\mathcal{F}}$ and $\ensuremath{\mathcal{U}}$ to be reproducing kernel Hilbert spaces (RKHSs) associated with positive definite kernels $k:\ensuremath{\mathcal{X}}\times\ensuremath{\mathcal{X}}\to\mathbb{R}$ and $l:\ensuremath{\mathcal{W}} \times \ensuremath{\mathcal{W}} \to\mathbb{R}$, respectively. Let $\phi: x\mapsto k(x,\cdot)$ and $\varphi: w\mapsto l(w,\cdot)$ be the canonical feature maps Scholkopf01:LKS. We assume both $\ensuremath{\mathcal{F}}$ and $\ensuremath{\mathcal{U}}$ are universal and hence are dense in the space of bounded continuous functions SteChr2008.
Then, for any $f\in\ensuremath{\mathcal{F}}$ and $u\in\ensuremath{\mathcal{U}}$, we can rewrite the objective in (ref) as
where $\mathbf{b} := \mathbb{E}_{\mathit{YZ}}[Y\varphi(Y,Z)]\in\ensuremath{\mathcal{U}}$, $\ensuremath{\mathcal{C}_{\mathit{W}}} := \mathbb{E}_{W}[\varphi(W)\otimes\varphi(W)] \in\ensuremath{\mathcal{U}}\otimes\ensuremath{\mathcal{U}}$ is a covariance operator, and $\ensuremath{\mathcal{C}_{\mathit{WX}}} := \mathbb{E}_{\mathit{WX}}[\varphi(W)\otimes\phi(X)] \in \ensuremath{\mathcal{U}}\otimes\ensuremath{\mathcal{F}}$ is a cross-covariance operator Baker1973,Fukumizu04:DRS (see Appendix (ref)). Since (ref) is quadratic in $u$, we have $\ensuremath{\mathcal{C}_{\mathit{W}}} u^* = \ensuremath{\mathcal{C}_{\mathit{WX}}} f - \mathbf{b}$. Substituting $u^*$ back into (ref) yields
We can view (ref) as a generalized least squares solution in RKHS. Since $\ensuremath{\mathcal{C}_{\mathit{W}}}^{-1}$ and $(\ensuremath{\mathcal{C}_{\mathit{XW}}}\ensuremath{\mathcal{C}_{\mathit{W}}}^{-1}\ensuremath{\mathcal{C}_{\mathit{WX}}})^{-1}$ may not exist in general, we replace them with regularized versions $(\ensuremath{\mathcal{C}_{\mathit{W}}} + \lambda_1\mathcal{I})^{-1}$ and $(\ensuremath{\mathcal{C}_{\mathit{XW}}}\ensuremath{\mathcal{C}_{\mathit{W}}}^{-1}\ensuremath{\mathcal{C}_{\mathit{WX}}} + \lambda_2\mathcal{I})^{-1}$ where $\mathcal{I}$ is the identity operator and $\lambda_1,\lambda_2 > 0$ are regularization parameters.
Given an i.i.d. sample $(x_i,y_i,z_i)_{i=1}^n$ from $\mathbb{P}(X,Y,Z)$, we define $\Phi := [\phi(x_1),\ldots,\phi(x_n)]$, $\Upsilon := [\varphi(y_1,z_1),\ldots,\varphi(y_n,z_n)]$, and $\ensuremath{\mathbf{y}}:=[y_1,\ldots,y_n]^\top$. Then, we can estimate $\mathbf{b}$, $\ensuremath{\mathcal{C}_{\mathit{W}}}$, and $\ensuremath{\mathcal{C}_{\mathit{XW}}}$ with their empirical counterparts $\hat{\mathbf{b}} := n^{-1}\sum_{i=1}^n y_i\varphi(y_i,z_i) = n^{-1}\Upsilon\mathbf{y}$, $\ensuremath{\widehat{\mathcal{C}}_{\mathit{XW}}} := n^{-1}\sum_{i=1}^n\phi(x_i)\otimes\varphi(y_i,z_i) = n^{-1}\Phi\Upsilon^\top$ and $\ensuremath{\widehat{\mathcal{C}}_{\mathit{W}}} = n^{-1}\sum_{i=1}^n\varphi(y_i,z_i)\otimes\varphi(y_i,z_i) = n^{-1}\Upsilon\Upsilon^\top$. We denote the empirical version of (ref) by $\widehat{\Psi}(f,u)$ and the estimate of $f^*$ by $\hat{f}$.
Next, we show that the representer theorem Schoelkopf01:Representer for $\widehat{\Psi}(f,u)$ holds for both $f$ and $u$.
By virtue of Lemma (ref), the solution to (ref) can be expressed as $f(x) = \sum_{i=1}^n\beta_ik(x_i,x)$ where the coefficients $\bm{\beta}$ are given by the following proposition.
Compared to previous work which involved conditional density estimation Hall05:NIV,Darolles11:NIV,Hartford17:DIV and vector-valued regression Singh19:KIV as first-stage regression, estimating the dual function $u$, a real-valued function, is arguably easier. This is especially so when $X$ and $Z$ are high-dimensional.
\paragraph{Hyperparameter selection.} Our estimator depends on two hyper-parameters, $\lambda_1$ and $\lambda_2$. Given a dataset $(x_i,y_i,z_i)_{i=1}^{2n}$ of size $2n$, we provide a simple heuristic to determine the values of $(\lambda_1,\lambda_2)$. Ideally, if we know the optimal dual function $u^*$, we can interpret $u^*(y,z)^2$ as a loss function of $f^*$ at $(y,z)$, as discussed in Section (ref). To this end, we first estimate $\hat{f}$ via Proposition (ref) and $\hat{u}_{\lambda} := (\ensuremath{\widehat{\mathcal{C}}_{\mathit{W}}} + \lambda\mathcal{I})^{-1}(\ensuremath{\widehat{\mathcal{C}}_{\mathit{WX}}}\hat{f} - \hat{\mathbf{b}})$ on the first half of the data $(x_i,y_i,z_i)_{i=1}^n$. Next, the out-of-sample loss of $\hat{f}$ is evaluated on the second half $(x_i,y_i,z_i)_{i=n+1}^{2n}$ by
Note that $\hat{u}_{\lambda} = \Upsilon(\mathbf{L} + n\lambda I)^{-1}(\mathbf{K}\bm{\beta} - \ensuremath{\mathbf{y}}) = \Upsilon\bm{\alpha}$ where $\bm{\alpha} := (\mathbf{L} + n\lambda I)^{-1}(\mathbf{K}\bm{\beta} - \ensuremath{\mathbf{y}})$ and $\mathbf{K},\mathbf{L} \in\mathbb{R}^{n\times n}$ are kernel matrices evaluated on $(x_i,y_i,z_i)_{i=1}^n$. Hence, $\widehat{R}(\hat{f}) \approx \bm{\alpha}^\top\widetilde{\mathbf{L}}\mathbf{1}/n$ where $\widetilde{\mathbf{L}}_{ij}=l((y_i,z_i),(y_j,z_j))$ for $i=1,\ldots,n$ and $j=n+1,\ldots,2n$ and $\mathbf{1}$ is the all-ones column vector. In practice, we fix $\lambda$ in $\hat{u}_{\lambda}$ to a small constant to stabilize the loss (ref) and only optimize $(\lambda_1,\lambda_2)$ that appear in $\bm{\beta}$. Note that this procedure differs from the two-stage causal validation procedures used in Hartford17:DIV,Singh19:KIV. Alternatively, one may choose the hyperparameters by cross-validation with respect to (ref).
Algorithm (ref) outlines the kernelized DualIV method whose consistency is studied in Appendix (ref). We note that above algorithm involves matrix inversions ($\mathcal{O}(n^3)$) which become the primary computational bottlenecks when scaling to large datasets. To improve the scalability of our algorithm, we can leverage the rich literature on large-scale kernel machines such as random Fourier features and Nystr\"om method; see, e.g., NIPS2012_4588 and references therein. Alternatively, we can employ stochastic gradient descent-based (SGD) algorithms similar to those proposed in Dai17:Dual to solve the dual formulation (ref) directly. This would also allow us to employ flexible models such as neural networks to parameterize the function classes $\mathcal{F}$ and $\mathcal{U}$. Recently, Liao20:NeuralSEM has taken an important step in this direction and provided convergence analysis for neural networks under a similar formulation.
In this section, we compare kernelized DualIV\footnote{Our implementation is available at \url{https://github.com/krikamol/DualIV-NeurIPS2020}.} with:
To provide a fair comparison, we adhered to the provided hyperparameter settings. Given the low-dimensional nature of our experiments, we used DeepGMM's settings for low-dimensional scenarios in Bennett19:DeepGMM. We ran 20 simulations of each algorithm for sample sizes of $50$ and $1000$ and calculated the $\log_{10}$ mean squared error and its standard deviations w.r.t. the true function $f$ for 2800 out-of-sample test points.
\paragraph{Demand design.} We consider the same simulation as in Hartford17:DIV,Singh19:KIV: $Y = f(X) + \varepsilon$ where $Y$ is outcome, $X=(P,T,S)$ are inputs, and $Z=(C,T,S)$ are instruments. Specifically, $Y$ is sales, $P$ is price, which is endogeneous, $C$ is a supply cost shifter (instrument), and $(T,S)$ are time of year and customer sentiment acting as exogeneous variables. The aim is to estimate the demand function $f(p,t,s) = 100 + (10+p)s\psi(t)-2p$ where $\psi(t) = 2[(t-5)^4/600 + \exp(-4(t-5)^2)+t/10 -2]$. Training data is sampled according to $S\sim \text{Unif}\{1,\ldots,7\}$, $T\sim\text{Unif}[0,10]$, $(C,V)\sim\mathcal{N}(\mathbf{0},I_2)$, $\varepsilon \sim\mathcal{N}(\rho V, 1-\rho^2)$, and $P = 25 + (C+3)\psi(T)+V$. The parameter $\rho\in\{0.1,0.25,0.5,0.75,0.9\}$ controls the extent to which price $P$ is confounded with outcome $Y$ by supply-side market forces. In our notation, $X=(P,T,S)$, $Z=(C,T,S)$ and $W=(Y,Z)=(Y,C,T,S)$.
For DualIV, we used the Gaussian RBF kernel for both $k$ and $l$. In the experiments, the kernels on $X$ and $Z$ are product kernels, i.e., $k(x_i,x_j) = k_p(p_i,p_j)k_t(t_i,t_j)k_s(s_i,s_j)$ and $k(z_i,z_j) = k_c(c_i,c_j)k_t(t_i,t_j)k_s(s_i,s_j)$, and $l([y_i, z_i],[y_j, z_j])=\exp([y_i-y_j, z_i - z_j]^\top V_{yz}^{-1}[y_i-y_j, z_i - z_j])$ where $V_{yz}$ is a symmetric bandwidth matrix. The values of all bandwidth parameters are determined via the median heuristic. We choose $(\lambda_1, \lambda_2)$ from $\{10^{-10},10^{-9},\ldots,10^{-1}\}$ via cross-validation. Once $(\lambda_1, \lambda_2)$ is chosen, we refit $\hat{f}$ on the entire training set.
Table (ref) reports the results of different methods evaluated on the test data. First, we observe that 2SLS achieves the largest MSE in both regimes as expected because the linearity assumption is violated here. Second, in the small sample size regime, DeepIV achieves relatively larger MSE than the other non-linear methods. KernelIV, DeepGMM, and DualIV, on the other hand, have comparable performance, with DeepGMM having the lowest MSE. However, we note that the results attained by DeepGMM were unstable out of the box and we had to reduce the variance of the initialization of the neural networks to $0.1$ to obtain some degree of stability which is reflected in the standard deviations. We can fully attribute this variability to initialization as DeepGMM's default batch size of 1024 is larger than that of both training datasets so there is no sampling variability in the optimization process. This suggests that DeepGMM, like DeepIV, is relatively brittle compared to kernel-based methods in the small sample size regime. Furthermore, DeepGMM comes with an extensive hyperparameter selection process, which highlights its need for fine-tuning. Last but not least, DualIV is competitive to KernelIV across $\rho$ with slightly smaller MSE, which lends weight to our hypothesis that estimating the real-valued dual function is easier than vector-valued regression.
In the medium sample size regime, we observe that performance of DeepIV is in the same ballpark as the rest of the non-linear IV regression methods and the variance of DeepGMM is reduced, albeit still highest among the non-linear methods. The results of DualIV, KernelIV and DeepGMM are almost indistinguishable with DualIV having an edge as $\rho$ increases. This could mean accounting for both $Y$ and $Z$ is perhaps slightly more effective than $Z$ alone in the presence of greater confounding.
This paper proposes a general framework for non-linear IV regression called DualIV. Unlike previous work, DualIV does not require the first-stage regression which is the critical bottleneck of modern two-stage procedures. By exploiting tools in stochastic programming, we were able to reformulate the two-stage problem as the convex-concave saddle-point problem which is relatively simpler to solve. Instead of first-stage regression, DualIV requires the dual function $u(y,z)$ to be estimated, which is arguably easier than first-stage regression, especially when the instruments and treatments are high-dimensional. We demonstrate the validity of our framework with a kernel-based algorithm. Results show the competitiveness of our algorithm with respect to existing ones. Finally, potential directions for future work include
This work provides a new framework for non-linear instrumental variable regression which allows one to perform causal analysis under the presence of unobserved confounders. This could have a profound impact in other fields such as economics, social science, and epidemiology, among others. Understanding the role of instruments in the context of learning theory may also pave the way towards creating more robust and trustworthy machine learning algorithms that are capable of surviving in the world full of hidden biases.