Unconditional Quantile Regression with High Dimensional Data

Unconditional Quantile Regression with High-Dimensional Data

Introduction

Analysis of an outcome response to a counterfactual shift in the covariate distribution is of interest in policy studies. Such a counterfactual analysis requires accounting for the Oaxaca-Blinder decomposition of heterogeneous outcome distributions into structural heterogeneity ($F_{Y|X}$) and distributional heterogeneity ($F_X$); see fortin2011decomposition for a review. To conduct a credible counterfactual analysis, it is crucial to control a structure ($F_{Y|X}$) with rich information about $X$ while applying a counterfactual shift in the distribution of $X$. In this light, a researcher ideally wants to use high-dimensional $X$ in data.

Motivated by this feature of causal inference and the recently increasing availability of high-dimensional data, we develop a novel theory and method of estimation and inference for heterogeneous counterfactual effects with high-dimensional controls. The existing literature features a number of alternative approaches and frameworks of counterfactual analysis. Among others, we focus on the unconditional quantile partial effect \citep*[UQPE;][]{firpo2009unconditional} in the unconditional quantile regression based on the re-centered influence function (RIF) of firpo2009unconditional for two reasons: (i) its advantage of providing “a simple way of performing detailed decompositions” fortin2011decomposition and (ii) its popularity.\footnote{As of February 17, 2022, firpo2009unconditional have attracted 2275 Google Scholar citations.} This parameter measures the marginal effect of counterfactually shifting the distribution of a coordinate of $X$ on population quantiles of an outcome.

The UQPE is defined with the conditional distribution $F_{Y|X}$ and the marginal distribution $F_X$. Let $X=(X_1,X_{-1})$ denote the status quo, where $X_1$ is a scalar treatment variable of interest and $X_{-1}$ consists of controls. We focus on the change from $X$ to $(X_1+\varepsilon,X_{-1})$ throughout this paper, while our analysis can be generalized to the change in any fixed direction. The counterfactual distribution of $Y$ after this change is $$ F_{Y}^\varepsilon(y)=\int F_{Y|X=(x_1+\varepsilon,x_{-1})}(y)dF_X(x). $$ The UQPE with respect to the first coordinate, $X_1$, of $X$ is defined by

where $Q_\tau$ is the $\tau$-th quantile operator. The UQPE measures the change in the outcome quantile when the distribution of $X$ changes infinitesimally in the direction of the first coordinate. Under the assumption of conditional exogeneity, as in chernozhukov2013inference, the UQPE can be interpreted as the causal effect of changing the distribution of $X$ infinitesimally. Without such an assumption, UQPE may still be of interest as a summary statistic of the counterfactual distributional relationship between $Y$ and $X_1$.

To fix ideas, suppose that a policy maker is interested in analyzing the counterfactual effects of extending the duration $X_1$ of an exposure to the Job Corps training program on the outcome $Y$ of hourly wages, controlling for a large number of demographic, socioeconomic and behavioral attributes $X_{-1}$. The median of the actual distribution of hourly wages is given by $Q_{0.5}(F_Y^0)$. On the other hand, if the exposures are extended by $\varepsilon$ days for every participant, then the median of the counterfactual distribution of hourly wages becomes $Q_{0.5}(F_Y^\varepsilon)$. In this case, $Q_{0.5}(F_Y^\varepsilon) - Q_{0.5}(F_Y^0)$ measures the counterfactual effect on the median of the wage distribution, and $UQPE(\tau) = \lim_{\varepsilon \rightarrow 0} \big(Q_{0.5}(F_Y^\varepsilon) - Q_{0.5}(F_Y^0)\big)/\varepsilon$ measures its marginal effect.

While the RIF regression approach is indeed simpler to implement than alternative methods of counterfactual analysis as emphasized by fortin2011decomposition, an estimation of the UQPE still requires a three-step procedure. The first step is an estimation of unconditional quantiles. The second step implements the RIF regression. The third step integrates the RIF regression estimates to in turn estimate the UQPE. firpo2009unconditional provide an estimation procedure for the case of low-dimensional data. If we allow for high-dimensional controls with the aforementioned motivation, then the second step will require some estimation of the high-dimensional RIF regression, and the traditional techniques to incorporate estimation errors of the second step into the third step no longer apply. To overcome this challenge, we construct a novel doubly/locally robust score for estimation of the UQPE. The key insight for the construction is the identification result in firpo2009unconditional that the UQPE has the same structure as the average derivative estimator, whose influence function in the presence of nonparametric preliminary estimation has been well studied in the existing literature (e.g., newey1994asymptotic). With this doubly/locally robust score, we obtain a Z-estimation criterion with robustness against perturbations in functional nuisance parameters as in belloni2014uniform and belloni2018uniformly, and can thereby use the debiased estimation approach belloni2014uniform,chernozhukov2017double,chernozhukov2018double,chernozhukov2016locally, which allows one to obtain the asymptotic distribution of a UQPE estimator, independently of the second-step estimation as far as it satisfies some convergence rate conditions satisfied by major nonparametric estimators and machine learners. To provide a readily applicable method for practitioners, we focus on a specific method of estimation and bootstrap inference in the main text, but we also provide a generic method in the online supplement.

\noindentNotations. In this paper, we will use the following mathematical symbols and notations. $\mathcal{X}$ denotes the support of $X$. For a vector $v$, we define $\mathrm{Supp}(v)=\{i: v_i\ne 0\}$, $\|v\|_1=\sum_{i}|v_i|$, $\|v\|_2=(\sum_{i}v_i^2)^{1/2}$, and $\|v\|_\infty=\max_{i}|v_i|$. Denote the cardinality of $\mathrm{Supp}(v)$ by $\|v\|_0$. For a matrix $A$, we define $\|A\|_\infty = \max_{ij}|A_{ij}|$. We let $\Lambda$ be the standard logistic CDF and $\Phi$ be the standard normal CDF.

