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.
94,062 characters · 18 sections · 73 citation commands
Model Selection Results for the International Growth Regressions
Quantile regression is an important statistical method for analyzing the impact of regressors on the conditional distribution of a response variable (cf. Laplace, KB78). It captures the heterogeneous impact of regressors on different parts of the distribution Buchinsky1994, exhibits robustness to outliers K2005, has excellent computational properties PornoyKoenker1997, and has wide applicability K2005. The asymptotic theory for quantile regression has been developed under both a fixed number of regressors and an increasing number of regressors. The asymptotic theory under a fixed number of regressors is given in KB78, Portnoy:mess, GJ1992, Knight1998, C2005 and others. The asymptotic theory under an increasing number of regressors is given in HS00 and BC-MCMC,BC-qrID, covering the case where the number of regressors $p$ is negligible relative to the sample size $n$ (i.e., $p=o(n)$).
In this paper, we consider quantile regression in high-dimensional sparse models (HDSMs). In such models, the overall number of regressors $p$ is very large, possibly much larger than the sample size $n$. However, the number of significant regressors for each conditional quantile of interest \comment{-- those having a non-zero impact on the response variable --} is at most $s$, which is smaller than the sample size, that is, $s = o(n)$. HDSMs (CandesTao2007,MY2007,BickelRitovTsybakov2009) have emerged to deal with many new applications arising in biometrics, signal processing, machine learning, econometrics, and other areas of data analysis where high-dimensional data sets have become widely available.
A number of papers have begun to investigate estimation of HDSMs, focusing primarily on penalized mean regression, with the $\ell_1$-norm acting as a penalty function BickelRitovTsybakov2009,CandesTao2007,Koltchinskii2009,MY2007,vdGeer,ZhangHuang2006. BickelRitovTsybakov2009,CandesTao2007,Koltchinskii2009,MY2007,ZhangHuang2006 demonstrated the fundamental result that $\ell_1$-penalized least squares estimators achieve the rate $\sqrt{s/n} \sqrt{\log p}$, which is very close to the oracle rate $\sqrt{s/n}$ achievable when the true model is known. vdGeer demonstrated a similar fundamental result on the excess forecasting error loss under both quadratic and non-quadratic loss functions. Thus the estimator can be consistent and can have excellent forecasting performance even under very rapid, nearly exponential, growth of the total number of regressors $p$. See FanLv2006,BickelRitovTsybakov2009,BuneaTsybakovWegkamp2006,BuneaTsybakovWegkamp2007,BuneaTsybakovWegkamp2007b,LouniciPontilTsybakovvandeGeer2009,RosenbaumTsybakov2008 for many other interesting developments and a detailed review of the existing literature.
Our paper's contribution is to develop a set of results on model selection and rates of convergence for quantile regression within the HDSM framework. Since ordinary quantile regression is inconsistent in HDSMs, we consider quantile regression penalized by the $\ell_1$-norm of parameter coefficients, denoted $\ell_1$-QR. First, we show that under general conditions $\ell_1$-QR estimates of regression coefficients and regression functions are consistent at the near-oracle rate $\sqrt{s/n} \sqrt{\log (p\vee n)}$, uniformly in a compact interval $\mathcal{U} \subset (0,1)$ of quantile indices.\footnote{Under $s \to \infty$, the oracle rate, uniformly over a proper compact interval $\mathcal{U}$, is $\sqrt{(s/n) \log n}$, cf. BC-qrID; the oracle rate for a single quantile index is $\sqrt{s/n}$, cf. HS00.} (This result is different from and hence complementary to vdGeer's fundamental results on the rates for excess forecasting error loss.) Second, in order to make $\ell_1$-QR practical, we propose a partly pivotal, data-driven choice of the penalty level, and show that this choice leads to the same sharp convergence rate. Third, we show that $\ell_1$-QR correctly selects the true model as a valid submodel when the non-zero coefficients of the true model are well separated from zero. Fourth, we also propose and analyze the post-penalized estimator (post-$\ell_1$-QR), which applies ordinary, unpenalized quantile regression to the model selected by the penalized estimator, and thus aims at reducing the regularization bias of the penalized estimator. We show that under similar conditions post-$\ell_1$-QR can perform as well as $\ell_1$-QR in terms of the rate of convergence, uniformly over $\mathcal{U}$, even if the $\ell_1$-QR-based model selection misses some components of the true models. This occurs because $\ell_1$-QR-based model selection only misses those components that have relatively small coefficients. Moreover, post-$\ell_1$-QR can perform better than $\ell_1$-QR if the $\ell_1$-QR-based model selection correctly includes all components of the true model as a subset. (Obviously, post-$\ell_1$-QR can perform as well as the oracle if the $\ell_1$-QR perfectly selects the true model, which is, however, unrealistic for many designs of interest.) Fifth, we illustrate the use of $\ell_1$-QR and post-$\ell_1$-QR with a Monte Carlo experiment and an international economic growth example. To the best of our knowledge, all of the above results are new and contribute to the literature on HDSMs. Our results on post-penalized estimators and some proof techniques could also be of interest in other problems. We provide further technical comparisons to the literature in Section (ref).
In what follows, we implicitly index all parameter values by the sample size $n$, but we omit the index whenever this does not cause confusion. We use the empirical process notation as defined in vdV-W. In particular, given a random sample $Z_1,...,Z_n$, let $\mathbb{G}_n (f) = \mathbb{G}_n (f(Z_i)) := n^{-1/2} \sum_{i=1}^n (f(Z_i) - \mathrm{E}[f(Z_i)])$ and $\mathbb{E}_n f = \mathbb{E}_n f(Z_i) := \sum_{i=1}^n f(Z_i)/n$. We use the notation $a \lesssim b$ to denote $a = O(b)$, that is, $a \leq c b$ for some constant $c>0$ that does not depend on $n$; and $a\lesssim_P b$ to denote $a=O_P(b)$. We also use the notation $a \vee b = \max\{ a, b\}$ and $a \wedge b = \min\{ a , b \}$. We denote the $\ell_2$-norm by $\|\cdot\|$, $\ell_1$-norm by $\|\cdot\|_1$, $\ell_{\infty}$-norm by $\|\cdot\|_\infty$, and the $\ell_0$-“norm" by $\|\cdot\|_0$ (i.e., the number of non-zero components). We denote by $\|\beta\|_{1,n} = \sum_{j=1}^p \widehat \sigma_j |\beta_j|$ the $\ell_1$-norm weighted by $\widehat\sigma_j$'s. Finally, given a vector $\delta \in {\rm I\kern-0.18em R}^p$, and a set of indices $T \subset \{1,\ldots,p\}$, we denote by $\delta_T$ the vector in which $\delta_{Tj} = \delta_j$ if $j\in T$, $\delta_{Tj}=0$ if $j \notin T$.
In this section we formulate the setting and the estimator, and state primitive regularity conditions. We also provide an overview of the main results.
The setting of interest corresponds to a parametric quantile regression model, where the dimension $p$ of the underlying model increases with the sample size $n$. Namely, we consider a response variable $y$ and $p$-dimensional covariates $x$ such that the $u$-th conditional quantile function of $y$ given $x$ is given by
where $\mathcal{U}\subset(0,1)$ is a compact set of quantile indices. Recall that the $u$-th conditional quantile $F^{-1}_{y_i|x_i}(u|x)$ is the inverse of the conditional distribution function $F_{y_i|x_i}(y|x_i)$ of $y_i$ given $x_i$. We consider the case where the dimension $p$ of the model is large, possibly much larger than the available sample size $n$, but the true model $\beta(u)$ has a sparse support $$T_u = {\rm support}(\beta(u)) = \{ j \in \{1,\ldots,p\} \ : \ |\beta_j(u)|>0 \} $$ having only $s_u \leq s \leq n/\log(n \vee p)$ non-zero components for all $u\in\mathcal{U}$.
The population coefficient $\beta(u)$ is known to minimize the criterion function
where $\rho_{u} (t) = (u - 1\{t\leq 0\})t$ is the asymmetric absolute deviation function KB78. Given a random sample $(y_1,x_1),\ldots,(y_n,x_n)$, the quantile regression estimator of $\beta(u)$ is defined as a minimizer of the empirical analog of ((ref)):
In high-dimensional settings, particularly when $p \geq n$, ordinary quantile regression is generally inconsistent, which motivates the use of penalization in order to remove all, or at least nearly all, regressors whose population coefficients are zero, thereby possibly restoring consistency. A penalization that has proven quite useful in least squares settings is the $\ell_1$-penalty leading to the Lasso estimator T1996.
The $\ell_1$-penalized quantile regression estimator $\widehat \beta (u)$ is a solution to the following optimization problem:
where $\widehat\sigma_j^2 = {\mathbb{E}_n}[x_{ij}^2]$. The criterion function in ((ref)) is the sum of the criterion function ((ref)) and a penalty function given by a scaled $\ell_1$-norm of the parameter vector. The overall penalty level $\lambda\sqrt{u(1-u)}$ depends on each quantile index $u$, while $\lambda$ will depend on the set $\mathcal{U}$ of quantile indices of interest. The $\ell_1$-penalized quantile regression has been considered in KnightFu2000 under small (fixed) $p$ asymptotics. It is important to note that the penalized quantile regression problem ((ref)) is equivalent to a linear programming problem (see Appendix (ref)) with a dual version that is useful for analyzing the sparsity of the solution. When the solution is not unique, we define $\widehat \beta(u)$ as any optimal basic feasible solution (see, e.g., BertsimasTsitsiklis). Therefore, the problem ((ref)) can be solved in polynomial time, avoiding the computational curse of dimensionality. Our goal is to derive the rate of convergence and model selection properties of this estimator.
The post-penalized estimator (post-$\ell_1$-QR) applies ordinary quantile regression to the model $\widehat T_u$ selected by the $\ell_1$-penalized quantile regression. Specifically, set $$\widehat T_u = {\rm support}( \widehat \beta(u) ) = \{ j \in \{1,\ldots,p\} \ : \ |\widehat\beta_j(u)| > 0\},$$ and define the post-penalized estimator $\widetilde \beta(u)$ as
which removes from further estimation the regressors that were not selected. If the model selection works perfectly -- that is, $\widehat T_u = T_u$ -- then this estimator is simply the oracle estimator, whose properties are well known. However, perfect model selection might be unlikely for many designs of interest. Rather, we are interested in the more realistic scenario where the first-step estimator $\widehat \beta(u)$ fails to select some components of $ \beta(u)$. Our goal is to derive the rate of convergence for the post-penalized estimator and show it can perform well under this scenario.
In order to describe our choice of the penalty level $\lambda$, we introduce the random variable
where $u_1,\ldots, u_n$ are i.i.d. uniform $(0,1)$ random variables, independently distributed from the regressors, $x_1,\ldots, x_n$. The random variable $\Lambda$ has a known, that is, pivotal, distribution conditional on $X= [x_1,\ldots, x_n]'$. We then set
where $\Lambda(1-\alpha|X) := (1-\alpha)$-quantile of $\Lambda$ conditional on $X$, and the constant $c>1$ depends on the design.\footnote{$c$ depends only on the constant $c_0$ appearing in condition D.4; when $c_0 \geq 9$, it suffices to set $c=2$.} Thus the penalty level depends on the pivotal quantity $\Lambda(1-\alpha|X)$ and the design. Under assumptions D.1-D.4 we can set $c=2$, similar to BickelRitovTsybakov2009's choice for least squares. Furthermore, we recommend computing $\Lambda(1-\alpha|X)$ using simulation of $\Lambda$.\footnote{We also provide analytical bounds on $\Lambda(1-\alpha|X)$ of the form $ C(\alpha,\mathcal{U}) \sqrt{n\log p}$ for some numeric constant $C(\alpha,\mathcal{U})$. We recommend simulation because it accounts for correlation among the columns of $X$ in the sample.} Our concrete recommendation for practice is to set $1-\alpha =0.9$.
The parameter $1-\alpha$ is the confidence level in the sense that, as in BickelRitovTsybakov2009, our (non-asymptotic) bounds on the estimation error will contract at the optimal rate with this probability. We refer the reader to Koenker Koenker2010 for an implementation of our choice of penalty level and practical suggestions concerning the confidence level. In particular, both here and in Koenker Koenker2010, the confidence level $1-\alpha =0.9$ gave good performance results in terms of balancing regularization bias with estimation variance. Cross-validation may also be used to choose the confidence level $1-\alpha$. Finally, we should note that, as in BickelRitovTsybakov2009, our theoretical bounds allow for any choice of $1-\alpha$ and are stated as a function of $1-\alpha$.
The formal rationale behind the choice ((ref)) for the penalty level $\lambda$ is that this choice leads precisely to the optimal rates of convergence for $\ell_1$-QR. (The same or slightly higher choice of $\lambda$ also guarantees good performance of post-$\ell_1$-QR.) Our general strategy for choosing $\lambda$ follows BickelRitovTsybakov2009, who recommend selecting $\lambda$ so that it dominates a relevant measure of noise in the sample criterion function, specifically the supremum norm of a suitably rescaled gradient of the sample criterion function evaluated at the true parameter value. In our case this general strategy leads precisely to the choice ((ref)). Indeed, a (sub)gradient $\widehat S_u(\beta(u)) = {\mathbb{E}_n}[( u - 1\{y_i \leq x_i'\beta(u)\})x_i] \in \partial \widehat Q_u (\beta(u))$ of the quantile regression objective function evaluated at the truth has a pivotal representation, namely $\widehat S_u(\beta(u))= {\mathbb{E}_n}[( u - 1\{u_i \leq u\})x_i]$ for $u_1,\ldots,u_n$ i.i.d. uniform $(0,1)$ conditional on $X$, and so we can represent $\Lambda$ as in ((ref)), and, thus, choose $\lambda$ as in ((ref)).
We consider the following conditions on a sequence of models indexed by $n$ with parameter dimension $p=p_n \to \infty$. In these conditions all constants can depend on $n$, but we omit the explicit indexing by $n$ to ease exposition.
D.1. Sampling and Smoothness. Data $(y_i,x_i')', i=1,\ldots,n,$ are an i.i.d. sequence of real $(1+p)$-vectors, with the conditional $u$-quantile function given by ((ref)) for each $u \in \mathcal{U}$, with the first component of $x_i$ equal to one, and $n \wedge p \geq 3$. For each value $x$ in the support of $x_i$, the conditional density $f_{y_i|x_i}(y|x)$ is continuously differentiable in $y$ at each $y \in \Bbb{R}$, and $f_{y_i|x_i}(y|x)$ and $\frac{\partial}{\partial y} f_{y_i|x_i}(y|x)$ are bounded in absolute value by constants $\bar{f}$ and $\bar{f'}$, uniformly in $y\in \mathbb{R}$ and $x$ in the support $x_i$. Moreover, the conditional density of $y_i$ evaluated at the conditional quantile $x_i'\beta(u)$ is bounded away from zero uniformly in $\mathcal{U}$, that is $f_{y_i|x_i}(x'\beta(u)|x)> \underline{f} >0$ uniformly in $u \in \mathcal{U}$ and $x$ in the support of $x_i$.
Condition D.1 imposes only mild smoothness assumptions on the conditional density of the response variable given regressors, and does not impose any normality or homoscedasticity assumptions\comment{ commonly made in the literature on HDSMs}. The assumption that the conditional density is bounded below at the conditional quantile is standard, but we can replace it by the slightly more general condition $\inf_{u \in \mathcal{U}} \inf_{\delta \neq 0} (\delta' J_u \delta)/(\delta'\mathrm{E}[x_i x_i'] \delta) \geq \underline{f}>0,$ on the Jacobian matrices $$J_u = \mathrm{E}[f_{y_i|x_i}(x_i'\beta(u)|x_i)x_ix_i'] \ \ \mbox{for all} \ \ u \in \mathcal{U},$$ throughout the paper.
D.2. Sparsity and Smoothness of $u \mapsto \beta(u)$. Let $\mathcal{U}$ be a compact subset of $(0,1)$. The coefficients $\beta(u)$ in ((ref)) are sparse and smooth with respect to $u\in \mathcal{U}$: $$ \sup_{u \in \mathcal{U}} \| \beta(u) \|_0 \leq s \ \ \ \ \mbox{and} \ \ \ \ \|\beta(u) - \beta(u')\| \leq L|u-u'|, \ \ \ \mbox{for all} \ \ u, u' \in \mathcal{U}$$ where $ s \geq 1$, and $\log L \leq C_L \log (p\vee n)$ for some constant $C_L$.
Condition D.2 imposes sparsity and smoothness on the behavior of the quantile regression coefficients $\beta(u)$ as we vary the quantile index $u$.
D.3. Well-behaved Covariates. Covariates are normalized such that $\sigma_j^2 = \mathrm{E}[x_{ij}^2] = 1$ for all $j = 1,\ldots,p$, and $\widehat \sigma_j^2 = {\mathbb{E}_n}[x^2_{ij}]$ obeys $P(\max_{1\leq j \leq p} \left| \widehat \sigma_{j} - 1 \right| \leq 1/2) \geq 1-\gamma \to 1$ as $n \to \infty$.
Condition D.3 requires that $\widehat \sigma_j$ does not deviate too much from $\sigma_j$ and normalizes $\sigma_j^2 = 1$.
In order to state the next assumption, for some $c_0 \geq 0$ and each $u \in \mathcal{U}$, define $$ A_u:=\{ \delta \in {\rm I\kern-0.18em R}^p : \ \|\delta_{T_u^c}\|_1 \leq c_0 \|\delta_{T_u}\|_1, \|\delta_{T_u^c}\|_0\leq n \},$$ which will be referred to as the restricted set. Define $\overline T_u(\delta,m) \subset\{1,...,p\}\setminus T_u$ as the support of the $m$ largest in absolute value components of the vector $\delta$ outside of $T_u={\rm support}(\beta(u))$, where $\overline T_u(\delta,m)$ is the empty set if $ m =0$.
D.4. Restricted Identifiability and Nonlinearity. For some constants $m \geq 0$ and $c_0 \geq 9$, the matrix $\mathrm{E}[x_i x_i']$ satisfies
and $\log (\underline{f} \kappa_0^2) \leq C_f \log (n \vee p)$ for some constant $C_f$. Moreover,
}
The restricted eigenvalue (RE) condition is analogous to the condition in BickelRitovTsybakov2009 and CandesTao2007; see BickelRitovTsybakov2009 and CandesTao2007 for different sufficient primitive conditions that yield bounds on $\kappa_m$. \comment{prove that the RE condition on the design matrix $\mathrm{E}[x_ix_i']$ is quite general, and is weaker than nearly all other design conditions used in the least squares literature;} Also, since $\kappa_m$ is non-increasing in $m$, $RE(c_0,m)$ for any $m>0$ implies $RE(c_0,0)$. The restricted non-linear impact (RNI) coefficient $q$ appearing in D.4 is a new concept, which controls the quality of minoration of the quantile regression objective function by a quadratic function over the restricted set.
Finally, we state another condition needed to derive results on the post-model selected estimator. In order to state the condition, define the sparse set $ \widetilde A_u(\widetilde m) = \{ \delta \in {\rm I\kern-0.18em R}^p : \|\delta_{T_u^c}\|_0 \leq \widetilde m \}$ for $\widetilde m \geq 0$ and $u\in\mathcal{U}$.
D.5. Sparse Identifiability and Nonlinearity. The matrix $\mathrm{E}[x_i x_i']$ satisfies for some $\widetilde m \geq 0$:
and
We invoke the sparse eigenvalue (SE) condition in order to analyze the post-penalized estimator ((ref)). This assumption is similar to the conditions used in MY2007 and ZhangHuang2006 to analyze Lasso. Our form of the SE condition is neither less nor more general than the RE condition. The SNI coefficient $\widetilde q_{\widetilde m}$ controls the quality of minoration of the quantile regression objective function by a quadratic function over sparse neighborhoods of the true parameter.
In order to highlight the nature and usefulness of conditions D.1-D.5 it is instructive to state some simple sufficient conditions (note that D.1-D.5 allow for much more general conditions). We relegate the proofs of this section to the Supplementary Material Appendix (ref) for brevity.
{\sc Design 1: Location Model with Correlated Normal Design}. Let us consider estimating a standard location model $$ y = x'\beta^o + \varepsilon, $$ where $\varepsilon \sim N(0,\sigma^2)$, $\sigma>0$ is fixed, $x=( 1, z' )'$, with $z \sim N(0, \Sigma),$ where $\Sigma$ has ones in the diagonal, a minimum eigenvalue bounded away from zero by a constant $\kappa^2>0$, and a maximum eigenvalue bounded from above, uniformly in $n$.
Note that the normality of errors can be easily relaxed by allowing for the disturbance $\varepsilon$ to have a smooth density that obeys the conditions stated in D.1. The conditions on the population design matrix can also be replaced by more general primitive conditions specified in Remark 2.1.
\comment{
}
{\sc Design 2: Location-scale model with bounded regressors}. Let us consider estimating a standard location-scale model $$ y = x'\beta^o + x'\eta \cdot \varepsilon, $$ where $\varepsilon \sim F$ independent of $x$, with a continuously differentiable probability density function $f$. We assume that the population design matrix $\mathrm{E}[xx']$ has ones in the diagonal and has eigenvalues uniformly bounded away from zero and from above, $x_1 =1$, $\max_{1\leq j\leq p} |x_{j}| \leq K_B$. Moreover, the vector $\eta$ is such that $0 < \upsilon \leq x'\eta \leq \Upsilon < \infty$ for all values of $x$.
\comment{
}
\comment{
}
Here we discuss our results under the simple setup of Design 1 and under $1/p \leq \alpha \to 0$ and $\gamma \to 0$. These simple assumptions allow us to straightforwardly compare our rate results to those obtained in the literature. We state our more general non-asymptotic results under general conditions in the subsequent sections. Our first main rate result is that $\ell_1$-QR, with our choice ((ref)) of parameter $\lambda$, satisfies
provided that the upper bound on the number of non-zero components $s$ satisfies
Note that $\kappa_0$, $\kappa_s$, $\underline{f}$, and $q$ are bounded away from zero in this example. Therefore, the rate of convergence is $\sqrt{s/n} \cdot \sqrt{\log (n \vee p)}$ uniformly in the set of quantile indices $u \in \mathcal{U}$, which is very close to the oracle rate when $p$ grows polynomially in $n$. Further, we note that our resulting restriction ((ref)) on the dimension $s$ of the true models is very weak; when $p$ is polynomial in $n$, $s$ can be of almost the same order as $n$, namely $s = o( n / \log n ) $.
Our second main result is that the dimension $\|\widehat \beta(u)\|_0 $ of the model selected by the $\ell_1$-penalized estimator is of the same stochastic order as the dimension $s$ of the true models, namely
Further, if the parameter values of the minimal true model are well separated from zero, then with a high probability the model selected by the $\ell_1$-penalized estimator correctly nests the true minimal model:
Moreover, we provide conditions under which a hard-thresholded version of the estimator selects the correct support.
Our third main result is that the post-penalized estimator, which applies ordinary quantile regression to the selected model, obeys
where $\widehat m= \sup_{u\in\mathcal{U}}\|\widehat\beta_{T_u^c}(u)\|_0$ is the maximum number of wrong components selected for any quantile index $u\in\mathcal{U}$, provided that the bound on the number of non-zero components $s$ obeys the growth condition ((ref)) and
(Note that when $\mathcal{U}$ is a singleton, the $s \log n$ factor in ((ref)) becomes $s$.)
We see from ((ref)) that post-$\ell_1$-QR can perform well in terms of the rate of convergence even if the selected model $\widehat T_u$ fails to contain the true model $T_u$. Indeed, since in this design $\widehat m \lesssim_P s$, post-$\ell_1$-QR has the rate of convergence $\sqrt{s/n} \cdot \sqrt{\log (n \vee p)}$, which is the same as the rate of convergence of $\ell_1$-QR. The intuition for this result is that the $\ell_1$-QR based model selection can only miss covariates with relatively small coefficients, which then permits post-$\ell_1$-QR to perform as well or even better due to reductions in bias, as confirmed by our computational experiments.
We also see from ((ref)) that post-$\ell_1$-QR can perform better than $\ell_1$-QR in terms of the rate of convergence if the number of wrong components selected obeys $\widehat m = o_P(s)$ and the selected model contains the true model, $\{T_u \subseteq \widehat T_u\}$ with probability converging to one. In this case post-$\ell_1$-QR has the rate of convergence $\sqrt{ (o_P(s)/n) \log (n \vee p) + (s/n) \log n },$ which is faster than the rate of convergence of $\ell_1$-QR. In the extreme case of perfect model selection, that is, when $\widehat m = 0$, the rate of post-$\ell_1$-QR becomes $\sqrt{(s/n) \log n}$ uniformly in $\mathcal{U}$. (When $\mathcal{U}$ is a singleton, the $\log n$ factor drops out.) Note that inclusion $\{T_u \subseteq \widehat T_u\}$ necessarily happens when the coefficients of the true models are well separated from zero, as we stated above. Note also that the condition $\widehat m = o(s)$ or even $\widehat m = 0$ could occur under additional conditions on the regressors (such as the mutual coherence conditions that restrict the maximal pairwise correlation of regressors). Finally, we note that our second restriction ((ref)) on the dimension $s$ of the true models is very weak in this design; when $p$ is polynomial in $n$, $s$ can be of almost the same order as $n$, namely $s = o( n / \log n ) $.
To the best of our knowledge, all of the results presented above are new, both for the single $\ell_1$-penalized quantile regression problem as well as for the infinite collection of $\ell_1$-penalized quantile regression problems. These results therefore contribute to the rate results obtained for $\ell_1$-penalized mean regression and related estimators in the fundamental papers of BickelRitovTsybakov2009,CandesTao2007,Koltchinskii2009,MY2007,vdGeer,ZhangHuang2006. The results on post-$\ell_1$ penalized quantile regression had no analogs in the literature on mean regression, apart from the rather exceptional case of perfect model selection, in which case the post-penalized estimator is simply the oracle. Building on the current work these results have been extended to mean regression in BC-PostLASSO. Our results on the sparsity of $\ell_1$-QR and model selection also contribute to the analogous results for mean regression MY2007. Also, our rate results for $\ell_1$-QR are different from, and hence complementary to, the fundamental results in vdGeer on the excess forecasting loss under possibly non-quadratic loss functions, which also specializes the results to density estimation, mean regression, and logistic regression. In principle we could apply theorems in vdGeer to the single quantile regression problem to derive the bounds on the excess loss $\mathrm{E}[\rho_u(y_i - x_i'\widehat \beta(u))] - \mathrm{E}[\rho_u(y_i - x_i'\beta(u))]$.\footnote{Of course, such a derivation would entail some difficult work, since we must verify some high-level assumptions made directly on the performance of the oracle and penalized estimators in population and others (cf. vdGeer's conditions I.1 and I.2, where I.2 assumes uniform in $x_i$ consistency of the penalized estimator in the population, and does not hold in our main examples, e.g., in Design 1 with normal regressors.)} However, these bounds would not imply our results ((ref)), ((ref)), ((ref)), ((ref)), and ((ref)), which characterize the rates of estimating coefficients $\beta(u)$ by $\ell_1$-QR and post-$\ell_1$-QR, sparsity and model selection properties, and the data-driven choice of the penalty level.
In this section we derive rates of convergence for $\ell_1$-QR and post-$\ell_1$-QR, sparsity bounds, and model selection results.
We start with a characterization of $\Lambda$ and its $(1-\alpha)$-quantile, $\Lambda(1-\alpha | X)$, which determines the magnitude of our suggested penalty level $\lambda$ via equation ((ref)).
In this section we establish the rate of convergence of $\ell_1$-QR. We start with the following preliminary result which shows that if the penalty level exceeds the specified threshold, for each $u \in \mathcal{U}$, the estimator $\widehat \beta(u) - \beta(u)$ will belong to the restricted set $ A_u:=\{ \delta \in {\rm I\kern-0.18em R}^p : \ \|\delta_{T_u^c}\|_1 \leq c_0 \|\delta_{T_u}\|_1, \|\delta_{T_u^c}\|_0\leq n \}$.
This result is inspired BickelRitovTsybakov2009's analogous result for least squares.
This second preliminary result derives identifiability relations over $A_u$. It shows that the coefficients $\underline{f}$, $\kappa_0$, and $\kappa_m$ control moduli of continuity between various norms over the restricted set $A_u$, and the RNI coefficient $q$ controls the quality of minoration of the objective function by a quadratic function over $A_u$.
Finally, the third preliminary result derives bounds on the empirical error over $A_u$:
In order to prove the lemma we use a combination of chaining arguments and exponential inequalities for contractions LedouxTalagrandBook. Our use of the contraction principle is inspired by its fundamentally innovative use in vdGeer; however, the use of the contraction principle alone is not sufficient in our case. Indeed, first we need to make some adjustments to obtain error bounds over the neighborhoods defined by the intrinsic norm $\|J^{1/2}_u \cdot \|$ instead of the $\|\cdot\|_1$ norm; and second, we need to use chaining over $u \in \mathcal{U}$ to get uniformity over $\mathcal{U}$.
Armed with Lemmas (ref)-(ref), we establish the first main result. The result depends on the constants $C_\Lambda$, $C_E$, $C_L$, and $C_f$ defined in Theorem (ref), Lemma (ref), D.2, and D.4.
This result derives the rate of convergence of the $\ell_1$-penalized quantile regression estimator in the intrinsic norm and other norms of interest uniformly in $u\in \mathcal{U}$ as well as uniformly in the penalty level $\lambda$ in the range specified by ((ref)), which includes our recommended choice of $\lambda_0$. We see that the rates of convergence for $\ell_1$-QR generally depend on the number of significant regressors $s$, the logarithm of the number of regressors $p$, the strength of identification summarized by $\kappa_0$, $\kappa_m$, $\underline{f}$, and $q$, and the quantile indices of interest $\mathcal{U}$ (as expected, extreme quantiles can slow down the rates of convergence). These rate results parallel the results of BickelRitovTsybakov2009 obtained for $\ell_1$-penalized mean regression. Indeed, the role of the parameter $\underline{f}$ is similar to the role of the standard deviation of the disturbance in mean regression. It is worth noting, however, that our results do not rely on normality and homoscedasticity assumptions, and our proofs have to address the non-quadratic nature of the objective function, with parameter $q$ controlling the quality of quadratization. This parameter $q$ enters the results only through the growth restriction ((ref)) on $s$. At this point we refer the reader to Section (ref) for a further discussion of this result in the context of the correlated normal design. Finally, we note that our proof combines the star-shaped geometry of the restricted set $A_u$ with classical convexity arguments; this insight may be of interest in other problems.
Next, we derive sparsity properties of the solution to $\ell_1$-penalized quantile regression. Fundamentally, sparsity is linked to the first order optimality conditions of ((ref)) and therefore to the (sub)gradient of the criterion function. In the case of least squares, the gradient is a smooth (linear) function of the parameters. In the case of quantile regression, the gradient is a highly non-smooth (piece-wise constant) function. To control the sparsity of $\widehat\beta(u)$ we rely on empirical process arguments to approximate gradients by smooth functions. In particular, we crucially exploit the fact that the entropy of all $m$-dimensional submodels of the $p$-dimensional model is of order $m \log p$, which depends on $p$ only logarithmically.
The statement of the results will depend on the maximal $k$-sparse eigenvalue of $\mathrm{E}\[x_ix_i'\right]$ and ${\mathbb{E}_n}\[x_ix_i'\right]$:
In order to establish our main sparsity result, we need two preliminary lemmas.
This lemma establishes an initial bound on the number of non-zero components $\widehat s$ as a function of $\lambda$ and $\phi(\widehat s)$. Restricting $\lambda \geq 2\sqrt{2}W_{\mathcal{U}}\sqrt{ n \log(n \vee p) \phi( n/\log(n\vee p))}$ makes the term $\phi\left(n/\log(n\vee p)\right)$ appear in subsequent bounds instead of the term $\phi(n)$, which in turn weakens some assumptions. Indeed, not only is the first term smaller than the second, but also there are designs of interest where the second term diverges while the first does not; for instance, in Design 1, if $p \geq 2n$, we have $\phi(n/\log (n\vee p)) \lesssim_P 1$ while $\phi(n) \gtrsim_P \sqrt{\log p}$ by the Supplementary Material Appendix (ref).
The following lemma establishes a bound on the sparsity as a function of the rate of convergence.
Finally, we combine these results to establish the main sparsity result. In what follows, we define $\bar \phi_{\varepsilon}$ as a constant such that $\phi(n/\log(n\vee p)) \leq \bar \phi_{\varepsilon}$ with probability $1-\varepsilon.$
The theorem states that by setting the penalty level $\lambda$ to be possibly higher than our initial recommended choice $\lambda_0$, we can control $\widehat s$, which will be crucial for good performance of the post-penalized estimator. As a corollary, we note that if (a) ${\mu} \lesssim 1$, (b) $1/(\underline{f} ^{1/2}\kappa_0) \lesssim 1$, and (c) $\bar \phi_{\varepsilon} \lesssim 1$ for each $\varepsilon>0$, then $\widehat s \lesssim s$ with a high probability, so the dimension of the selected model is about the same as the dimension of the true model. Conditions (a), (b), and (c) easily hold for the correlated normal design in Design 1. In particular, (c) follows from the concentration inequalities and from results in classical random matrix theory; see the Supplementary Material Appendix (ref) for proofs. Therefore the possibly higher $\lambda$ needed to achieve the stated sparsity bound does not slow down the rate of $\ell_1$-QR in this case. The growth condition ((ref)) on $s$ is also weak in this case.
Next we turn to the model selection properties of $\ell_1$-QR.
These results parallel analogous results in MY2007 for mean regression. The first result says that if non-zero coefficients are well separated from zero, then the support of $\ell_1$-QR includes the support of the true model. The inclusion of the true support in ((ref)) is in general one-sided; the support of the estimator can include some unnecessary components having true coefficients equal to zero. The second result states that if the further conditions are satisfied, additional hard thresholding can eliminate inclusions of such unnecessary components. The value of the hard threshold must explicitly depend on the unknown value $\min_{ j\in T_u }|\beta_j(u)|$, characterizing the separation of non-zero coefficients from zero. The additional conditions stated in this theorem are strong and perfect model selection appears quite unlikely in practice. Certainly it does not work in all real empirical examples we have explored. This motivates our analysis of the post-model-selected estimator under conditions that allow for imperfect model selection, including cases where we miss some non-zero components or have additional unnecessary components.
In this section we establish a bound on the rate of convergence of the post-penalized estimator. The proof relies crucially on the identifiability and control of the empirical error over the sparse sets $\widetilde A_u(\widetilde m):= \{ \delta \in {\rm I\kern-0.18em R}^p \ : \ \| \delta_{T_u^c} \|_0 \leq \widetilde m \}.$
In order to prove this lemma we exploit the crucial fact that the entropy of all $m$-dimensional submodels of the $p$-dimensional model is of order $m \log p$, which depends on $p$ only logarithmically. The following theorem establishes the properties of post-model-selection estimators.
This theorem describes the performance of post-$\ell_1$-QR. However, an inspection of the proof reveals that it can be applied to any post-model selection estimator. From Theorem (ref) we can conclude that in many interesting cases the rates of post-$\ell_1$-QR could be the same or faster than the rate of $\ell_1$-QR. Indeed, first consider the case where the model selection fails to contain the true model, i.e., $\sup_{u \in \mathcal{U}}1\{T_u \not\subseteq \widehat T_u\}=1$ with a non-negligible probability. If (a) $\widehat m \leq \widehat s \lesssim_P s$, (b) $\phi(\widehat m + s) \lesssim_P 1$, and (c) the constants $\underline{f}$ and $\widetilde \kappa_{\widehat m}^2$ are of the same order as $\underline{f}$ and $\kappa_0\kappa_{m}$, respectively, then the rate of convergence of post-$\ell_1$-QR is the same as the rate of convergence of $\ell_1$-QR. Recall that Theorem (ref) provides sufficient conditions needed to achieve (a), which hold in Design 1. Recall also that in Design 1, (b) holds by concentration of measure and classical results in random matrix theory, as shown in the Supplementary Material Appendix (ref), and (c) holds by the calculations presented in Section 2. This verifies our claim regarding the performance of post-$\ell_1$-QR in the overview, Section 2.4. The intuition for this result is that even though $\ell_1$-QR misses true components, it does not miss very important ones, allowing post-$\ell_1$-QR still to perform well. Second, consider the case where the model selection succeeds in containing the true model, i.e., $\sup_{u \in \mathcal{U}}1\{T_u \not\subseteq \widehat T_u\}=0$ with probability approaching one, and that the number of unnecessary components obeys $\widehat m = o_P(s)$. In this case the rate of convergence of post-$\ell_1$-QR can be faster than the rate of convergence of $\ell_1$-QR. In the extreme case of perfect model selection, when $\widehat m = 0$ with a high probability, post-$\ell_1$-QR becomes the oracle estimator with a high probability. We refer the reader to Section (ref) for further discussion, and note that this result could be of interest in other problems.
In order to access the finite sample practical performance of the proposed estimators, we conducted a Monte Carlo study and an application to international economic growth.
In order to assess the finite sample practical performance of the proposed estimators, we conducted a Monte Carlo study. We will compare the performance of the $\ell_1$-penalized, post-$\ell_1$-penalized, and the ideal oracle quantile regression estimators. Recall that the post-penalized estimator applies canonical quantile regression to the model selected by the penalized estimator. The oracle estimator applies canonical quantile regression to the true model. (Of course, such an estimator is not available outside Monte Carlo experiments.) We focus our attention on the model selection properties of the penalized estimator and biases and empirical risks of these estimators.
We begin by considering the following regression model: $$ y = x'\beta(0.5) + \varepsilon, \ \ \beta(0.5) =(1,1,1/2,1/3,1/4,1/5,0,\ldots,0)', $$ where as in Design 1, $x = (1,z')'$ consists of an intercept and covariates $z \sim N(0,\Sigma)$, and the errors $\varepsilon$ are independently and identically distributed $\varepsilon \sim N(0,\sigma^2)$. The dimension $p$ of covariates $x$ is $500$, and the dimension $s$ of the true model is $6$, and the sample size $n$ is $100$. We set the regularization parameter $\lambda$ equal to the $0.9$-quantile of the pivotal random variable $\Lambda$, following our proposal in Section (ref). The regressors are correlated with $\Sigma_{ij} = \rho^{|i-j|}$ and $\rho = 0.5$. We consider two levels of noise, namely $\sigma = 1$ and $\sigma=0.1$.
We summarize the model selection performance of the penalized estimator in Figures (ref) and (ref). In the left panels of the figures, we plot the frequencies of the dimensions of the selected model; in the right panels we plot the frequencies of selecting the correct regressors. From the left panels we see that the frequency of selecting a much larger model than the true model is very small in both designs. In the design with a larger noise, as the right panel of Figure (ref) shows, the penalized quantile regression never selects the entire true model correctly, always missing the regressors with small coefficients. However, it almost always includes the three regressors with the largest coefficients. (Notably, despite this partial failure of the model selection, post-penalized quantile regression still performs well, as we report below.) On the other hand, we see from the right panel of Figure (ref) that in the design with a lower noise level penalized quantile regression rarely misses any component of the true support. These results confirm the theoretical results of Theorem (ref), namely, that when the non-zero coefficients are well separated from zero, the penalized estimator should select a model that includes the true model as a subset. Moreover, these results also confirm the theoretical result of Theorem (ref), namely, that the dimension of the selected model should be of the same stochastic order as the dimension of the true model. In summary, the model selection performance of the penalized estimator agrees very well with our theoretical results.
We summarize results on estimation performance in Table (ref), which records for each estimator $\tilde \beta$ the norm of the bias $\|\mathrm{E}[\tilde \beta] - \beta_0\|$ and also the empirical risk $[\mathrm{E}[x_i' ( \tilde \beta - \beta_0 )]^2]^{1/2}$ for recovering the regression function. Penalized quantile regression has a substantial bias, as we would expect from the definition of the estimator which penalizes large deviations of coefficients from zero. We see that the post-penalized quantile regression drastically improves upon the penalized quantile regression, particularly in terms of reducing the bias, which results in a much lower overall empirical risk. Notably, despite that under the higher noise level the penalized estimator never recovers the true model correctly the post-penalized estimator still performs well. This is because the penalized estimator always manages to select the most important regressors. We also see that the empirical risk of the post-penalized estimator is within a factor of $\sqrt{\log p}$ of the empirical risk of the oracle estimator, as we would expect from our theoretical results. Under the lower noise level, the post-penalized estimator performs almost identically to the ideal oracle estimator. We would expect this since in this case the penalized estimator selects the model especially well, making the post-penalized estimator nearly the oracle. In summary, we find the estimation performance of the penalized and post-penalized estimators to be in close agreement with our theoretical results.
In this section we apply $\ell_1$-penalized quantile regression to an international economic growth example, using it primarily as a method for model selection. We use the Barro and Lee data consisting of a panel of 138 countries for the period of 1960 to 1985. We consider the national growth rates in gross domestic product (GDP) per capita as a dependent variable $y$ for the periods 1965-75 and 1975-85.\footnote{The growth rate in GDP over a period from $t_1$ to $t_2$ is commonly defined as $\log (GDP_{t_2}/GDP_{t_1})-1$.} In our analysis, we will consider a model with $p=60$ covariates, which allows for a total of $n=90$ complete observations. Our goal here is to select a subset of these covariates and briefly compare the resulting models to the standard models used in the empirical growth literature (Barro and Sala-i-Martin BarroLee1994, Koenker and Machado KoenkerMachado1999).
One of the central issues in the empirical growth literature is the estimation of the effect of an initial (lagged) level of GDP per capita on the growth rates of GDP per capita. In particular, a key prediction from the classical Solow-Swan-Ramsey growth model is the hypothesis of convergence, which states that poorer countries should typically grow faster and therefore should tend to catch up with the richer countries. Thus, such a hypothesis states that the effect of the initial level of GDP on the growth rate should be negative. As pointed out in Barro and Sala-i-Martin BarroSala1995, this hypothesis is rejected using a simple bivariate regression of growth rates on the initial level of GDP. (In our case, median regression yields a positive coefficient of $0.00045$.) In order to reconcile the data and the theory, the literature has focused on estimating the effect conditional on the pertinent characteristics of countries. Covariates that describe such characteristics can include variables measuring education and science policies, strength of market institutions, trade openness, savings rates and others BarroSala1995. The theory then predicts that for countries with similar other characteristics the effect of the initial level of GDP on the growth rate should be negative (BarroSala1995).
Given that the number of covariates we can condition on is comparable to the sample size, covariate selection becomes an important issue in this analysis (OneMillion, Two million). In particular, previous findings came under severe criticism for relying on ad hoc procedures for covariate selection. In fact, in some cases, all of the previous findings have been questioned (OneMillion). Since the number of covariates is high, there is no simple way to resolve the model selection problem using only classical tools. Indeed the number of possible lower-dimensional models is very large, although OneMillion and Two million attempt to search over several millions of these models. Here we use the Lasso selection device, specifically $\ell_1$-penalized median regression, to resolve this important issue.
Let us now turn to our empirical results. We performed covariate selection using $\ell_1$-penalized median regression, where we initially used our data-driven choice of penalization parameter $\lambda$. This initial choice led us to select no covariates, which is consistent with the situations in which the true coefficients are not well-separated from zero. We then proceeded to slowly decrease the penalization parameter in order to allow for some covariates to be selected. We present the model selection results in Table (ref). With the first relaxation of the choice of $\lambda$, we select the black market exchange rate premium (characterizing trade openness) and a measure of political instability. With a second relaxation of the choice of $\lambda$ we select an additional set of educational attainment variables, and several others reported in the table. With a third relaxation of $\lambda$ we include yet another set of variables also reported in the table. We refer the reader to BarroLee1994 and BarroSala1995 for a complete definition and discussion of each of these variables.
We then proceeded to apply ordinary median regression to the selected models and we also report the standard confidence intervals for these estimates. Table (ref) shows these results. We should note that the confidence intervals do not take into account that we have selected the models using the data. (In an ongoing companion work, we are working on devising procedures that will account for this.) We find that in all models with additional selected covariates, the median regression coefficients on the initial level of GDP is always negative and the standard confidence intervals do not include zero. Similar conclusions also hold for quantile regressions with quantile indices in the middle range. In summary, we believe that our empirical findings support the hypothesis of convergence from the classical Solow-Swan-Ramsey growth model. Of course, it would be good to find formal inferential methods to fully support this hypothesis. Finally, our findings also agree and thus support the previous findings reported in Barro and Sala-i-Martin BarroLee1994 and Koenker and Machado KoenkerMachado1999.
{
}