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.
80,994 characters · 13 sections · 56 citation commands
Kernel Methods for Causal Functions: Dose, Heterogeneous, and Incremental Response Curves
Program evaluation aims to measure the counterfactual relationship between treatment $D$ and outcome $Y$, which may vary for different subpopulations: if we intervened on treatment, setting $D=d$, what would be the expected counterfactual outcome $Y^{(d)}$ for individuals with characteristics $V=v$? When treatment is binary, the causal parameter is a function $\theta_0(v)=E\{Y^{(1)}-Y^{(0)} \mid V=v\}$ called the heterogeneous treatment effect; when treatment is continuous, it is a function $\theta_0(d,v)=E\{Y^{(d)} \mid V=v\}$ that we call a heterogeneous response curve. Assuming selection on observable covariates $(V,X)$, the causal function $\theta_0(d,v)$ can be recovered by integrating the regression function $\gamma_0(d,v,x)=E(Y \mid D=d,V=v,X=x)$ according to the conditional distribution $\text{\normalfont pr}(x \mid v)$: $\theta_0(d,v)=\int \gamma_0(d,v,x)\mathrm{d}\text{\normalfont pr}(x \mid v)$ rosenbaum1983central,robins1986new, which may be complex when there are many covariates.
The same is true for other causal functions such as dose and incremental response curves, and even counterfactual distributions, albeit with different regressions and reweightings. Therefore nonparametric estimation of a causal function involves three challenging steps: estimating a nonlinear regression, with possibly many covariates; estimating the distribution for reweighting, which may be conditional; and using the nonparametric distribution to integrate the nonparametric regression. For this reason, flexible estimation of nonparametric causal functions, such as $\theta_0(d,v)$, is often deemed too computationally demanding to be practical for program evaluation.
Our key insight is that the reproducing kernel Hilbert space (RKHS), a popular nonparametric setting in machine learning, is precisely the class of functions for which the steps of nonparametric causal estimation can be separated. This decomposition follows almost immediately from the definition of the RKHS, and it is a specific strength of our framework; random forests, for example, do not allow such decoupling. Our key insight follows from a more fundamental one. Evaluation of a causal function is generally not a bounded functional over all of $\mathbb{L}^2$ van1991differentiable,newey1994asymptotic. We prove that the evaluation of a causal function is a bounded functional over the RKHS $\mathcal{H}$, which is a subset of $\mathbb{L}^2$. By the classic Riesz representation theorem of functional analysis, a bounded functional over a Hilbert space admits a decoupled inner product representation within the Hilbert space. We show how to use this representation to separate the steps of nonparametric causal estimation. This insight appears to be original.
We adapt kernel ridge regression, a classic machine learning algorithm that generalizes splines wahba1990spline, to address the computational challenges of estimating causal functions such as dose, heterogeneous, and incremental response curves. Nonparametric estimation with kernels is quite simple: the nonlinear regression with many covariates can be estimated by simple matrix operations; the conditional distribution can be expressed as a regression problem and estimated by simple matrix operations as well; and the step of integration can be performed by taking the product of the results. The final nonparametric estimator for the causal function has a one line, closed form solution, unlike previous work. This simplicity makes the family of estimators highly practical. The proposed estimators are substantially simpler yet outperform some leading alternatives in nonlinear simulations with many covariates; see Supplement (ref). As extensions, we generalize our new algorithmic techniques to counterfactual distributions in Supplement (ref) as well as causal functions and distributions identified by front and back door criteria in Supplement (ref).
Theoretically, our statistical guarantees rely on smoothness of the causal function and spectral decay of the covariance operator rather than the explicit dimension of treatment and covariates. In economic modelling, many variables may matter for labor market decisions, yet economic theory suggests that the effect of different intensities of job training should be well approximated by smooth functions. The emphasis on smoothness in the causal interpretation of RKHS assumptions generalizes standard Sobolev assumptions, and it differs from the emphasis on sparsity in lasso-type assumptions. Our causal function estimators are uniformly consistent with rates that combine minimax optimal rates for smooth nonparametric regressions. En route to our main results, we prove an improved rate for conditional expectation operators. Our main results are nonasymptotic and imply asymptotic uniform validity.
Conceptually, we illustrate how to use RKHS techniques in order to separate the steps of nonparametric causal estimation. In doing so, we provide a template for researchers to develop simple kernel estimators for complex causal estimands. Specifically, we clarify five assumptions under which we derive our various results: (i) identification, from the social scientific problem at hand; (ii) basic regularity conditions on the kernels, which are satisfied by all of the kernels typically used in practice; (iii) basic regularity on the outcome, treatment, and covariates, allowing them to be discrete or continuous variables that take values in general spaces (even texts, images, or graphs); (iv) smoothness of the causal estimand; and (v) spectral decay of the covariance operator. We combine these five assumptions to estimate causal functions, providing insight into the meaning and applicability of RKHS approximation assumptions for causal inference.
Statistically, we prove uniform consistency: our estimators converge to causal functions in $\sup$ norm, which encodes caution about worst case scenarios when informing policy decisions. Our finite sample rates of convergence explicitly account for each source of error at any finite sample size. Our rates do not directly depend on the data dimension, but rather the smoothness of the causal estimand and spectral decay of the covariance operator. The rates may indirectly depend on dimension; see Section (ref) for discussion in the context of the Sobolev space, which is a special case of an RKHS. Of independent interest, we provide a technical innovation to justify our main results: relative to previous work, we prove faster rates of convergence in Hilbert--Schmidt norm for conditional expectation operators. We generalize our main results to prove convergence in distribution for counterfactual distributions. The analysis of uniform confidence bands for our causal function estimators is an open question that we pose for future research, since uniform inference for kernel ridge regression remains an open question in statistics.
Computationally, we demonstrate state-of-the-art performance in nonlinear simulations with many covariates, despite the relative simplicity of our proposal compared to existing machine learning approaches. In order to simplify the causal estimation problem, we assume that underlying conditional expectation functions are elements in an RKHS. We propose a family of global estimators with closed form solutions, avoiding density estimation and sampling even for complex integrals. Throughout, the only hyperparameters are kernel hyperparameters and ridge regression penalties. The former have well established tuning procedures, and the latter are easily tuned using the closed form solution for generalized cross validation (which is asymptotically optimal) or leave-one-out cross validation (which we derive). In practice, the tunings are similar, and the asymptotically optimal choice aligns with our statistical theory.
Empirically, our kernel ridge regression approach allows for simple yet flexible estimation of nuanced causal estimands. Such estimands provide meaningful insights about the Job Corps, the largest job training program for disadvantaged youth in the US. Our key statistical assumption is that different intensities of job training have smooth effects on counterfactual employment, and those effects are smoothly modified by age. In our program evaluation in Supplement (ref), we find that the effect of job training on employment substantially varies by class hours and by age; a targeted policy will be more effective. Our program evaluation confirms earlier findings while also uncovering meaningful heterogeneity. We demonstrate how kernel methods for causal functions are a practical addition to the empirical economic toolkit.
We view nonparametric causal functions as reweightings of an underlying regression, synthesizing the $g$ formula robins1986new and partial means newey1994kernel frameworks. To express causal functions in this way, we build on canonical identification theorems under the assumption of selection on observables rosenbaum1983central,robins1986new,altonji2005cross. We propose simple, global estimators that combine kernel ridge regressions. Previous works that take a global view include van2003unified,luedtke2016super,diaz2013targeted,kennedy2020optimal, and references therein. A broad literature instead views causal functions as collections of localized treatment effects and proposes local estimators with Nadaraya--Watson smoothing, e.g. imai2004causal,rubin2005general,rubin2006extending,galvao2015uniformly,luedtke2016statistical,kennedy2017nonparametric,semenova2021debiased,kallus2018policy,chernozhukov2022debiased,fan2019estimation,zimmert2019nonparametric,colangelo2020double, and references therein. By taking a global view rather than a local view, we propose simple estimators that can be computed once and evaluated at any value of a continuous treatment, rather than a computationally intensive procedure that must be reimplemented at any treatment value.
Our work appears to be the first to reduce estimation of dose, heterogeneous, and incremental response curves to kernel ridge regressions. Previous works incorporating the RKHS into nonparametric estimation focus on different causal functions: nonparametric instrumental variable regression carrasco2007linear,darolles2011nonparametric,singh2019kernel, and heterogeneous treatment effect conditional on the full vector of covariates nie2021quasi. nie2021quasi propose the R learner to estimate the heterogeneous treatment effect $\theta_0(x)=E\{Y^{(1)}-Y^{(0)} \mid X=x\}$. nie2021quasi reviews the extensive literature that considers this estimand. The R learner minimizes a loss that contains inverse propensities and different regularization nie2021quasi, and it does not appear to have a closed form solution. The authors prove oracle mean square error rates. By contrast, we pursue a more general heterogeneous response curve with discrete or continuous treatment, conditional on some interpretable subvector $V$ abrevaya2015estimating: $\theta_0(d,v)=E\{Y^{(d)} \mid V=v\}$. Unlike previous work on nonparametric causal functions in the RKHS, we (i) consider dose, heterogeneous, and incremental response curves; (ii) propose estimators with closed form solutions; and (iii) prove uniform consistency, which is an important norm for policy evaluation.
We extend the framework from causal functions to counterfactual distributions. Existing work focuses on distributional generalizations of average treatment effect (ATE) or average treatment on the treated (ATT) for binary treatment firpo2007efficient,cattaneo2010efficient,chernozhukov2013inference, e.g. $\theta_0=\text{\normalfont pr}\{Y^{(1)}\}-\text{\normalfont pr}\{Y^{(0)}\}$. muandet2021counterfactual propose an RKHS approach for distributional ATE and ATT with binary treatment using inverse propensity scores and an assumption on the smoothness of a ratio of densities, which differs from our approach. Unlike previous work, we (i) allow treatment to be continuous; (ii) avoid inversion of propensity scores and densities; and (iii) study a broad class of counterfactual distributions for the full population, subpopulations, and alternative populations, e.g. $\theta_0(d,v)=\text{\normalfont pr}\{Y^{(d)} \mid V=v\}$.
We provide a detailed comparison with kernel methods for binary treatment effects in Section (ref). Whereas we study causal functions, these works study causal scalars kallus2020generalized,hirshberg2019minimax,singh2021debiased. We clarify the sense in which our causal function estimators generalize known estimators for treatment effects to new estimators for causal functions. Previous work is inherently tied to the $\mathbb{L}^2$ bounded functional perspective. However, evaluation of a causal function is not a bounded functional over all of $\mathbb{L}^2$ van1991differentiable,newey1994asymptotic. Therefore our algorithms extend the conceptual framework of kernel methods for causal inference in a new direction. Our statistical contribution is a new, uniform analysis of response curves that goes beyond pointwise approximation of response curves by local treatment effects.
This paper subsumes our previous draft singh2020kernel.
A causal function summarizes the expected counterfactual outcome $Y^{(d)}$ given a hypothetical intervention on continuous treatment that sets $D=d$. The causal inference literature studies a rich variety of causal functions with nuanced interpretation, which we define below. Unless otherwise noted, expectations are with respect to the population distribution $\text{\normalfont pr}$.
The superscript of each nonparametric causal function corresponds to its familiar parametric analogue. Results for means of potential outcomes immediately imply results for differences thereof. See Supplement (ref) for counterfactual distributions and Supplement (ref) for graphical models.
The dose response curves $\theta_0^{ATE}(d)$ and $\theta_0^{DS}(d,\tilde{\text{\normalfont pr}})$ are causal functions for entire populations. The second argument of $\theta_0^{DS}(d,\tilde{\text{\normalfont pr}})$ gets to the heart of external validity: though our data were drawn from population $\text{\normalfont pr}$, what would be the dose response curve for a different population $\tilde{\text{\normalfont pr}}$? For example, a job training study may be conducted in Virginia, yet we may wish to inform policy in Arkansas, a state with different demographics hotz2005predicting. Predictive questions of this nature are widely studied in machine learning under the names of transfer learning, distribution shift, and covariate shift quinonero2009dataset.
$\theta_0^{ATE}(d)$ and $\theta_0^{DS}(d,\tilde{\text{\normalfont pr}})$ are dose response curves for entire populations, but causal functions may vary for different subpopulations. Towards the goal of personalized or targeted interventions, an analyst may ask another nuanced counterfactual question: what would have been the effect of treatment $D=d'$ for the subpopulation who actually received treatment $D=d$? When treatment is continuous, we may define the conditional response $ \theta_0^{ATT}(d,d')=E\{Y^{(d')} \mid D=d\}$.
In $\theta_0^{ATT}(d,d')$, heterogeneity is indexed by treatment $D$. Heterogeneity may instead be indexed by some interpretable covariate subvector $V$, e.g. age, race, or gender abrevaya2015estimating. An analyst may therefore prefer to measure heterogeneous effects for subpopulations characterized by different values of $V$. For simplicity, we will write covariates as $(V,X)$ for this setting, where $X$ are additional identifying covariates besides the interpretable covariates $V$. While many works focus on the special case where treatment is binary, our definition of heterogeneous response curve $\theta_0^{CATE}(d,v)=E\{Y^{(d)} \mid V=v\}$ allows for continuous treatment.
Lemma (ref) clarifies the data requirements for estimating each causal function. The dose response $\theta_0^{ATE}(d)$ and conditional response $\theta_0^{ATT}(d,d')$ require observations of outcome $Y$, treatment $D$, and covariates $X$ drawn from the population $\text{\normalfont pr}$. The dose response with distribution shift $\theta_0^{DS}(d,\tilde{\text{\normalfont pr}})$ additionally requires observations of covariates $\tilde{X}$ drawn from the alternative population $\tilde{\text{\normalfont pr}}$. For the heterogeneous response $\theta_0^{CATE}(d,v)$, we abuse notation by denoting the covariates by $(V,X)$, where $V$ is the subcovariate of interest and selection is with respect to the union $(V,X)$. An analyst requires observations of $(Y,D,V,X)$ drawn from the population $\text{\normalfont pr}$.
In particular, Lemma (ref) expresses each causal function as an integral of the regression function $\gamma_0$ according to a marginal or conditional distribution. As previewed in Section (ref), nonparametric estimation of $\theta_0^{CATE}(d,v)$ involves three steps: estimating a nonlinear regression $\gamma_0(d,v,x)$, which may involve many covariates $X$; estimating the conditional distribution $\text{\normalfont pr}(x \mid v)$ for reweighting; and using the latter to integrate the former. In the next section, we propose original estimators that achieve all three steps in a one line, closed form solution.
To present the algorithm, we provide background on the RKHS. The essential property of a function $\gamma$ in an RKHS $\mathcal{H}$ is the eponymous reproducing property: $\gamma(w)=\langle \gamma,\phi(w)\rangle_{\mathcal{H}}$ where $\phi(w)$ are features, formally defined below, that serve as the basis functions for $\mathcal{H}$. Our key algorithmic insight is to interpret the reproducing property as a way to separate the function $\gamma$ from the features $\phi(w)$. We use this defining property of the RKHS to decouple the three steps of nonparametric causal estimation. After providing RKHS background material, we prove an inner product representation that formalizes the decoupling, then introduce the causal estimators.
A scalar-valued RKHS $\mathcal{H}$ is a Hilbert space with elements that are functions $\gamma:\mathcal{W}\rightarrow\mathbb{R}$, on which the operator of evaluation is bounded berlinet2011reproducing. Polynomial, spline, and Sobolev spaces are widely used examples of RKHSs. $\mathcal{W}$ can be any Polish space, so a value $w\in \mathcal{W}$ can be discrete or continuous. An RKHS is fully characterized by its feature map, which takes a point $w$ in the original space $\mathcal{W}$ and maps it to a feature $\phi(w)$ in the RKHS $\mathcal{H}$. The closure of $span\{\phi(w)\}_{w\in \mathcal{W}}$ is the RKHS $\mathcal{H}$. In other words, $\{\phi(w)\}_{w\in\mathcal{W}}$ can be viewed as the dictionary of basis functions for the RKHS $\mathcal{H}$. The kernel $k:\mathcal{W}\times \mathcal{W}\rightarrow \mathbb{R}$ is the inner product of features $\phi(w)$ and $\phi(w')$: $ k(w,w')=\langle \phi(w),\phi(w') \rangle_{\mathcal{H}} $. A real-valued kernel $k$ is continuous, symmetric, and positive definite. Though we have constructed the kernel from the feature map, the Moore--Aronszajn Theorem states that, for any positive definite kernel $k$, there exists a unique RKHS $\mathcal{H}$ with feature map $\phi:w\mapsto k(w,\cdot)$. We have already seen that if $\gamma\in\mathcal{H}$, then $\gamma:\mathcal{W}\rightarrow \mathbb{R}$. With the additional notation of the feature map, we write $ \gamma(w)= \langle \gamma,\phi(w) \rangle_{\mathcal{H}} $. If $\mathcal{W}$ is separable and $\phi$ is continuous, then $\mathcal{H}$ is separable and may be infinite dimensional.
The RKHS is a practical hypothesis space for nonparametric regression. Consider output $Y\in\mathbb{R}$, input $W\in\mathcal{W}$, and the goal of estimating the conditional expectation function $\gamma_0(w)=E(Y \mid W=w)$. A kernel ridge regression estimator of $\gamma_0$ is
$\lambda>0$ is a hyperparameter on the ridge penalty $\|\gamma\|^2_{\mathcal{H}}$, which imposes smoothness in estimation. The solution to the optimization problem has a well known closed form kimeldorf1971some, which we exploit and generalize throughout this work:
The closed form solution involves the kernel matrix $K_{WW}\in\mathbb{R}^{n\times n}$ with $(i,j)$th entry $k(W_i,W_j)$, and the kernel vector $K_{Ww}\in\mathbb{R}^n$ with $i$th entry $k(W_i,w)$. To tune the ridge hyperparameter $\lambda$, both generalized cross validation and leave-one-out cross validation have closed form solutions, and the former is asymptotically optimal craven1978smoothing,li1986asymptotic.
We have seen that the feature map takes a value in the original space $w\in\mathcal{W}$ and maps it to a feature in the RKHS $\phi(w)\in\mathcal{H}$. Now we generalize this idea, from the embedding of a value $w$ to the embedding of a distribution $\text{\normalfont q}$. Just as a value $w$ in the original space is embedded as an element $\phi(w)$ in the RKHS, so too the distribution $\text{\normalfont q}$ over the original space can be embedded as an element $\mu=E_{\text{\normalfont q}}\{\phi(W)\}$ in the RKHS smola2007hilbert,berlinet2011reproducing. Boundedness of the kernel implies existence of the mean embedding as well as Bochner integrability, which permits us to exchange the expectation and inner product. Mean embeddings facilitate the evaluation of expectations of RKHS functions: for $\gamma \in \mathcal{H}$, $E_{\text{\normalfont q}}\{\gamma(W)\}=E_{\text{\normalfont q}}\left\{\langle \gamma,\phi(W) \rangle_{\mathcal{H}}\right\}=\langle \gamma,\mu \rangle_{\mathcal{H}}$. The final expression foreshadows how we will use the technique of mean embeddings to decouple the nonparametric regression step from the nonparametric reweighting step in the estimation of causal functions. A natural question is whether the embedding $\text{\normalfont q}\mapsto E_{\text{\normalfont q}}\{\phi(W)\}$ is injective, i.e. whether the RKHS element representation is unique. This is called the characteristic property of the kernel $k$, and it holds for commonly used RKHSs e.g. the exponentiated quadratic kernel sriperumbudur2010relation.
The tensor product RKHS is one way to construct an RKHS for functions with multiple arguments. Consider the RKHSs $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ with positive definite kernels $k_{1}:\mathcal{W}_1\times\mathcal{W}_1\rightarrow \mathbb{R}$ and $k_{2}:\mathcal{W}_2\times\mathcal{W}_2\rightarrow \mathbb{R}$, respectively. An element $\gamma_1\in \mathcal{H}_{1}$ is a function $\gamma_1:\mathcal{W}_1\rightarrow \mathbb{R}$ and an element $\gamma_2\in \mathcal{H}_{2}$ is a function $\gamma_2:\mathcal{W}_2\rightarrow \mathbb{R}$. The tensor product RKHS $\mathcal{H}=\mathcal{H}_{1}\otimes \mathcal{H}_{2}$ is the RKHS with the product kernel $ k:(\mathcal{W}_1\times \mathcal{W}_2) \times (\mathcal{W}_1\times \mathcal{W}_2)\rightarrow \mathbb{R},\; \{(w_1,w_2),(w'_1,w'_2)\}\mapsto k_{1}(w_1,w_1') k_{2}(w_2,w_2') $. Equivalently, the tensor product RKHS $\mathcal{H}$ has feature map $\phi(w_1)\otimes \phi(w_2)$ such that $\|\phi(w_1)\otimes \phi(w_2)\|_{\mathcal{H}}=\|\phi(w_1)\|_{\mathcal{H}_1}\|\phi(w_2)\|_{\mathcal{H}_2}$. Formally, tensor product notation means $(a\otimes b)c=a \langle b,c\rangle$. An element of the tensor product RKHS $\gamma \in\mathcal{H}$ is a function $\gamma:\mathcal{W}_1\times \mathcal{W}_2 \rightarrow \mathbb{R}$. We assume that the regression function $\gamma_0(w_1,w_2)=E(Y \mid W=w_1,w_2=w_2)$ is an element of a tensor product RKHS, i.e. $\gamma_0\in \mathcal{H}$. As such, the different arguments of $\gamma_0$ are decoupled, which we exploit when calculating partial means.
Finally, we introduce the RKHS $\mathcal{L}_2(\mathcal{H}_{1},\mathcal{H}_{2})$ that we employ for conditional expectation operators. Rather than being a space of real-valued functions, it is a space of Hilbert--Schmidt operators from one RKHS to another. If the operator $E$ is an element of $ \mathcal{L}_2(\mathcal{H}_{1},\mathcal{H}_{2})$, then $E:\mathcal{H}_{1} \rightarrow \mathcal{H}_{2}$. Formally, it can be shown that $\mathcal{L}_2(\mathcal{H}_{1},\mathcal{H}_{2})$ is an RKHS in its own right with an appropriately defined kernel and feature map. $\mathcal{L}_2(\mathcal{H}_{1},\mathcal{H}_{2})$ is an example of a vector-valued RKHS; see micchelli2005learning for a more general discussion. In the present work, we assume the conditional expectation operator $E_0:\gamma_1(\cdot)\mapsto E\{\gamma_1(W_1) \mid W_2=\cdot\}$ is an element of this RKHS, i.e. $E_0\in \mathcal{L}_2(\mathcal{H}_{1},\mathcal{H}_{2})$. We estimate $E_0$ by a kernel ridge regression in $\mathcal{L}_2(\mathcal{H}_{1},\mathcal{H}_{2})$, which coincides with estimating the conditional mean embedding $\mu_{w_1}(w_2)=E\{\phi(W_1) \mid W_2=w_2\}$ via the kernel ridge regression of $\phi(W_1)$ on $\phi(W_2)$; see the derivation of Algorithm (ref) below.
Lemma (ref) makes precise how each causal function is identified as a partial mean of the form $\int \gamma_0(d,x)\mathrm{d}\text{\normalfont q}$ for some distribution $\text{\normalfont q}$. To facilitate estimation, we now assume that $\gamma_0$ is an element of an RKHS. In our construction, we define scalar valued RKHSs for treatment $D$ and covariates $(V,X)$, then assume that the regression is an element of the tensor product space. Let $k_{\mathcal{D}}:\mathcal{D}\times \mathcal{D}\rightarrow \mathbb{R}$, $k_{\mathcal{V}}:\mathcal{V}\times \mathcal{V}\rightarrow \mathbb{R}$, and $k_{\mathcal{X}}:\mathcal{X}\times \mathcal{X}\rightarrow \mathbb{R}$ be measurable positive definite kernels corresponding to scalar valued RKHSs $\mathcal{H}_{\mathcal{D}}$, $\mathcal{H}_{\mathcal{V}}$, and $\mathcal{H}_{\mathcal{X}}$. Denote the feature maps $ \phi_{\mathcal{D}}:\mathcal{D}\rightarrow \mathcal{H}_{\mathcal{D}}, \; d\mapsto k_{\mathcal{D}}(d,\cdot ); \; \phi_{\mathcal{V}}:\mathcal{V}\rightarrow \mathcal{H}_{\mathcal{V}}, \; v\mapsto k_{\mathcal{V}}(v,\cdot ); \; \phi_{\mathcal{X}}:\mathcal{X}\rightarrow \mathcal{H}_{\mathcal{X}}, \; x\mapsto k_{\mathcal{X}}(x,\cdot ) $. To lighten notation, we suppress subscripts when arguments are provided.
For $\theta_0^{ATE}$, $\theta_0^{DS}$, and $\theta_0^{ATT}$, we assume the regression $\gamma_0$ is an element of the RKHS $\mathcal{H}$ with the kernel $k(d,x;d',x')=k_{\mathcal{D}}(d,d')k_{\mathcal{X}}(x,x')$. We appeal to the fact that the product of positive definite kernels for $\mathcal{H}_{\mathcal{D}}$ and $\mathcal{H}_{\mathcal{X}}$ defines a new positive definite kernel for $\mathcal{H}$. The product construction provides a rich composite basis; $\mathcal{H}$ has the tensor product feature map $\phi(d)\otimes \phi(x)$ and $\mathcal{H}=\mathcal{H}_{\mathcal{D}}\otimes \mathcal{H}_{\mathcal{X}}$. In this RKHS, $ \gamma_0(d,x)=\langle \gamma_0, \phi(d)\otimes \phi(x)\rangle_{\mathcal{H}} $. Likewise for $\theta_0^{CATE}$ we assume $\gamma_0\in \mathcal{H}=\mathcal{H}_{\mathcal{D}}\otimes\mathcal{H}_{\mathcal{V}}\otimes \mathcal{H}_{\mathcal{X}}$. We place regularity conditions on this RKHS construction in order to represent causal functions as inner products in $\mathcal{H}$. In anticipation of counterfactual distributions in Supplement (ref), we also include conditions for an outcome RKHS in parentheses.
Commonly used kernels are continuous and bounded. Measurability is a similarly weak condition. The characteristic property ensures injectivity of the mean embeddings.
See Supplement (ref) for the full proof. $\mu_x(v)=\int\phi(x) \text{\normalfont pr}(x \mid v)$ is the mean embedding of the conditional distribution $\text{\normalfont pr}(x \mid v)$. It encodes the distribution $\text{\normalfont pr}(x \mid v)$ as a function $\mu_x(v)\in\mathcal{H}_{\mathcal{X}}$ such that the causal function $\theta_0^{CATE}(d,v)$ can be expressed as an inner product in $\mathcal{H}$.
The representation in Theorem (ref) is essential to the algorithm derivation. In particular, the representation cleanly separates the three steps necessary to estimate a causal function: estimating a nonlinear regression, which may involve many covariates; estimating the distribution for reweighting; and using the nonparametric distribution to integrate the nonparametric regression. For example, for $\theta_0^{CATE}(d,v)$, our estimator is $\hat{\theta}^{CATE}(d,v)=\langle \hat{\gamma}, \phi(d)\otimes \phi(v)\otimes \hat{\mu}_x(v)\rangle_{\mathcal{H}}$. The nonlinear regression estimator $\hat{\gamma}$ is a standard kernel ridge regression of $Y$ on $\phi(D)\otimes \phi(V)\otimes \phi(X)$; the reweighting distribution estimator $\hat{\mu}_x(v)$ is a generalized kernel ridge regression of $\phi(X)$ on $\phi(V)$; and the latter can be used to integrate the former by simply multiplying the two. This algorithmic insight is a key innovation of the present work, and the reason why our estimators have simple closed form solutions despite complicated causal integrals.
See Supplement (ref) for the full derivation and a comparison to series estimation. We give theoretical values for $(\lambda,\lambda_1,\lambda_2)$ that optimally balance bias and variance in Theorem (ref) below. Supplement (ref) gives practical tuning procedures based on generalized and leave-one-out cross validation to empirically balance bias and variance, the former of which is asymptotically optimal.
We now connect our kernel methods for causal functions with related kernel methods for treatment effects. Recall the definition $\theta_0^{ATE}(d)=E\{Y^{(d)}\}$. We allow treatment to be continuous, so $\theta_0^{ATE}$ is a causal function called the dose response. In related work, treatment is binary, so $\theta_0^{ATE}$ is a vector of two causal scalars $\theta^{ATE}_0(1),\theta^{ATE}_0(0)$ whose difference is the treatment effect.
We clarify three points. (i) There is a sense in which our algorithms generalize known estimators for treatment effects to new estimators for causal functions. (ii) A treatment effect is a bounded functional over $\mathbb{L}^2$ with a balancing weight representation, while a response curve is not. Our key insight is that a response curve is a bounded functional over the RKHS $\mathcal{H}$, which is a subset of $\mathbb{L}^2$. (iii) Our theoretical contribution is a new, uniform analysis of response curves. The analysis goes beyond pointwise approximation of response curves by local treatment effects.
We begin by reviewing the theory of balancing weights, which are popular in causal inference with binary treatments. For clarity, in this section we emphasize a fixed treatment value by writing $d^*\in \mathcal{D}$. The following representation is well known.
In summary, a treatment effect has two representations: the primal representation of Lemma (ref) as a partial mean of the regression $\gamma_0(d,x)=E(Y \mid D=d,X=x)$, and the dual representation of Proposition (ref) as a reweighting of the outcome $Y$ using the balancing weight $\alpha_0(d,x)=1(d=d^*)/\text{\normalfont pr}(D=d^* \mid x)$. Clearly, the two representations are related by the law of iterated expectations. Moreover, from the closed form of $\alpha_0$, we require $\text{\normalfont pr}(D=d^* \mid X)>0$ for $\alpha_0$ to exist. This property keenly relies on the treatment being discrete. Indeed, it is well known that a balancing weight representation does not exist for response curves.
Whereas a binary treatment effect is a bounded functional over $\mathbb{L}^2$ with a balancing weight representation, a dose response is not a bounded functional over $\mathbb{L}^2$ and does not have a balancing weight representation in the classic sense. From a functional analytic perspective, this discrepancy is the reason why the problems we study are nonparametric whereas previous work on kernel methods for treatment effects are semiparametric. See Supplement (ref) for discussion.
Our key insight is that the dose response is a bounded functional over the RKHS $\mathcal{H}$, which is a subset of $\mathbb{L}^2$. This fact follows from three simple observations: (i) the dose response is a partial mean; (ii) in the RKHS, a partial mean can be reformulated as a kind of evaluation; and (iii) the RKHS $\mathcal{H}$ is the subset of $\mathbb{L}^2$ for which evaluation is a bounded functional. Through this lens, Theorem (ref) shows that there can exist a function $\tilde{\alpha}_0\in\mathcal{H}$ such that $\theta_0^{ATE}(d)=\langle \gamma_0,\tilde{\alpha}_0 \rangle_{\mathcal{H}}$ even when there does not exist a function $\alpha_0\in \mathbb{L}^2$ such that $\theta_0^{ATE}(d)=\langle \gamma_0,\alpha_0 \rangle_{\mathbb{L}^2}$.
What is the relationship between between our kernel methods for causal functions and existing kernel methods for treatment effects? There is a sense in which our dose response estimator, which is the simplest case of our framework, is a relaxation of kernel balancing weight estimators from binary treatment to continuous treatment. We formalize this connection as follows.
See Supplement (ref) for the proof. The balancing weight estimator $\hat{\alpha}$ minimizes a generalized balancing weight loss with ridge regularization; see kallus2020generalized, hirshberg2019minimax, and singh2021debiased for various formulations. Corollary (ref) provides intuition for our tensor product RKHS construction. Our product kernel construction ensures that using the binary treatment kernel amounts to subsetting, which recovers previous algorithms. The tensor product RKHS provides a natural way to relax binary treatment to continuous treatment while retaining computational tractability.
As argued in Proposition (ref), the balancing weight $\alpha_0\in \mathbb{L}^2$ does not exist for the dose response. Nonetheless, our key insight in Theorem (ref) is that a function $\tilde{\alpha}_0\in\mathcal{H}$ does exist to serve a similar purpose. By combining the partial mean perspective with the technique of kernel mean embedding, we demonstrate that our framework easily extends to conditional nonparametric causal functions, e.g. the heterogeneous response curve $\theta_0^{CATE}(d,v)$, which are substantially more challenging than unconditional nonparametric causal functions, e.g. the dose response $\theta^{ATE}_0(d)$.
Perhaps the most surprising consequence of our construction is the closed form solution for causal functions. In particular, each closed form solution is a reweighting of the observed outcomes with empirical weights that we characterize even though a population balancing weight in $\mathbb{L}^2$ does not exist. In sum, previous work kallus2020generalized,hirshberg2019minimax,singh2021debiased on kernel methods for treatment effects is inherently tied to the $\mathbb{L}^2$ population balancing weight perspective; our algorithms apply the conceptual framework of kernel methods to new classes of causal functions. The following corollary reinterprets Algorithm (ref) through this lens.
Each of our proposed causal function estimators is global. In particular, within Corollary (ref), the weights $(\hat{\alpha}_j^{ATE},\hat{\alpha}_j^{DS},\hat{\alpha}_j^{ATT},\hat{\alpha}_j^{CATE})$ $(j=1,...,n)$ depend on all of the observations as refracted through the ridge regularized empirical covariance and the kernel evaluations $k(D_i,d)$. This approach departs from a localization approach to causal functions whereby the weight assigned to each observation is determined by Nadaraya--Watson smoothing kennedy2017nonparametric,kallus2018policy,colangelo2020double,chernozhukov2022debiased. In the localization approach, the weight is $k^{NW}\{(D_i-d)/h\}$ where $k^{NW}$ is a Nadaraya--Watson kernel and $h$ is a vanishing bandwidth. By contrast, we consider a fixed kernel and vanishing ridge regularization.
The global perspective has three main advantages. First, our estimators can be computed once and evaluated at any value of a continuous treatment. By contrast, a localized estimator is a computationally intensive procedure that must be reimplemented at any treatment value. Second, our estimators are constructed from function classes with designed-in smoothness properties, which leads to smoother and therefore more plausible response curves. We compare our smooth estimate with a jagged localizing estimate in the program evaluation of Supplement (ref). Third, we prove uniform consistency of response curves, whereas localizations of previous results would only lead to pointwise consistency. These uniform guarantees are the focus of the next section.
In Section (ref), we defined the causal functions of interest, and identified them as partial means. In Section (ref), we introduced the tensor product RKHS as the function space in which the three steps of nonparametric causal estimation may be decoupled. We then proposed estimators based on kernel ridge regression with closed form solutions. In Section (ref), we demonstrated that our estimators generalize known estimators for the binary treatment case. In this section, we prove uniform consistency of the estimators, with finite sample rates that combine minimax optimal rates. To do so, we define our key approximation assumptions, which are standard in RKHS learning theory: smoothness and spectral decay.
To state our key assumptions, we must introduce a certain eigendecomposition. Recall the example of a generic RKHS $\mathcal{H}$ with kernel $k:\mathcal{W}\times \mathcal{W}\rightarrow \mathbb{R}$ consisting of functions $\gamma:\mathcal{W}\rightarrow \mathbb{R}$. Let $\nu$ be any Borel measure on $\mathcal{W}$. We denote by $\mathbb{L}^2_{\nu}(\mathcal{W})$ the space of square integrable functions with respect to measure $\nu$. Given the kernel, define the integral operator $ L:\mathbb{L}_{\nu}^2(\mathcal{W})\rightarrow \mathbb{L}_{\nu}^2(\mathcal{W}),\; \gamma \mapsto \int k(\cdot,w)\gamma(w)\mathrm{d}\nu(w) $. If the kernel $k$ is defined on $\mathcal{W}\subset\mathbb{R}^d$ and shift invariant, then $L$ is a convolution of $k$ and $\gamma$. If $k$ is smooth, then $L\gamma$ is a smoothed version of $\gamma$. $L$ is a self adjoint, positive, compact operator, so by the spectral theorem we can denote its countable eigenvalues by $(\eta_j)$ and its countable eigenfunctions, which are equivalence classes, by $\{(\varphi_j)_{\nu}\}$: $$ L\gamma =\sum_{j=1}^{\infty} \eta_j\langle (\varphi_j)_{\nu},\gamma \rangle_{\mathbb{L}^2_{\nu}(\mathcal{W})} (\varphi_j)_{\nu},\quad (\varphi_j)_{\nu}=\{f:\nu(f\neq \varphi_j)=0\}. $$ Without loss of generality, $\eta_j\geq \eta_{j+1}$, and these are also the eigenvalues of the feature covariance operator $T=E\{\phi(W)\otimes \phi(W)\}$. For simplicity, we assume $(\eta_j)>0$ in this discussion; see cucker2002mathematical for the more general case. $\{(\varphi_j)_{\nu}\}$ form an orthonormal basis of $\mathbb{L}_{\nu}^2(\mathcal{W})$. By the generalized Mercer's Theorem for Polish spaces steinwart2012mercer, we can express the kernel as $k(w,w')=\sum_{j=1}^{\infty}\eta_j \varphi_j(w)\varphi_j(w')$, where $(w,w')$ are in the support of $\nu$, $\varphi_j$ is a continuous element in the equivalence class $(\varphi_j)_{\nu}$, and the convergence is absolute and uniform.
With this notation, we express $\mathbb{L}_{\nu}^2(\mathcal{W})$ and the RKHS $\mathcal{H}$ in terms of the series $\{(\varphi_j)_{\nu}\}$. If $\gamma\in \mathbb{L}_{\nu}^2(\mathcal{W})$, then $\gamma$ can be uniquely expressed as $ \gamma=\sum_{j=1}^{\infty}\gamma_j(\varphi_j)_{\nu} $ and the partial sums $\sum_{j=1}^J \gamma_j (\varphi_j)_{\nu}$ converge to $\gamma$ in $\mathbb{L}^2_{\nu}(\mathcal{W})$. Indeed, for $\gamma=\sum_{j=1}^{\infty}\gamma_j(\varphi_j)_{\nu}$ and $\gamma'=\sum_{j=1}^{\infty}\gamma_j'(\varphi_j)_{\nu}$, $$ \mathbb{L}^2_{\nu}(\mathcal{W})=\left\{\gamma=\sum_{j=1}^{\infty}\gamma_j(\varphi_j)_{\nu}:\; \sum_{j=1}^{\infty}\gamma_j^2<\infty\right\},\quad \langle \gamma,\gamma' \rangle_{\mathbb{L}^2_{\nu}(\mathcal{W})}=\sum_{j=1}^{\infty} \gamma_j\gamma_j'. $$ By cucker2002mathematical, the RKHS $\mathcal{H}$ can be explicitly represented as $$ \mathcal{H}=\left(\gamma=\sum_{j=1}^{\infty}\gamma_j\varphi_j:\;\sum_{j=1}^{\infty} \frac{\gamma_j^2}{\eta_j}<\infty\right),\quad \langle \gamma,\gamma' \rangle_{\mathcal{H}}=\sum_{j=1}^{\infty} \frac{\gamma_j\gamma_j'}{\eta_j}. $$ To interpret this result, recall that $(\eta_j)$ is a weakly decreasing sequence. The RKHS $\mathcal{H}$ is the subset of functions in $\mathbb{L}^2_{\nu}(\mathcal{W})$ which are continuous and for which higher order terms in the series $\{(\varphi_j)_{\nu}\}$ have a smaller contribution. The RKHS inner product penalizes higher order coefficients, and the magnitude of the penalty corresponds to how small the eigenvalue is.
We have seen how to conduct kernel ridge regression with the RKHS $\mathcal{H}$. To analyze the bias from ridge regularization, we place a smoothness assumption called the source condition on the regression function $\gamma_0(w)=E(Y \mid W=w)$ smale2007learning,caponnetto2007optimal,carrasco2007linear. Formally, we place assumptions of the form
While $c=1$ recovers correct specification $\gamma_0\in \mathcal{H}$, $c\in(1,2]$ is a stronger condition: $\gamma_0$ is a particularly smooth element of $\mathcal{H}$, well approximated by the leading terms in the series $\{(\varphi_j)_{\nu}\}$. Smoothness delivers uniform consistency. A larger value of $c$ corresponds to a smoother target $\gamma_0$ and a faster convergence rate for $\hat{\gamma}$. Rates do not further improve for $c>2$, which is known as the saturation effect for ridge regularization.
To analyze the variance of kernel ridge regression, we place a spectral decay assumption called the effective dimension of the basis $(\varphi_j)$ for the RKHS $\mathcal{H}$. To obtain faster convergence rates, we place a direct assumption on the rate at which the eigenvalues $(\eta_j)$, and hence the importance of the eigenfunctions $(\varphi_j)$, decay: we assume there exists some constant $C$ such that for all $j$
A bounded kernel, which we have already assumed, implies $b=1$ fischer2017sobolev. The limit $b\rightarrow \infty$ may be interpreted as a finite dimensional RKHS caponnetto2007optimal. For intermediate values of $b$, the polynomial rate of spectral decay quantifies the effective dimension of the RKHS $\mathcal{H}$ in light of the measure $\nu$. Intuitively, a higher value of $b$ corresponds to a lower effective dimension and a faster convergence rate for $\hat{\gamma}$.
For intuition, we relate the source condition and effective dimension to a familiar notion of smoothness in the Sobolev space. The restriction that defines an RKHS generalizes higher order smoothness in a Sobolev space. Indeed, certain Sobolev spaces are RKHSs. Let $\mathcal{W}\subset \mathbb{R}^p$. Denote by $\mathbb{H}_2^s$ the Sobolev space with $s>p/2$ derivatives that are square integrable. This space can be generated by the Mat\`ern kernel, which converges to the popular exponentiated quadratic kernel as $s\rightarrow \infty$. Suppose $\mathcal{H}=\mathbb{H}_2^s$ is chosen as the RKHS for estimation. Suppose the measure $\nu$ supported on $\mathcal{W}$ is absolutely continuous with respect to the uniform distribution and bounded away from zero. If $\gamma_0\in \mathbb{H}_2^{s_0}$, then $c=s_0/s$ pillaud2018statistical. Written another way, $(\mathbb{H}_2^{s})^c=\mathbb{H}_2^{s_0}$. In this sense, $c$ precisely quantifies the additional smoothness of $\gamma_0$ relative to $\mathcal{H}$. Moreover, in this Sobolev space, $b=2s/p>1$ fischer2017sobolev. The effective dimension is increasing in the input dimension $p$ and decreasing in the degree of smoothness $s$. The minimax optimal rate in Sobolev norm is $n^{-(c-1)/\{2(c+1/b)\}}=n^{-(s_0-s)/(2s_0+p)}$, which is achieved by kernel ridge regression with the rate optimal regularization $\lambda=n^{-1/(c+1/b)}=n^{-2s/(2s_0+p)}$. Our analysis applies to Sobolev spaces over $\mathbb{R}^p$ as a special case; our results are much more general, allowing treatment and covariates to be in Polish spaces.
Towards a guarantee of uniform consistency, we place regularity conditions on the original spaces. In anticipation of counterfactual distributions in Supplement (ref), we also include conditions for the outcome space in parentheses.
A Polish space is a separable and completely metrizable topological space. Random variables with support in a Polish space may be discrete or continuous and may even be infinite dimensional. Bounded $Y$ implies the moment condition.
Next, we assume the regression $\gamma_0$ is smooth in the sense of (ref), and $\mathcal{H}$ has low effective dimension in the sense of (ref). Denote the $j$th eigenvalue of the convolution operator for $\mathcal{H}$ by $\eta_j(\mathcal{H})$. Recall that $\eta_j(\mathcal{H})$ is also the $j$th eigenvalue of the feature covariance operator.
See Supplement (ref) for alternative ways of writing and interpreting Assumption (ref). We place similar smoothness and spectral decay conditions on the conditional mean embeddings $\mu_x(d)$ and $\mu_x(v)$, which are generalized conditional expectation functions. We articulate this assumption abstractly for the conditional mean embedding $\mu_{a}(b)=\int \phi(a)\mathrm{d}\text{\normalfont pr}(a \mid b)$ where $a\in\mathcal{A}_{\ell}$ and $b\in\mathcal{B}_{\ell}$. All one has to do is specify $\mathcal{A}_{\ell}$ and $\mathcal{B}_{\ell}$ to specialize the assumption. For $\mu_x(d)$, $\mathcal{A}_1=\mathcal{X}$ and $\mathcal{B}_1=\mathcal{D}$; for $\mu_x(v)$, $\mathcal{A}_2=\mathcal{X}$ and $\mathcal{B}_2=\mathcal{V}$. For fixed $\mathcal{A}_{\ell}$ and $\mathcal{B}_{\ell}$, we parametrize smoothness by $c_{\ell}$ and spectral decay by $b_{\ell}$.
Formally, define the conditional expectation operator $E_{\ell}:\mathcal{H}_{\mathcal{A}_{\ell}}\rightarrow\mathcal{H}_{\mathcal{B}_{\ell}}$, $f(\cdot)\mapsto E\{f(A_{\ell}) \mid B_{\ell}=\cdot\}$. By construction, $E_{\ell}$ encodes the same information as $\mu_{a}(b)$ since $$ \{\mu_{a}(b)\}(\cdot)=\int \phi(a)\mathrm{d}\text{\normalfont pr}(a \mid b) =\{E_{\ell}\phi(\cdot)\}(b) =\{E_{\ell}^* \phi(b)\}(\cdot),\quad a\in\mathcal{A}_{\ell}, \quad b\in\mathcal{B}_{\ell}, $$ where $E_{\ell}^*$ is the adjoint of $E_{\ell}$. We denote the space of Hilbert--Schmidt operators between $\mathcal{H}_{\mathcal{A}_{\ell}}$ and $\mathcal{H}_{\mathcal{B}_{\ell}}$ by $\mathcal{L}_2(\mathcal{H}_{\mathcal{A}_{\ell}},\mathcal{H}_{\mathcal{B}_{\ell}})$. grunewalder2013smooth,singh2019kernel prove that $\mathcal{L}_2(\mathcal{H}_{\mathcal{A}_{\ell}},\mathcal{H}_{\mathcal{B}_{\ell}})$ is an RKHS in its own right, for which we can assume smoothness in the sense of (ref) and spectral decay in the sense of (ref).
Just as we place approximation assumptions for $\gamma_0$ in terms of $\mathcal{H}$, which provides the features onto which we project $Y$, we place approximation assumptions for $E_{\ell}$ in terms of $\mathcal{H}_{\mathcal{B}_{\ell}}$, which provides the features $\phi(B_{\ell})$ onto which we project $\phi(A_{\ell})$. Under these conditions, we arrive at our main theoretical guarantee.
Explicit constants hidden by the $O(\cdot)$ notation, as well as explicit specializations of Assumption (ref), are indicated in Appendices (ref) and (ref). These rates approach $n^{-1/4}$ when $(c,c_1,c_2)=2$ and $(b,b_1,b_2)\rightarrow \infty$, i.e. when the regressions are smooth and when the effective dimensions are finite. Interestingly, each rate combines minimax optimal rates in RKHS norm: $n^{-(c-1)/\{2(c+1/b)\}}$ for standard nonparametric regression fischer2017sobolev; $\tilde{n}^{-1/2}$ for unconditional mean embeddings tolstikhin2017minimax; and, in contemporaneous work, $n^{-(c_{\ell}-1)/\{2(c_{\ell}+1/b_{\ell})\}}$ for conditional mean embeddings li2022optimal.
Overall, rates slower than $n^{-1/4}$ reflect the challenge of a $\sup$ norm guarantee, which is stronger than a mean square error guarantee and encodes caution about worst case scenarios when informing policy decisions. For comparison, the minimax optimal Sobolev norm rate for learning an $s_0$-smooth regression, using $\mathbb{H}_2^s$ over $\mathbb{R}^p$, is $n^{-(c-1)/\{2(c+1/b)\}}=n^{-(s_0-s)/(2s_0+p)}$.