Robust Score, Estimation, and Inference for $UQPE(\tau)$

In this section, we develop a new score for doubly/locally robust estimation of the UQPE. We then present specific estimation and inference procedures in Sections (ref) and (ref), respectively. It is worthwhile to mention here that our analysis allows the dimensionality of $X$ to depend on the sample size $N$ and to diverge as $N\rightarrow\infty$.

Following firpo2009unconditional, we can rewrite our parameter of interest, defined in (ref), as a function of identifiable objects. Namely, $$ UQPE(\tau)=-\frac{\theta(\tau)}{f_Y(q_\tau)}, $$ where $f_Y$ is the density function of $Y$, $q_\tau$ is the $\tau$-th quantile of $Y$, and

This equation is shown in firpo2009unconditional.

Doubly/Locally Robust Score

We could estimate $\theta(\tau)$ based on (ref) and some estimator for $F_{Y\mid X}(\cdot)$. When $X$ is high-dimensional, this direct estimation of $\theta(\tau)$ can result in a large bias, a large variance, or both. Instead, we propose to construct an estimator for $\theta(\tau)$ based on an alternative representation:

where $\omega(x)=\partial\log f_{X_1|X_{-1}=x_{-1}}(x_1)/\partial x_1$, $m_0(x,q)=F_{Y\mid X=x}(q)$ and $m_1(x,q)={\partial m_0(x,q)}/{\partial x_1}$. This representation in (ref) comes from the influence adjustment term for the average derivative estimator newey1994asymptotic. Namely, $\int \omega(x)(1\{y\leq q_\tau\}-{m}_0(x,q_\tau))dF_{Y,X}(y,x)$ in (ref) adjusts the estimation error from the regularized preliminary estimation.

The advantage of (ref) over (ref) is that (ref) is doubly robust in the sense that

and

hold for a set of values that the high-dimensional nuisance parameters $(\tilde\omega(x),$ $\tilde{m}_{0}(x,q),$ $\tilde{m}_{1}(x,q))$ take as far as $\tilde{m}_1(x,q)={\partial\tilde{m}_0(x,q)}/{\partial x_1}$ and some regularity conditions to be formally stated below are satisfied. Note that $(\tilde\omega(x),\tilde{m}_{0}(x,q),\tilde{m}_{1}(x,q))$ in (ref) and (ref) can be different from the true value $(\omega(x),{m}_{0}(x,q),{m}_{1}(x,q))$. Our doubly/locally robust moment involves two nuisance parameters, $\omega(x)$ and $m_0(x,q)$, which are based on the conditional distribution of $X_1$ on $X_{-1}$ and the conditional distribution of $Y$ on $X$, respectively. They are analogous to the propensity score and the conditional mean function in the doubly robust moment for the average treatment effect. The equality in (ref) (resp. (ref)) implies that, even if we mis-specify the conditional distribution of $Y$ given $X$ (resp. the conditional distribution of $X_1$ given $X_{-1}$), (ref) provides the correct parameter value $\theta(\tau)$. The construction of our doubly/locally robust moment comes from the fact that ${\theta}(\tau)$ can be represented in two different ways: $$ {\theta}(\tau)=\int {m}_1(x,q_\tau)dF_X(x) \quad\mbox{ and }\quad {\theta}(\tau)=-\int \omega(x)1\{y\leq q_\tau\}dF_{Y,X}(y,x). $$ Each of the two representations only requires one of the two functions, ${m}_1(x,q_\tau)$ and $\omega(x)$, to be specified correctly. A precise statement for the double robustness and its proof are found in Appendix (ref) in the online supplement.

Estimation Procedure

With the sample $\{(Y_i,X_i): i=1,\ldots,N\}$ and the moment condition (ref), we propose to estimate ${\theta}(\tau)$ by a plug-in method. Let $(\hat\omega(x),\hat{m}_{0}(x,q),\hat{m}_{1}(x,q))$ denote an estimator of $(\omega(x),{m}_{0}(x,q),{m}_{1}(x,q))$ -- a concrete procedure to construct $(\hat\omega(x),\hat{m}_{0}(x,q),\hat{m}_{1}(x,q))$ is provided below. Letting $\hat{q}_\tau$ denote the sample $\tau$-th empirical quantile of $Y$, we estimate $\theta(\tau)$ by

With this estimator for $\theta(\tau)$, our proposed estimator for $UQPE(\tau)$ is in turn defined by $$ \widehat{UQPE}(\tau) = -\frac{\hat{ \theta}(\tau)}{\hat{f}_Y(\hat{q}_\tau)}, $$ where $\hat{f}_Y(y)$ is the kernel density estimator defined by $$ \hat{f}_Y(y) = \frac{1}{N}\sum_{i=1}^N\frac{1}{h_1}K_1\left(\frac{Y_i - y}{h_1}\right) $$ for a kernel function $K_1$ and a bandwidth parameter $h_1$.

We use the logistic Lasso regression BCFH13 to construct $\hat{m}_{0}(x,q)$. Once we construct $\hat{m}_{0}(x,q)$, we in turn define $\hat{m}_{1}(x,q)$ by $$ \hat{m}_{1}(x,q) = \frac{\partial\hat{m}_{0}(x,q)}{\partial x_1}. $$ Consider the approximately sparse logistic regression model for $m_0(x,q)$: $$ m_0(X,q) = \Lambda(b(X)^\top\beta_{q})+(\mathrm{approximation\ error}), $$ where $b(X)$ is a $p_b$-dimensional observed vector and $\beta_{q}$ is an unknown parameter. Assumption (ref) (to be stated below in Section (ref)) specifies the conditions for the sparsity and formalizes the approximation error. In our numerical examples, we define $b(X)$ by including powers of $X$ up to the third degree and standardize each component of $b(X)$ so that the variance is one. We estimate $\beta_{q}$ by the Lasso penalized logistic regression

