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.
63,723 characters · 13 sections · 98 citation commands
Kernel Instrumental Variable Regression
Instrumental variable regression is a method in causal statistics for estimating the counterfactual effect of input $X$ on output $Y$ using observational data stock2003retrospectives. If measurements of $(X,Y)$ are confounded, the causal relationship--also called the structural relationship--can nonetheless be identified if an instrumental variable $Z$ is available, which is independent of $Y$ conditional on $X$ and the unmeasured confounder. Intuitively, $Z$ only influences $Y$ via $X$, identifying the counterfactual relationship of interest.
Economists and epidemiologists use instrumental variables to overcome issues of strategic interaction, imperfect compliance, and selection bias. The original application is demand estimation: supply cost shifters ($Z$) only influence sales ($Y$) via price ($X$), thereby identifying counterfactual demand even though prices reflect both supply and demand market forces wright1928tariff,blundell2012measuring. Randomized assignment of a drug ($Z$) only influences patient health ($Y$) via actual consumption of the drug ($X$), identifying the counterfactual effect of the drug even in the scenario of imperfect compliance angrist1996identification. Draft lottery number ($Z$) only influences lifetime earnings ($Y$) via military service ($X$), identifying the counterfactual effect of military service on earnings despite selection bias in enlistment angrist1990lifetime.
The two-stage least squares algorithm (2SLS), widely used in economics, simplifies the IV estimation problem by assuming linear relationships: in stage 1, perform linear regression to obtain the conditional means $\bar{x}(z):=\mathbb{E}_{X|Z=z}(X)$; in stage 2, linearly regress outputs $Y$ on these conditional means. 2SLS works well when the underlying assumptions hold. In practice, the relation between $Y$ and $X$ may not be linear, nor may be the relation between $X$ and $Z$.
In the present work, we introduce kernel instrumental variable regression (KIV), an easily implemented nonlinear generalization of 2SLS (Sections (ref) and (ref)).\footnote{Code: \url{https://github.com/r4hu1-5in9h/KIV}} In stage 1 we learn a conditional mean embedding, which is the conditional expectation $\mu(z):=\mathbb{E}_{X|Z=z}\psi(X)$ of features $\psi$ which map $X$ to a reproducing kernel Hilbert space (RKHS) song2009hilbert. For a sufficiently rich RKHS, called a characteristic RKHS, the mean embedding of a random variable is injective sriperumbudur2010relation. It follows that the conditional mean embedding characterizes the full distribution of $X$ conditioned on $Z$, and not just the conditional mean. We then implement stage 2 via kernel ridge regression of outputs $Y$ on these conditional mean embeddings, following the two-stage distribution regression approach described by szabo2015two,szabo2016learning. As in our work, the inputs for szabo2015two,szabo2016learning are distribution embeddings. Unlike our case, the earlier work uses unconditional embeddings computed from independent samples.
As a key contribution of our work, we provide consistency guarantees for the KIV algorithm for an increasing number of training samples in stages 1 and 2 (Section (ref)). To establish stage 1 convergence, we note that the conditional mean embedding song2009hilbert is the solution to a regression problem grunewalder2012conditional,grunewalder2012modelling,grunewalder2013smooth, and thus equivalent to kernel dependency estimation ciliberto2016consistent,cortes2005general. We prove that the kernel estimator of the conditional mean embedding (equivalently, the conditional expectation operator) converges in RKHS-norm, generalizing classic results by smale2005shannon,smale2007learning. We allow the conditional mean embedding RKHS to be infinite-dimensional, which presents specific challenges that we carefully address in our analysis. We also discuss previous approaches to establishing consistency in both finite-dimensional grunewalder2012modelling and infinite-dimensional song2009hilbert,song2010nonparametric,fukumizu2013kernel,hefny2015supervised,ciliberto2016consistent settings.
We embed the stage 1 rates into stage 2 to get end-to-end guarantees for the two-stage procedure, adapting caponnetto2007optimal,szabo2015two,szabo2016learning. In particular, we provide a ratio of stage 1 to stage 2 samples required for minimax optimal rates in the second stage, where the ratio depends on the difficulty of each stage. We anticipate that these proof strategies will apply generally in two-stage regression settings.
Several approaches have been proposed to generalize 2SLS to the nonlinear setting, which we will compare in our experiments (Section (ref)). A first generalization is via basis function approximation newey2003instrumental, an approach called sieve IV, with uniform convergence rates in chen2018optimal. The challenge in chen2018optimal is how to define an appropriate finite dictionary of basis functions. In a second approach, carrasco2007linear,darolles2011nonparametric implement stage 1 by computing the conditional distribution of the input $X$ given the instrument $Z$ using a ratio of Nadaraya-Watson density estimates. Stage 2 is then ridge regression in the space of square integrable functions. The overall algorithm has a finite sample consistency guarantee, assuming smoothness of the $(X,Z)$ joint density in stage 1 and the regression in stage 2 darolles2011nonparametric. Unlike our bound, darolles2011nonparametric make no claim about the optimality of the result. Importantly, stage 1 requires the solution of a statistically challenging problem: conditional density estimation. Moreover, analysis assumes the same number of training samples used in both stages. We will discuss this bound in more detail in Appendix (ref) (we suggest that the reader first cover Section (ref)).
Our work also relates to kernel and IV approaches to learning dynamical systems, known in machine learning as predictive state representation models (PSRs) boots2013hilbert,hefny2015supervised,downey2017predictive and in econometrics as panel data models anderson1981estimation,arellano1991some. In this setting, predictive states (expected future features given history) are updated in light of new observations. The calculation of the predictive states corresponds to stage 1 regression, and the states are updated via stage 2 regression. In the kernel case, the predictive states are expressed as conditional mean embeddings boots2013hilbert, as in our setting. Performance of the kernel PSR method is guaranteed by a finite sample bound hefny2015supervised, however this bound is not minimax optimal. Whereas hefny2015supervised assume an equal number of training samples in stages 1 and 2, we find that unequal numbers of training samples matter for minimax optimality. More importantly, the bound makes strong smoothness assumptions on the inputs to the stage 1 and stage 2 regression functions, rather than assuming smoothness of the regression functions as we do. We show that the smoothness assumptions on the inputs made in hefny2015supervised do not hold in our setting, and we obtain stronger end-to-end bounds under more realistic conditions. We discuss the PSR bound in more detail in Appendix (ref).
Yet another recent approach is deep IV, which uses neural networks in both stages and permits learning even for complex high-dimensional data such as images hartford2017deep. Like darolles2011nonparametric, hartford2017deep implement stage 1 by estimating a conditional density. Unlike darolles2011nonparametric, hartford2017deep use a mixture density network Bishop06, i.e. a mixture model parametrized by a neural network on the instrument $Z$. Stage 2 is neural network regression, trained using stochastic gradient descent (SGD). This presents a challenge: each step of SGD requires expectations using the stage 1 model, which are computed by drawing samples and averaging. An unbiased gradient estimate requires two independent sets of samples from the stage 1 model hartford2017deep, though a single set of samples may be used if an upper bound on the loss is optimized hartford2017deep. By contrast, our stage 1 outputs--conditional mean embeddings--have a closed form solution and exhibit lower variance than sample averaging from a conditional density model. No theoretical guarantee on the consistency of the neural network approach has been provided.
In the econometrics literature, a few key assumptions make learning a nonparametric IV model tractable. These include the completeness condition newey2003instrumental: the structural relationship between $X$ and $Y$ can be identified only if the stage 1 conditional expectation is injective. Subsequent works impose additional stability and link assumptions blundell2007semi,chen2012estimation,chen2018optimal: the conditional expectation of a function of $X$ given $Z$ is a smooth function of $Z$. We adapt these assumptions to our setting, replacing the completeness condition with the characteristic property sriperumbudur2010relation, and replacing the stability and link assumptions with the concept of prior smale2007learning,caponnetto2007optimal. We describe the characteristic and prior assumptions in more detail below.
Extensive use of IV estimation in applied economic research has revealed a common pitfall: weak instrumental variables. A weak instrument satisfies Hypothesis 1 below, but the relationship between a weak instrument $Z$ and input $X$ is negligible; $Z$ is essentially irrelevant. In this case, IV estimation becomes highly erratic bound1995problems. In staiger1997instrumental, the authors formalize this phenomenon with local analysis. See murray2006avoiding,stock2002survey for practical and theoretical overviews, respectively. We recommend that practitioners resist the temptation to use many weak instruments, and instead use few strong instruments such as those described in the introduction.
Finally, our analysis connects early work on the RKHS with recent developments in the RKHS literature. In nashed1974generalized, the authors introduce the RKHS to solve known, ill-posed functional equations. In the present work, we introduce the RKHS to estimate the solution to an uncertain, ill-posed functional equation. In this sense, casting the IV problem in an RKHS framework is not only natural; it is in the original spirit of RKHS methods. For a comprehensive review of existing work and recent advances in kernel mean embedding research, we recommend muandet2017kernel,gretton2018notes.
Instrumental variable: We begin by introducing our causal assumption about the instrument. This prior knowledge, described informally in the introduction, allows us to recover the counterfactual effect of $X$ on $Y$. Let $(\mathcal{X},\mathcal{B}_{\mathcal{X}})$, $(\mathcal{Y},\mathcal{B}_{\mathcal{Y}})$, and $(\mathcal{Z},\mathcal{B}_{\mathcal{Z}})$ be measurable spaces. Let $(X,Y,Z)$ be a random variable on $\mathcal{X}\times \mathcal{Y}\times \mathcal{Z}$ with distribution $\rho$.
We call $h$ the structural function of interest. The error term $e$ is unmeasured, confounding noise. Hypothesis (ref).1, known as the exclusion restriction, was introduced by newey2003instrumental to the nonparametric IV literature for its tractability. Other hypotheses are possible, although a very different approach is then needed imbens2009identification. Hypothesis (ref).2, known as the relevance condition, ensures that $Z$ is actually informative. In Appendix (ref), we compare Hypothesis (ref) with alternative formulations of the IV assumption.
We make three observations. First, if $X=Z$ then Hypothesis (ref) reduces to the standard regression assumption of unconfounded inputs, and $h(X)=\mathbb{E}[Y|X]$; if $X=Z$ then prediction and counterfactual prediction coincide. The IV model is a framework that allows for causal inference in a more general variety of contexts, namely when $h(X)\neq \mathbb{E}[Y|X]$ so that prediction and counterfactual prediction are different learning problems. Second, Hypothesis (ref) will permit identification of $h$ even if inputs are confounded, i.e. $X\cancel\raisebox{0.05em}{\rotatebox[origin=c]{90}{$\models$}} e$. Third, this model includes the scenario in which the analyst has a combination of confounded and unconfounded inputs. For example, in demand estimation there may be confounded price $P$, unconfounded characteristics $W$, and supply cost shifter $C$ that instruments for price. Then $X=(P,W)$, $Z=(C,W)$, and the analysis remains the same.
Hypothesis (ref) provides the operator equation $ \mathbb{E}[Y|Z]=\mathbb{E}_{X|Z}h(X) $ newey2003instrumental. In the language of 2SLS, the LHS is the reduced form, while the RHS is a composition of stage 1 linear compact operator $\mathbb{E}_{X|Z}$ and stage 2 structural function $h$. In the language of functional analysis, the operator equation is a Fredholm integral equation of the first kind nashed1974generalized,kress1989linear,newey2003instrumental,florens2003inverse. Solving this operator equation for $h$ involves inverting a linear compact operator with infinite-dimensional domain; it is an ill-posed problem kress1989linear. To recover a well-posed problem, we impose smoothness and Tikhonov regularization.
RKHS model: We next introduce our RKHS model. Let $k_{\mathcal{X}}:\mathcal{X}\times \mathcal{X}\rightarrow\mathbb{R}$ and $k_{\mathcal{Z}}:\mathcal{Z}\times \mathcal{Z}\rightarrow\mathbb{R}$ be measurable positive definite kernels corresponding to scalar-valued RKHSs $\mathcal{H}_{\mathcal{X}}$ and $\mathcal{H}_{\mathcal{Z}}$. Denote the feature maps
Define the conditional expectation operator $E:\mathcal{H}_{\mathcal{X}}\rightarrow \mathcal{H}_{\mathcal{Z}}$ such that $ [Eh](z)=\mathbb{E}_{X|Z=z}h(X) $. $E$ is the natural object of interest for stage 1. We define and analyze an estimator for $E$ directly. The conditional expectation operator $E$ conveys exactly the same information as another object popular in the kernel methods literature, the conditional mean embedding $\mu:\mathcal{Z}\rightarrow \mathcal{H}_{\mathcal{X}}$ defined by $ \mu(z)=\mathbb{E}_{X|Z=z}\psi(X) $ song2009hilbert. Indeed, $\mu(z)=E^*\phi(z)$ where $E^*:\mathcal{H}_{\mathcal{Z}}\rightarrow\mathcal{H}_{\mathcal{X}}$ is the adjoint of $E$. Analogously, in 2SLS $\bar{x}(z)=\pi' z$ for stage 1 linear regression parameter $\pi$.
The structural function $h:\mathcal{X}\rightarrow \mathcal{Y}$ in Hypothesis (ref) is the natural object of interest for stage 2. For theoretical purposes, it is convenient to estimate $h$ indirectly. The structural function $h$ conveys exactly the same information as an object we call the structural operator $H:\mathcal{H}_{\mathcal{X}}\rightarrow\mathcal{Y}$. Indeed, $h(x)=H\psi(x)$. Analogously, in 2SLS $h(x)=\beta' x$ for structural parameter $\beta$. We define and analyze an estimator for $H$, which in turn implies an estimator for $h$. Figure (ref) summarizes the relationships among equivalent stage 1 objects $(E,\mu)$ and equivalent stage 2 objects $(H,h)$.
Our RKHS model for the IV problem is of the same form as the model in nashed1974convergence,nashed1974generalized, nashed1974regularization for general operator equations. We begin by choosing RKHSs for the structural function $h$ and the reduced form $\mathbb{E}[Y|Z]$, then construct a tensor-product RKHS for the conditional expectation operator $E$. Our model differs from the RKHS model proposed by carrasco2007linear,darolles2011nonparametric, which directly learns the conditional expectation operator $E$ via Nadaraya-Watson density estimation. The RKHSs of engl1996regularization,carrasco2007linear,darolles2011nonparametric for the structural function $h$ and the reduced form $\mathbb{E}[Y|Z]$ are defined from the right and left singular functions of $E$, respectively. They appear in the consistency argument, but not in the ridge penalty.
2SLS consists of two stages that can be estimated separately. Sample splitting in this context means estimating stage 1 with $n$ randomly chosen observations and estimating stage 2 with the remaining $m$ observations. Sample splitting alleviates the finite sample bias of 2SLS when instrument $Z$ weakly influences input $X$ angrist1995split. It is the natural approach when an analyst does not have access to a single data set with $n+m$ observations of $(X,Y,Z)$ but rather two data sets: $n$ observations of $(X,Z)$, and $m$ observations of $(Y,Z)$. We employ sample splitting in KIV, with an efficient ratio of $(n,m)$ given in Theorem (ref). In our presentation of the general two-stage learning problem, we denote stage 1 observations by $(x_i,z_i)$ and stage 2 observations by $(\tilde{y}_i,\tilde{z}_i)$.
We transform the problem of learning $E$ into a vector-valued kernel ridge regression following grunewalder2012conditional,grunewalder2013smooth,ciliberto2016consistent, where the hypothesis space is the vector-valued RKHS $\mathcal{H}_{\Gamma}$ of operators mapping $\mathcal{H}_{\mathcal{X}}$ to $\mathcal{H}_{\mathcal{Z}}$. In Appendix (ref), we review the theory of vector-valued RKHSs as it relates to scalar-valued RKHSs and tensor product spaces. The key result is that the tensor product space of $\mathcal{H}_{\mathcal{X}}$ and $\mathcal{H}_{\mathcal{Z}}$ is isomorphic to $\mathcal{L}_2(\mathcal{H}_{\mathcal{X}},\mathcal{H}_{\mathcal{Z}})$, the space of Hilbert-Schmidt operators from $\mathcal{H}_{\mathcal{X}}$ to $\mathcal{H}_{\mathcal{Z}}$. If we choose the vector-valued kernel $\Gamma$ with feature map $(x,z)\mapsto [\phi(z)\otimes \psi(x)](\cdot)=\phi(z)\langle \psi(x),\cdot \rangle_{\mathcal{H}_{\mathcal{X}}}$, then $\mathcal{H}_{\Gamma}= \mathcal{L}_2(\mathcal{H}_{\mathcal{X}},\mathcal{H}_{\mathcal{Z}})$ and it shares the same norm.
We now state the objective for optimizing $E\in \mathcal{H}_{\Gamma}$. The optimal $E$ minimizes the expected discrepancy
Both grunewalder2013smooth and ciliberto2016consistent refer to $\mathcal{E}_1$ as the surrogate risk. As shown in grunewalder2012conditional and grunewalder2013smooth, the surrogate risk upper bounds the natural risk for the conditional expectation, where the bound becomes tight when $\mathbb{E}_{X|Z=(\cdot)}f(X)\in \mathcal{H}_{\mathcal{Z}},\:\forall f\in\mathcal{H}_{\mathcal{X}}$. Formally, the target operator is the constrained solution $ E_{\mathcal{H}_{\Gamma}}=\operatorname*{\arg\!\min}_{E\in\mathcal{H}_{\Gamma}}\mathcal{E}_1(E) $. We will assume $E_{\rho}\in \mathcal{H}_{\Gamma}$ so that $E_{\rho}=E_{\mathcal{H}_{\Gamma}}$.
Next we impose Tikhonov regularization. The regularized target operator and its empirical analogue are given by
Our construction of a vector-valued RKHS $\mathcal{H}_{\Gamma}$ for the conditional expectation operator $E$ permits us to estimate stage 1 by kernel ridge regression. The stage 1 estimator of KIV is at once novel in the nonparametric IV literature and fundamentally similar to 2SLS. Basis function approximation newey2003instrumental,chen2018optimal is perhaps the closest prior IV approach, but we use infinite dictionaries of basis functions $\psi$ and $\phi$. Compared to density estimation carrasco2007linear, darolles2011nonparametric,hartford2017deep, kernel ridge regression is an easier problem.
Alternative stage 1 estimators in the literature estimate the singular system of $E$ to ensure that the adjoint of the estimator equals the estimator of the adjoint. These estimators differ in how they estimate the singular system: empirical distribution darolles2011nonparametric, Nadaraya-Watson density darolles2004kernel, or B-spline wavelets chen1997shape. The KIV stage 1 estimator has the desired property by construction; $(E_{\lambda}^n)^*=(E^*)_{\lambda}^n$. See Appendix (ref) for details.
Next, we transform the problem of learning $h$ into a scalar-valued kernel ridge regression that respects the IV problem structure. In Proposition (ref) of Appendix (ref), we show that under Hypothesis (ref) below, $$ \mathbb{E}_{X|Z=z}h(X)=[Eh](z)=\langle h,\mu(z)\rangle_{\mathcal{H}_{\mathcal{X}}}=H\mu(z) $$ where $h\in\mathcal{H}_{\mathcal{X}}$, a scalar-valued RKHS; $E\in \mathcal{H}_{\Gamma}$, the vector-valued RKHS described above; $\mu\in\mathcal{H}_{\Xi}$, a vector-valued RKHS isometrically isomorphic to $\mathcal{H}_{\Gamma}$; and $H\in \mathcal{H}_{\Omega}$, a scalar-valued RKHS isometrically isomorphic to $\mathcal{H}_{\mathcal{X}}$. It is helpful to think of $\mu(z)$ as the embedding into $\mathcal{H}_{\mathcal{X}}$ of a distribution on $\mathcal{X}$ indexed by the conditioned value $z$. When $k_\mathcal{X}$ is characteristic, $\mu(z)$ uniquely embeds the conditional distribution, and $H$ is identified. The kernel $\Omega$ satisfies $k_{\mathcal{X}}(x,x')=\Omega(\psi(x),\psi(x'))$. This expression establishes the formal connection between our model and szabo2015two,szabo2016learning. The choice of $\Omega$ may be more general; for nonlinear examples see szabo2016learning.
We now state the objective for optimizing $H\in \mathcal{H}_{\Omega}$. Hypothesis (ref) provides the operator equation, which may be rewritten as the regression equation $$ Y=\mathbb{E}_{X|Z}h(X)+e_Z=H\mu(Z)+e_Z,\quad \mathbb{E}[e_Z|Z]=0 $$ The unconstrained solution is
The target operator is the constrained solution $ H_{\mathcal{H}_{\Omega}}=\operatorname*{\arg\!\min}_{H\in\mathcal{H}_{\Omega}}\mathcal{E}(H) $. We will assume $H_{\rho}\in \mathcal{H}_{\Omega}$ so that $H_{\rho}=H_{\mathcal{H}_{\Omega}}$. With regularization,
The essence of the IV problem is this: we do not directly observe the conditional expectation operator $E$ (or equivalently the conditional mean embedding $\mu$) that appears in the stage 2 objective. Rather, we approximate it using the estimate from stage 1. Thus our KIV estimator is $\hat{h}^{m}_{\xi}=\hat{H}^{m}_{\xi}\psi$ where
and $\mu^n_{\lambda}=(E_{\lambda}^n)^*\phi$. The transition from $H_{\rho}$ to $H^{m}_{\xi}$ represents the fact that we only have $m$ samples. The transition from $H^{m}_{\xi}$ to $\hat{H}^{m}_{\xi}$ represents the fact that we must learn not only the structural operator $H$ but also the conditional expectation operator $E$. In this sense, the IV problem is more complex than the estimation problem considered by nashed1974convergence,nashed1974regularization in which $E$ is known.
We obtain a closed form expression for the KIV estimator. The apparatus introduced above is required for analysis of consistency and convergence rate. More subtly, our RKHS construction allows us to write kernel ridge regression estimators for both stage 1 and stage 2, unlike previous work. Because KIV consists of repeated kernel ridge regressions, it benefits from repeated applications of the representer theorem wahba1990spline,scholkopf2001generalized. Consequently, we have a shortcut for obtaining KIV's closed form; see Appendix (ref) for the full derivation.
Theorems (ref) and (ref) below theoretically determine efficient rates for the stage 1 regularization parameter $\lambda$ and stage 2 regularization parameter $\xi$, respectively. In Appendix (ref), we provide a validation procedure to empirically determine values for $(\lambda,\xi)$.
Integral operators: We use integral operator notation from the kernel methods literature, adapted to the conditional expectation operator learning problem. We denote by $L^2(\mathcal{Z},\rho_{\mathcal{Z}})$ the space of square integrable functions from $\mathcal{Z}$ to $\mathcal{Y}$ with respect to measure $\rho_{\mathcal{Z}}$, where $\rho_{\mathcal{Z}}$ is the restriction of $\rho$ to $\mathcal{Z}$.
$T_1=S_1\circ S_1^*$ is the uncentered covariance operator of fukumizu2004dimensionality. In Appendix (ref), we prove that $T_1$ exists and has finite trace even when $\mathcal{H}_{\mathcal{X}}$ and $\mathcal{H}_{\mathcal{Z}}$ are infinite-dimensional. In Appendix (ref), we compare $T_1$ with other covariance operators in the kernel methods literature.
Assumptions: We place assumptions on the original spaces $\mathcal{X}$ and $\mathcal{Z}$, the scalar-valued RKHSs $\mathcal{H}_{\mathcal{X}}$ and $\mathcal{H}_{\mathcal{Z}}$, and the probability distribution $\rho(x,z)$. We maintain these assumptions throughout the paper. Importantly, we assume that the vector-valued RKHS regression is correctly specified: the true conditional expectation operator $E_{\rho}$ lives in the vector-valued RKHS $\mathcal{H}_{\Gamma}$. In further research, we will relax this assumption.
Hypothesis 3.3 specializes the completeness condition of newey2003instrumental. Hypotheses 2-4 are sufficient to bound the sampling error of the regularized estimator $E_{\lambda}^n$. Bounding the approximation error requires a further assumption on the smoothness of the distribution $\rho(x,z)$. We assume $\rho(x,z)$ belongs to a class of distributions parametrized by $(\zeta_1,c_1)$, as generalized from smale2007learning to the space $\mathcal{H}_\Gamma$.
We use composition symbol $\circ$ to emphasize that $G_1:\mathcal{H}_{\mathcal{X}}\rightarrow\mathcal{H}_{\mathcal{Z}}$ and $T_1:\mathcal{H}_{\mathcal{Z}}\rightarrow \mathcal{H}_{\mathcal{Z}}$. We define the power of operator $T_1$ with respect to its eigendecomposition; see Appendix (ref) for formal justification. Larger $c_1$ corresponds to a smoother conditional expectation operator $E_{\rho}$. Proposition (ref) in Appendix (ref) shows $E_{\rho}^*\phi(z)=\mu(z)$, so Hypothesis 5 is an indirect smoothness condition on the conditional mean embedding $\mu$.
Estimation and convergence: The estimator has a closed form solution, as noted in grunewalder2012conditional and grunewalder2012modelling; ciliberto2016consistent use it in the first stage of the structured prediction problem. We present the closed form solution in notation similar to caponnetto2007optimal in order to elucidate how the estimator simply generalizes linear regression. This connection foreshadows our proof technique.
We prove an original, finite sample bound on the RKHS-norm distance of the estimator $E^n_{\lambda}$ from its target $E_{\rho}$. The proof is in Appendix (ref).
The efficient rate of $\lambda$ is $n^{\frac{-1}{c_1+1}}$. Note that the convergence rate of $E^n_{\lambda}$ is calibrated by $c_1$, which measures the smoothness of the conditional expectation operator $E_{\rho}$.
Integral operators: We use integral operator notation from the kernel methods literature, adapted to the structural operator learning problem. We denote by $L^2(\mathcal{H}_{\mathcal{X}},\rho_{\mathcal{H}_{\mathcal{X}}})$ the space of square integrable functions from $\mathcal{H}_{\mathcal{X}}$ to $\mathcal{Y}$ with respect to measure $\rho_{\mathcal{H}_{\mathcal{X}}}$, where $\rho_{\mathcal{H}_{\mathcal{X}}}$ is the extension of $\rho$ to $\mathcal{H}_{\mathcal{X}}$ steinwart2008support. Note that we present stage 2 analysis for general output space $\mathcal{Y}$ as in szabo2015two,szabo2016learning, though in practice we only consider $\mathcal{Y}\subset \mathbb{R}$ to simplify our two-stage RKHS model.
where $\Omega_{\mu(z)}: \mathcal{Y}\rightarrow\mathcal{H}_{\Omega}$ defined by $y\mapsto \Omega(\cdot,\mu(z))y$ is the point evaluator of micchelli2005learning,carmeli2006vector. Finally define $T_{\mu(z)}=\Omega_{\mu(z)}\circ \Omega^*_{\mu(z)}$ and covariance operator $T=S\circ S^*$.
Assumptions: We place assumptions on the original space $\mathcal{Y}$, the scalar-valued RKHS $\mathcal{H}_{\Omega}$, and the probability distribution $\rho$. Importantly, we assume that the scalar-valued RKHS regression is correctly specified: the true structural operator $H_{\rho}$ lives in the scalar-valued RKHS $\mathcal{H}_{\Omega}$.
Larger $\iota$ is interpretable as smoother kernel $\Omega$.
The convergence rate from stage 1 together with Hypotheses 6-8 are sufficient to bound the excess error of the regularized estimator $\hat{H}_{\xi}^{m}$ in terms of familiar objects in the kernel methods literature, namely the residual, reconstruction error, and effective dimension. We further assume $\rho$ belongs to a stage 2 prior to simplify these bounds. In particular, we assume $\rho$ belongs to a class of distributions parametrized by $(\zeta,b,c)$ as defined originally in caponnetto2007optimal, restated below.
We define the power of operator $T$ with respect to its eigendecomposition; see Appendix (ref) for formal justification. The latter condition is interpretable as polynomial decay of eigenvalues: $\lambda_k=\Theta(k^{-b})$. Larger $b$ means faster decay of eigenvalues of the covariance operator $T$ and hence smaller effective input dimension. Larger $c$ corresponds to a smoother structural operator $H_{\rho}$ szabo2016learning.
Estimation and convergence: The estimator has a closed form solution, as shown by szabo2015two,szabo2016learning in the second stage of the distribution regression problem. We present the solution in notation similar to caponnetto2007optimal to elucidate how the stage 1 and stage 2 estimators have the same structure.
We now present this paper's main theorem. In Appendix (ref), we provide a finite sample bound on the excess error of the estimator $\hat{H}^{m}_{\xi}$ with respect to its target $H_{\rho}$. Adapting arguments by szabo2016learning, we demonstrate that KIV is able to achieve the minimax optimal single-stage rate derived by caponnetto2007optimal. In other words, our two-stage estimator is able to learn the causal relationship with confounded data equally well as single-stage RKHS regression is able to learn the causal relationship with unconfounded data.
At $a=\frac{b(c+1)}{bc+1}<2$, the convergence rate $m^{-\frac{bc}{bc+1}}$ is minimax optimal while requiring the fewest observations szabo2016learning. This statistically efficient rate is calibrated by $b$, the effective input dimension, as well as $c$, the smoothness of structural operator $H_{\rho}$ caponnetto2007optimal. The efficient ratio between stage 1 and stage 2 samples is $n=m^{\frac{b(c+1)}{bc+1}\cdot \frac{(c_1+1)}{\iota(c_1-1)}}$, implying $n>m$. As far as we know, asymmetric sample splitting is a novel prescription in the IV literature; previous analyses assume $n=m$ angrist1995split,hefny2015supervised.
We compare the empirical performance of KIV (\verb|KernelIV|) to four leading competitors: standard kernel ridge regression (\verb|KernelReg|) saunders1998ridge, Nadaraya-Watson IV (\verb|SmoothIV|) carrasco2007linear,darolles2011nonparametric, sieve IV (\verb|SieveIV|) newey2003instrumental,chen2018optimal, and deep IV (\verb|DeepIV|) hartford2017deep. To improve the performance of sieve IV, we impose Tikhonov regularization in both stages with KIV's tuning procedure. This adaptation exceeds the theoretical justification provided by chen2018optimal. However, it is justified by our analysis insofar as sieve IV is a special case of KIV: set feature maps $\psi,\phi$ equal to the sieve bases.
We implement each estimator on three designs. The linear design chen2018optimal involves learning counterfactual function $h(x)=4x-2$, given confounded observations of continuous variables $(X,Y)$ as well as continuous instrument $Z$. The sigmoid design chen2018optimal involves learning counterfactual function $h(x)=\ln (|16x-8|+1)\cdot sgn(x-0.5)$ under the same regime. The demand design hartford2017deep involves learning demand function $h(p,t,s)=100+(10+p)\cdot s\cdot \psi(t)-2p$ where $ \psi(t)$ is the complex nonlinear function in Figure (ref). An observation consists of $(Y,P,T,S,C)$ where $Y$ is sales, $P$ is price, $T$ is time of year, $S$ is customer sentiment (a discrete variable), and $C$ is a supply cost shifter. The parameter $\rho\in\{0.9,0.75,0.5,0.25,0.1\}$ calibrates the extent to which price $P$ is confounded by supply-side market forces. In KIV notation, inputs are $X=(P,T,S)$ and instruments are $Z=(C,T,S)$.
For each algorithm, design, and sample size, we implement 40 simulations and calculate MSE with respect to the true structural function $h$. Figures (ref), (ref), and (ref) visualize results. In the sigmoid design, \verb|KernelIV| performs best across sample sizes. In the demand design, \verb|SmoothIV| performs best for sample size $n+m=1000$. Like hartford2017deep, we do not implement \verb|SmoothIV| for greater sample sizes due to running time. Among estimators that we are able to implement, \verb|KernelIV| performs best for sample sizes $n+m=5000$ and $n+m=10000$. \verb|KernelReg| ignores the instrument $Z$, and it is biased away from the structural function due to confounding noise $e$. This phenomenon can have counterintuitive consequences. Figure (ref) shows that in the highly nonlinear demand design, \verb|KernelReg| deviates further from the structural function as sample size increases because the algorithm is further misled by confounded data. Figure 2 of hartford2017deep documents the same effect when a feedforward neural network is used. The remaining algorithms make use of the instrument $Z$ to overcome this issue.
\verb|KernelIV| improves on \verb|SieveIV| in the same way that kernel ridge regression improves on ridge regression: by using an infinite dictionary of implicit basis functions rather than a finite dictionary of explicit basis functions. \verb|KernelIV| improves on \verb|SmoothIV| by using kernel ridge regression in not only stage 2 but also stage 1, avoiding costly density estimation. Finally, it improves on \verb|DeepIV| by directly learning stage 1 mean embeddings, rather than performing costly density estimation and sampling from the estimated density. The experiments show that \verb|KernelIV| works particularly well when the true structural function $h$ is smooth, confirming the theoretical guarantee of Theorem (ref). See Appendix (ref) for representative plots, implementation details, and a robustness study.
We introduce KIV, an algorithm for learning a nonlinear, causal relationship from confounded observational data. KIV is easily implemented and minimax optimal. As a contribution to the IV literature, we show how to estimate the stage 1 conditional expectation operator--an infinite by infinite dimensional object--by kernel ridge regression. As a contribution to the kernel methods literature, we show how the RKHS is well-suited to causal inference and ill-posed inverse problems. In simulations, KIV outperforms state of the art algorithms for nonparametric IV regression. The success of KIV suggests RKHS methods may be an effective bridge between econometrics and machine learning.
We are grateful to Alberto Abadie, Anish Agarwal, Michael Arbel, Victor Chernozhukov, Geoffrey Gordon, Jason Hartford, Motonobu Kanagawa, Anna Mikusheva, Whitney Newey, Nakul Singh, Bharath Sriperumbudur, and Suhas Vijaykumar. This project was made possible by the Marshall Aid Commemoration Commission.