where ${\Psi}_q$ is a diagonal matrix with penalty loadings defined in the next paragraph. We follow BCFH13 and set the regularization parameter as $$ \lambda_{L} = 1.1\Phi^{-1}(1-(0.1/\log(N))/(p_b\vee N))N^{1/2}. $$ We recommend using the post-Lasso estimator for $\beta_q$ defined by $$ \hat{\beta}_{q} = \operatorname*{arg\,min}_{\beta \in \mathbb{R}^p: \mathrm{Supp}(\beta)\subset(\mathrm{Supp}(\tilde{\beta}_{q}) \cup S_1)}-\frac{1}{N}\sum_{i=1}^N\log\left(\Lambda(b(X_i)^\top\beta)^{1\{Y_{i} \leq q\}}(1-\Lambda(b(X_i)^\top\beta))^{1\{Y_{i}>q\}}\right), $$ where $S_1 \subset \{1,\ldots,p_b\}$ denotes the coordinate set of covariates researchers want to include in the post-Lasso regression. For the UQPE with respect to $X_1$, it is natural to include $X_1$ in the regression. The post-Lasso estimator can do so by setting $1\in S_1$, whereas the Lasso estimator $\tilde{\beta}_{q}$ may exclude $X_1$ from the regression. With $\hat{\beta}_{q}$, we can estimate $\hat{m}_{0}(x,q)$ by $$ \hat{m}_{0}(x,q) = \Lambda(b(x)^\top \hat{\beta}_{q}). $$

The penalty loading matrix ${\Psi}_q=\text{diag}(\psi_{q,1},\cdots ,\psi_{q,p_b})$ needs to be estimated to implement (ref). Ideally, we would like to use the infeasible penalty loading $$ \bar{\psi}_{q,j}=\sqrt{\frac{1}{N}\sum_{i=1}^N\left(1\{Y_i\leq q\}-m_{0}(X_i,q)\right)^2b^2_{j}(X_i)}. $$ Since $m_{0}(X,q)$ is unknown, BCFH13 propose the following iterative algorithm to obtain the feasible version of the loading matrix:

Next, we consider a regularized estimation of $\omega(X)$ based on the Riesz representer approach chernozhukov2018automatic,chernozhukov2018double2. Suppose that $h(x)$ is a $p_h$-dimensional dictionary of approximating functions that are differentiable in $x_1$ and that $$ \omega(x)=h(x)^\top\overline{\rho}+(\mathrm{approximation\ error}) $$ holds, where Assumption (ref) (to be stated below in Section (ref)) formally describes this approximation. In our numerical examples, we define $h(X)$ by including powers of $X$ up to the third degree and standardize each component of $h(X)$ so that the variance is one. Since $\omega(x)={\partial \log f_{X_1|X_{-1}=x_{-1}}(x_1)}/{\partial x_1}$, the integration by parts yields $$ \mathbb{E}[h(X) \omega(X)] = - \mathbb{E}[\partial_{x_1}h(X)]. $$ Approximating $\omega(x)$ by $h(x)^\top\overline{\rho}$, we have $$ \mathbb{E}[h(X)h(X)^\top]\overline{\rho}=- \mathbb{E}[\partial_{x_1}h(X)]+(\mathrm{approximation\ error}). $$ Thus, $\overline{\rho}$ can be approximated by $\operatorname*{arg\,min}_{\rho}\left(-2M^\top \rho + \rho^\top G\rho\right)$, where $G = \mathbb{E}[h(X)h(X)^\top]$ and $M = -\mathbb{E}[\partial_{x_1}h(X)]$. To accommodate high-dimensional $h(x)$, we use the regularized minimizer $$ \operatorname*{arg\,min}_{\rho}\left(-2M^\top \rho + \rho^\top G\rho + \lambda_{R} \|\rho\|_1\right) $$ with $\lambda_{R}$ denoting a regularization parameter (cf. Assumption (ref).5). In the simulations and empirical application, we use $\lambda_{R} = 2\log(\log(N)) \sqrt{\log(p_h)/N}$. The Riesz representer approach uses the sample analog of this objective to estimate $\omega(x)$. Namely, we estimate $\omega(x)$ by $$\hat\omega(x) = h(x)^\top \hat{\rho},$$ where $\hat{G} = \frac{1}{N}\sum_{i=1}^Nh(X_i)h(X_i)^\top$, $\hat{M} = -\frac{1}{N}\sum_{i=1}^N\partial_{x_1}h(X_i)$, and $$ \hat{\rho} = \operatorname*{arg\,min}_{\rho}\left(-2 \hat{M}^\top \rho + \rho^\top \hat{G} \rho + \lambda_{R} \|\rho\|_1\right). $$

Bootstrap Inference

For an inference about $UQPE(\tau)$, we propose the multiplier bootstrap without requiring to recalculate the preliminary estimator $\hat\omega(x)$ in each bootstrap iteration. (More precisely, if we calculate $(\hat m_0(x,q), \hat m_1(x,q))$ on a grid of values of $q$ once, then we do not need to recalculate them in each bootstrap iteration either.) Using independent standard normal random variables $\{\eta_i\}_{i=1}^N$ that are independent of the data, we compute the bootstrap estimators $\hat{ \theta}^*(\tau)$ and $\widehat{UQPE}^*(\tau)$ in the following steps. The bootstrap estimator for $q_\tau$ is $\hat{q}^*_\tau$ defined by the $r_N^*$-th order statistic of $Y_i$, where $r_N^*$ is the integer part of $1+\sum_{i=1}^N\left(\tau+\eta_i(\tau - \mathbf{1}\{Y_i \leq \hat{q}_\tau\})\right)$.\footnote{It is equivalent to $\hat{q}^*_\tau = \arg\min_q \sum_{i=1}^N \rho_\tau(Y_i - q) - q \sum_{i=1}^N \eta_i(\tau - 1\{Y_i \leq \hat{q}_\tau\})$, which is the gradient bootstrap method chen2004quantile and directly perturbs the score for the quantile $q_\tau$. By the sub-gradient condition, we have that $\hat{q}^*_\tau$ equals the $r_N^*$th order statistic of $Y_i$, where $r_N^*$ is the integer that satisfies $N\tau + \sum_{i=1}^N \eta_i(\tau - 1\{Y_i \leq \hat{q}_\tau\}) + 1 \geq r_N^* \geq N\tau + \sum_{i=1}^N \eta_i(\tau - 1\{Y_i \leq \hat{q}_\tau\})$. This procedure gives us a closed-form solution for $\hat{q}^*_\tau$.} The bootstrap estimators for $f_Y(y)$ and $\theta(\tau)$ are $$ \hat{f}_Y^*(y) = \frac{1}{\sum_{i=1}^N(\eta_i+1)}\sum_{i=1}^N(\eta_i+1)\frac{1}{h_1}K_1\left(\frac{Y_i - y}{h_1}\right), $$ and $$ \hat{\theta}^*(\tau)= \frac{1}{\sum_{i=1}^N(\eta_i+1)}\sum_{i=1}^N (\eta_i+1)\left(\hat{m}_{1}(X_i,\hat{q}^*_\tau) - \hat\omega(X_i)(1\{Y_i \leq \hat{q}^*_\tau\} - \hat{m}_{0}(X_i,\hat{q}^*_\tau))\right), $$ respectively. With these components, the bootstrap estimator $\widehat{UQPE}^*(\tau)$ is given by $$ \widehat{UQPE}^*(\tau) = -\frac{\hat{ \theta}^*(\tau)}{\hat{f}^*_Y(\hat{q}^*_\tau)}. $$

We can use the above multiplier bootstrap method to conduct various types of statistical inference about the UQPE. First, testing the hypothesis of $UQPE(\tau) = 0, \forall\tau\in \Upsilon$ for some closed interval $\Upsilon\subset(0,1)$ is of main interest in many empirical applications. Because $f_Y(q_\tau)$ is assumed to be bounded away from zero, such a hypothesis is equivalent to $\theta(\tau) = 0, \forall\tau \in \Upsilon$, where $\theta(\tau)$ can be estimated in a parametric rate. We can thus test $UQPE(\tau) = 0, \forall\tau\in \Upsilon$ by constructing a confidence band for $\{\theta(\tau): \tau \in \Upsilon \}$ and checking whether the constant zero function belongs to this band. Specifically, let $$ \hat{\sigma}^{\theta}(\tau) = \frac{Q_{\hat\theta^*(\tau)}(0.75) - Q_{\hat\theta^*(\tau)}(0.25)}{\Phi^{-1}(0.75)-\Phi^{-1}(0.25)} $$ denote an estimator of the standard error of $\hat{\theta}(\tau)$ for $\tau \in \Upsilon$, where $Q_{\hat\theta^*(\tau)}(0.75)$ and $Q_{\hat\theta^*(\tau)}(0.25)$ denote the 75th and 25th percentiles of $\hat\theta^*(\tau)$ conditional on the data. Let $c_{\Upsilon}^{\theta}(1-\alpha)$ denote the $(1-\alpha)$ quantile of $$ \sup_{\tau \in \Upsilon} \left\vert \frac{\hat{\theta}^*(\tau)-\hat{\theta}(\tau)}{\hat\sigma^{\theta}(\tau)} \right\vert $$ conditional on the data. Let $CB^{\theta}_{\Upsilon}$ denote the confidence band of $\theta(\cdot)$ on $\Upsilon$ whose lower and upper bounds at $\tau \in \Upsilon$ are given by $\hat{\theta}(\tau) \pm \hat\sigma^{\theta}(\tau) c_{\Upsilon}^\theta(1-\alpha)$.

Second, we can similarly construct a confidence band for $\{UQPE(\tau):\tau\in \Upsilon\}$. Let $$ \hat\sigma(\tau) = \frac{Q_{\widehat{UQPE}^*(\tau)}(0.75) - Q_{\widehat{UQPE}^*(\tau)}(0.25)}{\Phi^{-1}(0.75)-\Phi^{-1}(0.25)} $$ denote an estimator of the standard error of $\widehat{UQPE}(\tau)$ for $\tau \in \Upsilon$, where $Q_{\widehat{UQPE}^*(\tau)}(0.75)$ and $Q_{\widehat{UQPE}^*(\tau)}(0.25)$ denote the 75th and 25th percentiles of $\widehat{UQPE}^*(\tau)$ conditional on the data. Let $c_\Upsilon(1-\alpha)$ denote the $(1-\alpha)$ quantile of $$ \sup_{\tau \in \Upsilon} \left\vert \frac{\widehat{UQPE}^*(\tau)-\widehat{UQPE}(\tau)}{\hat\sigma(\tau)}\right\vert $$ conditional on the data. Let $CB_\Upsilon$ denote the confidence band of $UQPE(\cdot)$ on $\Upsilon$ whose lower and upper bounds at $\tau \in \Upsilon$ are given by $\widehat{UQPE}(\tau) \pm \hat\sigma(\tau) c_\Upsilon(1-\alpha)$.

Asymptotic Theory

In this section, we investigate the asymptotic properties of the estimators $(\hat{ \theta},\widehat{UQPE})$ and their bootstrap counterparts $(\hat{ \theta}^*,\widehat{UQPE}^*)$ introduced in the previous section. The uniformity over $\tau$ is relevant to applications (e.g., analysis of heterogeneous counterfactual effects across $\tau$), and therefore, in this section, we aim to control the residuals for the linear expansion uniformly over $\tau\in \Upsilon$ for some closed interval $\Upsilon\subset(0,1)$. Let $\mathcal{Q} = \{{q}_\tau: \tau \in \Upsilon\}$, and let $\mathcal{Q}^\delta$ denote the $\delta$ enlargement of $\mathcal{Q}$.

This assumption is on the model primitives. Assumptions (ref).1 and (ref).3--(ref).6 impose regularity in terms of the smoothness of various functions representing features of the data. We use the enlargement $\mathcal{Q}^\delta$ instead of $\mathcal{Q}$ because with probability approaching one, $\hat{q}_\tau$ and $\hat{q}^*_\tau$ belong to the former but not necessarily the latter. Assumption (ref).2 is the standard moment condition.

We impose the following condition to bound the estimation error for $\hat{m}_{j}(x,q)$.

Several remarks are in order. First, Assumption (ref).1 is the common regularity condition. Second, Assumptions (ref).2--(ref).4 are common in the literature of logistic regressions with an $\ell_1$ penalty. See, for instance, BCFH13. Third, the various bounds for $m_1(x,q)$, $\beta_{q}$, $(m_0(X,q)-\Lambda(b(X)^\top\beta_{q}))$ need to hold uniformly over $q \in \mathcal{Q}^\delta$ because the estimator $\{\hat{q}_\tau: \tau \in \Upsilon\}$ belongs to $\mathcal{Q}^\delta$ for any fixed $\delta$ with probability approaching one. Fourth, as formally stated in Theorem (ref) in the Online Supplement, Assumption (ref) can bound the estimation error for the logistic Lasso estimation: $$ \sup_{q \in \mathcal{Q}^\delta}\int\left|\hat{m}_{j}(x,q)-{m}_j(x,q)\right|^2dF_X(x) = O_P\left(\frac{s_b \log(p_b)}{N}\right) $$ and $$ \sup_{q \in \mathcal{Q}^\delta,x \in \mathcal{X}}\left| \hat{m}_{j}(x,q) - m_j(x,q)\right|= O_P\left( \zeta_N s_b\sqrt{\frac{\log(p_b)}{N}}\right). $$

We provide the regularity condition for the Riesz representer estimation of $\omega(x)$.

Assumption (ref) follows chernozhukov2018automatic to which we refer readers for more discussion. Specifically, we have their $\varepsilon_n=\sqrt{ \log(p_h)/N}$ and $r = \lambda_{R}$ and Assumption (ref).4 implies chernozhukov2018automatic by BRT09. Theorem (ref) in the Online Supplement shows that Assumption (ref) can bound the estimation error for $\hat\omega(x)$: $$ \int \left(\hat\omega(x) - \omega(x)\right)^2dF_X(x) = o_P(N^c s_h \log(p_h)/N) \mbox{ and } \sup_{x \in \mathcal{X}}|\hat\omega(x) - \omega(x)| = o_P(1). $$ for all $c>0$. In Theorem (ref), we also show the Riesz representer $\hat\omega(x)$ belongs to a class of functions whose entropy or complexity level is well-controlled. Such a result is new to the literature and essential for our theory as we use all the observations to estimate $\omega(x)$, and thus, are subject to the model selection bias. chernozhukov2018automatic circumvent such bias via cross-fitting. In the Online Supplement, we also consider cross-fitting which can accommodate complicated general machine learning estimators for $\omega(x)$.

Testing $UQPE(\tau)=0, \forall\tau \in \Upsilon$

As mentioned earlier, testing $UQPE(\tau)=0, \forall \tau \in \Upsilon$ is equivalent to testing $\theta(\tau)=0, \forall\tau \in \Upsilon$. We can reject the null hypothesis if the constant zero function over $\Upsilon$ does not belong to $CB^\theta_\Upsilon$. In this section, we show that the proposed confidence band $\text{CB}^\theta_\Upsilon$ covers the true $\theta(\tau)$ uniformly over $\tau \in \Upsilon$ with the correct asymptotic size.

We impose an additional rate condition about upper bounds on $s_h$ and $s_b$.

When $\omega(\cdot)$ defined in Assumption (ref) is bounded so that $d = \infty$ and $\zeta_N$ defined in Assumption (ref) is also bounded, $\pi_N$ is roughly equal to $\sqrt{N^{2c}s_h \log(p_h)/N + s_b \log(p_b)/N}$, which is just the convergence rate for the first-stage estimators. In this case, this additional condition holds as long as $s_h \log(p_h) + s_b \log(p_b)= o(N^{1/2-c})$ for some $c>0$, which implies $$ \sqrt{\frac{s_h \log (p_h)}{N}} + \sqrt{\frac{s_b \log (p_b)}{N}} = o(N^{-1/4}). $$ It means the convergence rate of the nuisance functions should be faster than the rate of $N^{1/4}$. Such a rate is sufficient for the influence function representation for $\hat{ \theta}(\tau)$ and $\hat{ \theta}^*(\tau)$ (in Theorem (ref)) due to the use of doubly/locally robust moment.

Theorem (ref) provides a sufficient condition for the correct asymptotic size of the proposed confidence band $\text{CB}^\theta_\Upsilon$. It follows as a corollary of Theorem (ref) in the Online Supplement.

Confidence Band for $\{UQPE(\tau): \tau \in \Upsilon\}$

In this section, we consider the confidence band for $\{UQPE(\tau): \tau \in \Upsilon\}$, which can be used to infer the entire trajectory of $UQPE(\tau)$ over $\tau \in \Upsilon$. Recall that we use the kernel function $K_1(\cdot)$ in the kernel density estimation of $f_Y(\cdot)$. We impose the following conditions for the kernel function $K_1(\cdot)$ and the bandwidth parameter $h_1$.

We impose an additional rate condition about upper bounds on $s_h$ and $s_b$.

This assumption is weaker than Assumption (ref), as $\log(N)h_1=o(1)$. In other words, if $d$ is sufficiently large and $\zeta_N$ is bounded, such a condition holds as long as $\sqrt{\log(N)h_1}(s_h \log(p_h) + s_b \log(p_b))= o(N^{1/2-c})$ for some $c>0$.

The following theorem summarizes the validity for the bootstrap inference. The main takeaway from this theorem is that, by using the doubly robust method, we greatly relax the requirements on the sparsity and the number of effective covariates for $m_0(x,q)$ and $m_1(x,q)$, at the cost of imposing a sparsity condition on $\omega(x)$.\footnote{The leading term of the score function is $\frac{\theta(\tau)}{f_Y^2({q}_\tau)h_1}K_1\left(\frac{Y_i - {q}_\tau}{h_1}\right)$, but it does not imply that the doubly robust estimation for $\theta(\tau)$ is unnecessary. In fact, $\sup_{\tau \in \Upsilon}|\hat{ \theta}(\tau) - \theta(\tau)|$ is asymptotically negligible compared to $\frac{\theta(\tau)}{f_Y^2({q}_\tau)h_1}K_1\left(\frac{Y_i - {q}_\tau}{h_1}\right)$ if $\pi_N(\sqrt{s_h \log(p_h)/N} + \sqrt{s_b \log(p_b)/N}) + N^{-(1+d)/(2+d)}(s_h \log(p_h) + s_b \log(p_b))=o((\log(N)Nh_1)^{-1/2})$. If $h_1 = N^{-1/5}$, $d$ is sufficiently large, $c$ is arbitrarily small, and $\zeta_N$ is polylogarithmic, such a condition holds if $s_h \log(p_h) + s_b \log(p_b)= o(N^{3/5})$ up to some polylogarithmic factor. On the other hand, if we do not use the doubly robust method, the estimation error of $\theta(\tau)$ is $\sqrt{s_b\log(p_b)/N}$, which is asymptotically negligible if $\sqrt{s_b\log(p_b)/N} = o((\log(N)Nh_1)^{-1/2})$. Such a condition would require $s_h \log(p_h) = o(N^{1/5})$ up to some polylogarithmic factor.}

Theorem (ref) is a direct consequence of the linear expansions for $\widehat{UQPE}(\tau)$ and $\widehat{UQPE}^*(\tau)$ (formally stated in Theorem (ref) in the Online Supplement) and the strong approximation theory developed by CCK14-anti,CCK14. To compute $\hat\sigma(\tau)$, we can use either the plug-in method or the bootstrap method. For these methods, the convergence rate of $\sqrt{Nh_1}\hat\sigma(\tau)$ is polynomial in $N$, which implies $o_P(\log^{-1/2}(N))$.

Simulation Studies

In this section, we use Monte Carlo simulations to study the finite sample performance of the proposed method of estimation and inference for the UQPE.

Consider the following set of alternative data-generating designs. The outcome variable is generated according to the partial linear high-dimensional model $$ Y\mid X \sim N\left(g(X_1) + \sum_{j=2}^p \alpha_j X_j, \ \ 1\right), $$ where the function $g(x)$ is defined in the following three ways: $g(x) = x$ in DGP 1, $g(x) = x - 0.10 \cdot x^2$ in DGP 2, and $g(x) = x - 0.10 \cdot x^2 + 0.01 \cdot x^3$ in DGP 3. The high-dimensional controls $(X_1,...,X_p)^\top$ are generated by $$ X_1\mid (X_2,...,X_p) \sim N\left(\sum_{j=2}^p \gamma_j X_j ,\ 1\right) \mbox{\ \ and\ \ }(X_2,...,X_p) \sim N(0,\Sigma_{p-1}), $$ where $\Sigma_{p-1}$ is the $(p-1) \times (p-1)$ variance-covariance matrix whose $(r,c)$-element is $0.5^{2(|r-c|+1)}$. Note that this data-generating process induces dependence of the control $X_1$ of main interest on the rest of the $p-1$ controls $(X_2,\ldots,X_p)^\top$, as well as the dependence among the $p-1$ controls $(X_2,\ldots,X_p)^\top$. For the high-dimensional parameter vectors in the above data-generating model, we consider the following four cases of varying sparsity levels:

We follow the general estimation and inference approach outlined in Section (ref). We set $h(x) = (x^\top,(x^2)^\top,(x^3)^\top)^\top$ for estimation of $\omega(x)$, and set $b(x) = (x^\top,(x^2)^\top,(x^3)^\top)^\top$ for estimation of $m_0$ and $m_1$. For the choice of $h_1$, we under-smooth the rule-of-thumb optimal choice as $h_1 = 1.06\sigma(Y) N^{-1/5-0.01}$. For each design, we use 500 iterations of Monte Carlo simulations to compute the mean, bias, and root mean square error (RMSE) of the estimate, as well as the 95% uniform coverage over the set $[0.20,0.80]$ of quantiles. To evaluate the bias, RMSE, and the 95% uniform coverage, we first numerically approximate the true UQPE by large-sample Monte Carlo simulations. Across sets of Monte Carlo simulations, we vary the DGP $\in \{\text{DGP 1, DGP 2, DGP 3}\}$ and the sparsity design $\in \{\text{(i)},\text{(ii)},\text{(iii)},\text{(iv)}\}$, while we fix the sample size $N = 500$ and the dimension $p = 100$ throughout.

Table (ref) summarizes the simulation results under the sparsity designs (i) and (ii). We can make the following three observations in these results. First, the bias of our UQPE estimator is small, especially relative to the RMSE. This feature of the results supports the fact that our estimator mitigates the bias via the use of the doubly robust score. Second, the RMSE decreases as the sample size increases. Third, the 95% uniform coverage frequencies are close to the nominal probability, namely, 0.95. This feature of the results supports our theory on the asymptotic validity of the bootstrap inference. From these simulation results, we confirm the main theoretical properties of the proposed method of estimation and inference for the UQPE across alternative data-generating processes. Table (ref) shows the simulation results under the less sparse designs (iii) and (iv). While the bias and RMSE are slightly bigger here than those in Table (ref), the magnitudes of changes are modest. In addition to the simulation designs introduced above, we also experimented with other designs, and the simulation results are very similar and support the main theoretical properties of our proposed method as well -- see Appendix (ref) in the online supplement.

To highlight the value added by our proposed method to the existing literature, we also experiment with the RIF-Logit estimator from firpo2009unconditional as a benchmark. Table (ref) summarizes the simulation results based on the RIF-Logit for the sample size of $N=500$ under the most sparse design (i). Observe that, as the dimension $p$ increases from 25 to 50, the finite sample performance substantially degrades in terms of all of the displayed statistics, namely the bias, RMSE, and (pointwise and uniform) 95% coverage frequencies. In particular, the uniform coverage frequency drops to zero even for the dimension that is as small as $p=50$. With the same sample size of $N=500$, on the other hand, our proposed method produces accurate coverage frequencies as well as accurate estimates for the even larger dimension $p=100$ as presented in Tables (ref) and (ref). This comparison sheds lights on favorable finite sample performance of our proposed method when there are high-dimensional controls, in comparison with the existing alternative method. With that said, we would like to remark that we pay the costs of additional assumptions for nuisance parameter estimation, and hence there are tradeoffs between the existing procedure firpo2009unconditional and our proposed method.

In addition to the RIF-Logit estimator, we also experiment with a RIF-Lasso-Logit estimator. This estimator has not been formally investigated in the literature to our knowledge, but it coincides with a non-orthogonalized version of our procedure and hence serves as a useful benchmark to evaluate the benefits of our proposed doubly robust score. Table (ref) shows the simulation results based on this procedure without the doubly robust score (right panel) compared with the results of our proposed procedure with the doubly robust score (left panel) copied from Table (ref). While the coverage frequencies for our proposed method achieves the nominal probability of 95%, those for the counterpart without the doubly robust score fall short of 95%. These results show that the method without the doubly robust score incurs larger size distortions and demonstrate that it is useful to employ the doubly robust score as we do for our proposed method.

Finally, we consider the pointwise and uniform tests for $UQPE(\tau) =0$ by testing $\theta(\tau) = 0$. The corresponding simulation results are collected in Appendix (ref) in the online supplement.

Heterogeneous Counterfactual Marginal Effects of Job Corps Training

The UQPE identifies counterfactual effects that are heterogeneous across outcome levels $Y$. This feature of the UQPE is useful for evaluating economic policies designed to benefit targeted subpopulations of the economy that are identified in terms of economic outcomes such as wage and income. For instance, major job training programs are designed to benefit targeted subpopulations of individuals who are low wage earners, i.e., lower quantiles of $Y$. In redesigning a job training program, a policy maker may want to choose such changes in $X$ that particularly benefit these targeted subpopulations (with potentially lower wages) rather than the others (with potentially higher wages). Therefore, it is important for the policy maker to understand heterogeneous outcome gains (e.g., wage increase) of alternative counterfactual changes in $X$ across different subpopulations characterized by the levels of $Y$. The UQPE provides solutions to this goal.

While a rich set of empirical findings have been reported about the treatment effects of Job Corps, an analysis of heterogeneous counterfactual effects is missing in the literature to the best of our knowledge, despite its relevance to designing effective program policies and schemes as emphasized in the previous paragraph. Applying our proposed method, we analyze heterogeneous counterfactual marginal effects of Job Corps training on labor outcomes in this section. Specifically, it is important to find whether higher (respectively, lower) potential earners would benefit more (respectively, less) from counterfactually extending the duration of the training program. Since the entrance interview in Job Corps provides some information regarding the human capital of prospective trainees, answers to these empirical questions may possibly help the program designers to devise more efficient policies and schemes for the training programs. As such, we are interested in heterogeneous counterfactual marginal effects of the duration of the exposure to the program, as a continuous treatment variable, on labor outcomes measured by hourly wages.

We identify and estimate the counterfactual distributional change given a large set of observed controls by taking advantage of our machine-learning-based method. For the outcome variable, we consider hourly wages. For the continuous treatment variables, we consider two seemingly similar but different measures: the duration in days of participation in Job Corps and the duration in days of actually taking classes in Job Corps. As will be shown shortly, these two definitions lead to qualitatively different empirical findings. We use 42 observed controls (and their powers). Table (ref) shows the summary statistics of our data. Different sets of observations are missing across different variables, and hence we use the intersection of observations that are non-missing across all the variables in use for our analysis. After dropping the missing observations, we are left with $N=481$ when we define the duration of participation in Job Corps as the treatment, while we are left with $N=368$ when we define the duration of actually taking classes in Job Corps. Taking the intersection of these two samples, we use a subsample of size $N=347$. Note that the dimension of covariates is relatively large given this effective sample size, and hence high-dimensional econometric methods are indispensable in the current application.

Observe that our sample consists of high-dimensional controls and the sample size that results from the aforementioned sample selection is not sufficiently large for conventional econometric methods of estimation and inference for the UQPE. We therefore use our proposed method of estimation and inference for the UQPE that can accommodate a large dimension of controls via the use of the doubly robust score.

Using the same computer program as the one used for simulation studies presented in Section (ref), we obtain estimates, pointwise 95% confidence intervals, and uniform 95% confidence bands for $UQPE(\tau)$ for $\tau \in [0.20,0.80]$. Table (ref) summarizes the results. The row groups (I) and (II) report results for days in Job Corps as the treatment variable, while the row groups (III) and (IV) report results for days of taking classes in Job Corps as the treatment variable. The row groups (I) and (III) report results for the hourly wage as the outcome variable, while the row groups (II) and (IV) report results for the logarithm of the hourly wage as the outcome variable.

Overall, the magnitudes of the estimates are consistent with those from prior studies, and we also obtain the following new findings.\footnote{In the row group (I) in Table (ref) for instance, the daily marginal effects range from 0.0012 to 0.0046 dollars. This magnitude is consistent with the 0.22 difference in average hourly wages between the treatment and control groups \citep*[][Table 3]{SBM2008}, where the average number of days in Job Corps for the treated group is 153.4 (Table (ref)).} First, observe that the none of the uniform 95% confidence bands are contained in the negative reals. These results indicate that the counterfactual marginal effects of our interest are significantly negative for none of the heterogeneous subpopulations. We next look into the heterogeneity of these effects. Observe in row (I) that the uniform 95% confidence band for $\tau=0.2$ is contained in the positive reals while the uniform 95% confidence band for $\tau=0.8$ intersects with the zero. These results imply heterogeneous statistical significance across quantiles. Specifically, we predict significantly positive counterfactual effects for lower wage earners ($\tau=0.2$, $0.4$ and $0.6$) and insignificant effects for higher wage earners ($\tau=0.8$). On the other hand, the point estimate is smaller for $\tau=0.2$ than that for $\tau=0.8$ in row (I). The larger effects for the subpopulation of higher potential earners (i.e., higher quantiles) could simply result from the scale effect. Heterogeneity in causal effects across different quantiles often vanishes once we take the logarithm of the outcome variable.\footnote{The relationship ${\partial Q_{\tau}(F^\epsilon_{\log(Y)})}/{\partial \epsilon}\big|_{\epsilon=0} = {(\partial Q_{\tau}(F^\epsilon_{Y})/\partial \epsilon)}/{Q_{\tau}(F^\epsilon_{Y})}\big|_{\epsilon=0}$ implies that the sign of the level and the logarithm coincide at the population level, even though the signs of empirical estimates and statistical significance may not coincide as in our empirical results.} Therefore, we next consider the row group (II), where the outcome variable is defined as the logarithm of the hourly wage. Notice that, even in this row group, we continue to observe the same qualitative pattern as that in the row group (I). Namely, the point estimate is smaller for $\tau=0.2$ than that for $\tau=0.8$, but the uniform 95% confidence band for $\tau=0.2$ is contained in the positive reals while the uniform 95% confidence band for $\tau=0.8$ intersects with the zero. These results provide policy makers of confidence that extending the duration of exposures to the Job Corps program will benefit potential lower wage earners.

Once we turn to row groups (III) and (IV), where the treatment variable is now defined as days of taking classes in Job Corps, we no longer observe the aforementioned pattern of heterogeneous counterfactual marginal effects, and the uniform 95% confidence bands globally intersect with the zero. However, if we implement the test $UQPE(\tau)=0, \forall\tau\in [0.20,0.80]$ based on $\hat\theta(\tau)$ as in Section (ref), then we actually reject this hypothesis of uniformly zero counterfactual marginal effects with the 95% confidence.

In summary, we obtain the following three new findings about counterfactual marginal effects of the duration of exposure to Job Corps training on the hourly wage. First, the effects are significantly negative for none of the heterogeneous subpopulations under consideration, regardless of the definition of the treatment variable and the definition of the outcome variable. Second, the counterfactual marginal effects of days in Job Corps are significant for the subpopulation of lower wage earners, while they are insignificant for higher wage earners. This result holds robustly regardless of whether we define the outcome variable as the hourly wage or the logarithm of it. Third, we fail to detect the aforementioned pattern of counterfactual marginal effects once we define the treatment variable as days of taking classes in Job Corps, while we still reject the hypothesis of uniformly zero counterfactual marginal effects. These results contain the following policy implications. Extending the duration of exposures the Job Corps training program will be effective especially for the targeted subpopulations of lower potential wage earners. However, these benefits may come from sources other than the experience of taking classes in the Job Corps training program.

Finally, to get further insights about our proposed method, we conclude this section with discussions of more details about what was implemented in the black box to produce the results reported in Table (ref). When estimating the core functions $m_0(x,q_\tau)$ and $m_1(x,q_\tau)$ by the lasso logit, we in fact select different subvectors of $b(X)$ across different quantiles $\tau$. Table (ref) shows which controls and/or their powers are selected for each $\tau \in \{0.20,0.40,0.60,0.80\}$. While there are some controls that are common (such as the intercept and the number of people in household) across all $\tau$, the selections are by no means uniform across $\tau$. A researcher does not ex ante know which variables among many in the list should be included for each quantile $\tau$. Including all the potentially relevant controls would incur non-trivial size distortions -- recall the simulation results shown in Table (ref) in Section (ref). By our proposed method of inference that achieves the nominal size, on the other hand, it is those shown in Table (ref) that were selected to be relevant for each $\tau \in \{0.20,0.40,0.60,0.80\}$ in producing the estimation results reported in Table (ref).

Conclusion

Counterfactual analyses often involve high-dimensional controls. On the other hand, existing methods of estimation and inference for heterogeneous counterfactual changes are not compatible with high-dimensional settings. In this paper, we therefore propose a novel doubly/locally robust score for debiased estimation and inference for the UQPE as a measure of heterogeneous counterfactual marginal effects. A concrete implementation procedure is provide for estimation and multiplier bootstrap inference. The online supplement additionally presents a general class of estimation and inference procedures. Asymptotic theories are presented to guarantee that the bootstrap method works for size control. Simulation studies support our theoretical properties. Applying the proposed method of estimation and inference to survey data of Job Corps, the largest training program for disadvantaged youth in the United States, we obtain the following two policy implications. First, extending the duration of exposures the Job Corps training program will be effective especially for the targeted subpopulations of lower potential wage earners. Second, these benefits may come from sources other than the experience of taking classes in the Job Corps training program.