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.
85,429 characters · 15 sections · 61 citation commands
Oracle Estimation of a Change Point in High Dimensional Quantile Regression
\doublespacing
In this paper, we consider a high-dimensional quantile regression model where the sparsity structure (e.g., identities and effects of contributing regressors) may differ between two sub-populations, thereby allowing for a possible change point in the model. Let $Y \in \mathbb{R}$ be a response variable, $Q \in \mathbb{R}$ be a scalar random variable that determines a possible change point, and $X \in \mathbb{R}^{p}$ be a $p$-dimensional vector of covariates. Here, $Q$ can be a component of $X$, and $p$ is potentially much larger than the sample size $n$. Specifically, high-dimensional quantile regression with a change point is modelled as follows:
where $(\beta_0^T, \delta_0^T, \tau_0)$ is a vector of unknown parameters and the regression error $U$ satisfies $\mathbb{P}(U\leq 0|X,Q)=\gamma$ for some known $\gamma\in(0,1)$. Unlike mean regression, quantile regression analyzes the effects of active regressors on different parts of the conditional distribution of a response variable. Therefore, it allows the sparsity patterns to differ at different quantiles and also handles heterogeneity due to either heteroskedastic variance or other forms of non-location-scale covariate effects. By taking into account a possible change point in the model, we provide a more realistic picture of the sparsity patterns. For instance, when analyzing high-dimensional gene expression data, the identities of contributing genes may depend on the environmental or demographical variables (e.g., exposed temperature, age or weights).
Our paper is closely related to the literature on models with unknown change points (e.g., tong1990non, chan1993consistency, Hansen:1996, hansen2000sample, Pons2003, kosorok2007, Seijo:Sen:11a,Seijo:Sen:11b and Li:Ling:12 among many others). Recent papers on change points under high-dimensional setups include enikeeva2013high, chan2013group, Frick-et-al:14, cho2012multiple, Chan:et-al:2016, Callot:et-al:2016, and lee2012lasso among others; however, none of these papers consider a change point in high-dimensional quantile regression. The literature on high-dimensional quantile regression includes BC11, Bradic2011, Wang:Wu:Li:2012, Wang:2013, and FFB among others. All the aforementioned papers on quantile regression are under the homogeneous sparsity framework (equivalently, assuming that $\delta_0=0$ in (ref)). Ciuperca:13 considers penalized estimation of a quantile regression model with breaks, but the corresponding analysis is restricted to the case when $p$ is small.
In this paper, we consider estimating regression coefficients $\alpha_0 \equiv (\beta_0^T,\delta_0^T)^T$ as well as the threshold parameter $\tau_0$ and selecting the contributing regressors based on $\ell_1$-penalized estimators. One of the strengths of our proposed procedure is that it does not require to know or pretest whether $\delta_0=0$ or not, that is, whether the population's sparsity structure and covariate effects are invariant or not. In other words, we do not need to know whether the threshold $\tau_0$ is present in the model.
For a sparse vector $v\in\mathbb{R}^p$, we denote the active set of $v$ as $J(v) \equiv \{j: v_j\neq0\}$. One of the main contributions of this paper is that our proposed estimator of $\tau_0$ achieves an oracle property in the sense that its asymptotic distribution is the same as if the unknown active sets $J(\beta_{0})$ and $J(\delta_0)$ were known. Importantly, we establish this oracle property without assuming a perfect covariate selection, thereby avoiding the need for the minimum level condition on the signals of active covariates.
The proposed estimation method in this paper consists of three main steps: in the first step, we obtain the initial estimators of $\alpha_0$ and $\tau_0$, whose rates of convergence may be suboptimal; in the second step, we re-estimate $\tau_0$ to obtain an improved estimator of $\tau_0$ that converges at the rate of $O_P(n^{-1})$ and achieves the oracle property mentioned above; in the third step, using the second step estimator of $\tau_0$, we update the estimator of $\alpha_0$. In particular, we propose two alternative estimators of $\alpha_0$, depending on the purpose of estimation (prediction vs. variable selection).
The most closely related work is lee2012lasso. However, there are several important differences: first, lee2012lasso consider a high-dimensional mean regression model with a homoskedastic normal error and with deterministic covariates; second, their method consists of one-step least squares estimation with an $ \ell _1 $ penalty; third, they derive non-asymptotic oracle inequalities similar to those in Bickeletal but do not provide any distributional result on the estimator of the change point. Compared to lee2012lasso, dealing with high-dimensional quantile regression with an unknown change point calls for a new proof technique since the quantile loss function is different from the least squares objective function and is non-smooth. In addition, we allow for heteroskesdastic and non-normal regression errors and stochastic covariates. These changes coupled with the fact that the quantile regression objective function is non-convex with respect to the threshold parameter $\tau_0$ raise new challenges. It requires careful derivation and multiple estimation steps to establish the oracle property for the estimator of $\tau_0$ and also to obtain desirable properties of the estimator of $\alpha_0$. The technique developed in this paper is applicable to a general M-estimation framework with a change point, which may be of independent interest.
One particular application of (ref) comes from tipping in the racial segregation in social sciences card2008tipping. The empirical question addressed in card2008tipping is whether and the extent to which the neighborhood's white population decreases substantially when the minority share in the area exceeds a tipping point (or change point). In Section (ref), we use the US Census tract dataset constructed by card2008tipping and confirm that the tipping exists in the neighborhoods of Chicago.
The remainder of the paper is organized as follows. Section (ref) provides an informal description of our estimation methodology. In Section (ref), we derive the consistency of the estimators in terms of the excess risk. Further asymptotic properties of the proposed estimators are given in Sections (ref) and (ref). In Section (ref), we present the results of extensive Monte Carlo experiments. Section (ref) illustrates the usefulness of our method by applying it to tipping in the racial segregation. Section (ref) concludes and Appendix (ref) describes in detail regarding how to construct the confidence interval for $\tau_0$. In Appendix (ref), we provide a set of regularity assumptions to derive asymptotic properties of the proposed estimators in Sections (ref) and (ref). Online supplements are comprised of 6 appendices for all the proofs as well as additional theoretical and numerical results that are left out for the brevity of the paper.
Notation. Throughout the paper, we use $|v|_q$ for the $\ell_q$ norm for a vector $v$ with $q=0,1,2$. We use $|v|_\infty$ to denote the sup norm. For two sequences of positive real numbers $a_n$ and $b_n$, we write $a_n\ll b_n$ and equivalently $b_n\gg a_n$ if $a_n=o(b_n)$. If there exists a positive finite constant $c$ such that $a_n = c \cdot b_n$, then we write $a_n \propto b_n$. Let $\lambda_{\min}(A)$ denote the minimum eigenvalue of a matrix $A$. We use w.p.a.1 to mean “with probability approaching one.” We write $\theta_0 \equiv \beta_0+\delta_0$. For a $2p$ dimensional vector $\alpha$, let $\alpha_J$ and $\alpha_{J^c}$ denote its subvectors formed by indices in $J(\alpha_0)$ and $\{1,...,2p\}\setminus J(\alpha_0)$, respectively. Likewise, let $X_{J}(\tau)$ denote the subvector of $X(\tau)\equiv (X^{T},X^{T}1\{Q>\tau \})^{T}$ whose indices are in $J(\alpha_0)$. The true parameter vectors $\beta_0$, $\delta_0$ and $\theta_0$ (except $\tau_0$) are implicitly indexed by the sample size $n$, and we allow that the dimensions of $J(\beta_0)$, $J(\delta_0)$ and $J(\theta_0)$ can go to infinity as $n \rightarrow \infty$. For simplicity, we omit their dependence on $n$ in our notation. We also use the terms `change point' and `threshold' interchangeably throughout the paper.
In this section, we describe our estimation method. We take the check function approach of Koenker:Bassett:1978. Let $\rho(t_1,t_2) \equiv (t_1-t_2)(\gamma-1\{t_1-t_2\leq 0\})$ denote the loss function for quantile regression. Let $\mathcal{A}$ and $\mathcal{T}$ denote the parameter spaces for $\alpha_0 \equiv (\beta_0^T, \delta_0^T)^T$ and $\tau_0$, respectively. For each $\alpha \equiv (\beta,\delta)\in\mathcal{A}$ and $\tau\in\mathcal{T}$, we write $X^T\beta+X^T \delta1\{Q>\tau\}=X(\tau)^T\alpha$ with the shorthand notation that $X(\tau) \equiv (X^{T},X^{T}1\{Q>\tau \})^{T}$. We suppose that the vector of true parameters is defined as the minimizer of the expected loss:
By construction, $\tau_0$ is not unique when $\delta_0=0$. However, if $\delta_0 = 0$, then the model reduces to the linear quantile regression model in which $\beta_0$ is identifiable under the standard assumptions. In Appendix (ref), we provide sufficient conditions under which $\alpha_0$ and $\tau_0$ are identified when $\delta_0 \neq 0$.
Suppose we observe independent and identically distributed samples $\{Y_i, X_i, Q_i\}_{i\leq n}$. Let $X_{i}(\tau )$ and $X_{ij}\left( \tau \right) $ denote the $i$-th realization of $X(\tau )$ and $j$-th element of $ X_{i}\left( \tau \right) ,$ respectively, $i=1,\ldots,n$ and $j=1,\ldots,2p$, so that $X_{ij}(\tau) \equiv X_{ij}$ if $j\leq p$ and $X_{ij}(\tau) \equiv X_{i,j-p}1\{Q_i>\tau\}$ otherwise. Define $$ R_n(\alpha,\tau)\equiv \frac{1}{n} \sum_{i=1}^{n}\rho (Y_{i},X_{i}(\tau )^{T}\alpha )= \frac{1}{n}\sum_{i=1}^{n}\rho (Y_{i},X_{i}^T\beta+ X_i^T\delta 1\{Q_i>\tau\} ). $$ In addition, let $D_{j}(\tau ) \equiv \{ n^{-1} \sum_{i=1}^{n}X_{ij}(\tau )^{2} \}^{1/2}$, $j=1,\ldots,2p$.
We describe the main steps of our $\ell_1$-penalized estimation method. For some tuning parameter $\kappa_{n}$, define:
This step produces an initial estimator $(\breve\alpha,\breve\tau)$. The tuning parameter $\kappa_n$ is required to satisfy
Note that we take $\kappa_n$ that converges to zero at a rate slower than the standard $(\log p/n)^{1/2} $ rate in the literature. This modified rate of $\kappa_n$ is useful in our context to deal with an unknown $\tau_0$. A data-dependent method of choosing $\kappa _{n}$ is discussed in Section (ref).
The main purpose of the first step is to obtain an initial estimator of $\alpha_0$. The achieved convergence rates of this step might be suboptimal due to the uniform control of the score functions over the space $\mathcal{T}$ of the unknown $\tau_0$.
In the second step, we introduce our improved estimator of the change point $\tau_0$. It does not use a penalty term, while using the first step estimator of $\alpha_0$. Define:
where $\breve{\alpha}$ is the first step estimator of $\alpha_0$ in (ref). In Section (ref), we show that when $\tau_0$ is identifiable, $\widehat{\tau}$ is consistent for $\tau_0$ at a rate of $n^{-1}$. Furthermore, we obtain the limiting distribution of $n(\widehat{\tau} - \tau_0)$, and establish conditions under which its asymptotic distribution is the same as if the true $\alpha_0$ were known, without a perfect model selection on $\alpha_0$, nor assuming the minimum signal condition on the nonzero elements of $\alpha_0$.
In the third step, we update the Lasso estimator of $\alpha_0$ using a different value of the penalization tuning parameter and the second step estimator of $\tau_0$. In particular, we recommend two different estimators of $\alpha_0$\,: one for the prediction and the other for the variable selection, serving for different purposes of practitioners. For two different tuning parameters $\omega_n$ and $\mu_n$ whose rates will be specified later by (ref) and (ref), define:
where $\widehat{\tau}$ is the second step estimator of $\tau_0$ in (ref), and the “signal-adaptive" weight $w_j$ in (ref), motivated by the local linear approximation of the SCAD penalties Fan01, zouli, is calculated based on the Step 3a estimator $\widehat\alpha$ from ((ref)):
Here $a>1$ is some prescribed constant, and $a=3.7$ is often used in the literature. We take this as our choice of $a$.
Step 3 defines two estimators for $\alpha_0$. In this subsection we briefly explain their major differences and purposes. Step 3b is particularly useful when the variable selection consistency is the main objective, yet it often requires the minimum signal condition ($\min_{\alpha_{0j}\neq0}|\alpha_{0j}|$ is well separated from zero). In contrast, Step 3a does not require the minimum signal condition, and is recommended for prediction purposes. More specifically:
In this subsection, we provide details on how to choose tuning parameters in applications. Recall that our procedure involves three tuning parameters in the penalization: (1) $\kappa_n$ in Step 1 ought to dominate the score function uniformly over the range of $\tau$, and hence should be slightly larger than the others; (2) $\omega_n$ is used in Step 3a for the prediction, and (3) $\mu_n$ in Step 3b for the variable selection should be larger than $\omega_n$. Note that the tuning parameters in both Steps 3a and 3b are similar to those of the existing literature since the change point $\widehat\tau$ has been estimated.
We build on the data-dependent selection method in BC11. Define
where $U_i$ is simulated from the i.i.d.\ uniform distribution on the interval $[0,1]$; $\gamma$ is the quantile of interest (e.g. $\gamma = 0.5$ for median regression). Note that $\Lambda(\tau)$ is a stochastic process indexed by $\tau$. Let $\overline{\Lambda}_{1-\epsilon^\ast}$ be the $(1-\epsilon^\ast)$-quantile of $\sup_{\tau \in \mathcal{T}} \Lambda(\tau)$, where $\epsilon^\ast$ is a small positive constant that will be selected by a user. Then, we select the tuning parameter in Step 1 by $ \kappa_n = c_1 \cdot \overline{\Lambda}_{1-\epsilon^\ast}. $ Similarly, let $\Lambda_{1-\epsilon^\ast}(\widehat{\tau})$ be the $(1-\epsilon^\ast)$-quantile of $\Lambda(\widehat{\tau})$, where $\widehat{\tau}$ is chosen in Step 2. We select $\omega_n$ and $\mu_n$ in Step 3 by $ \omega_n = c_1 \cdot \Lambda_{1-\epsilon^\ast}(\widehat{\tau}) $ and $ \mu_n = c_2 \cdot \omega_n. $ It is also necessary to choose $\mathcal{T}$ in applications. In our Monte Carlo experiments in Section (ref), we take $\mathcal{T}$ to be the interval from the 15th percentile to the 85th percentile of the empirical distribution of the threshold variable $Q_i$. For example, Hansen:1996 employed the same range in his application to U.S. GNP dynamics.
Based on the suggestions of BC11 and some preliminary simulations, we choose to set $c_1=1.1$, $c_2=\log\log n$, and $\epsilon^\ast=0.1$. In addition, recall that we set $a=3.7$ when calculating the SCAD weights $w_j$ in Step 3b following the convention in the literature (e.g.\ Fan01 and loh2013regularized). In Step 1, we first solve the lasso problem for $\alpha$ given each grid point of $\tau\in \mathcal{T}$. Then, we choose $\breve{\tau}$ and the corresponding $\breve{\alpha}(\breve{\tau})$ that minimize the objective function. Step 2 can be solved simply by the grid search. Step 3 is a standard lasso quantile regression estimation given $\widehat{\tau}$, whose numerical implementation is well established. We use the rq() function of the R `quantreg' package with the method = "lasso" in each implementation of the standard lasso quantile regression estimation quantreg.
Given the loss function $\rho(t_1,t_2) \equiv (t_1-t_2)(\gamma-1\{t_1-t_2\leq 0\})$ for the quantile regression model, define the excess risk to be
By the definition of $(\alpha_0,\tau_0)$ in (ref), we have that $R(\alpha,\tau) \geq 0$ for any $\alpha \in \mathcal{A}$ and $\tau \in \mathcal{T}$. What we mean by the “risk consistency” here is that the excess risk converges in probability to zero for the proposed estimators. The other asymptotic properties of the proposed estimators will be presented in Sections (ref) and (ref).
In this section, we begin by stating regularity conditions that are needed to develop our first theoretical result. Recall that $X_{ij}$ denotes the $j$th element of $X_i$.
In addition to the random sampling assumption, condition (ref) imposes mild moment restrictions on $X$. Condition (ref) imposes a weak restriction that the probability that $Q \in (\tau_1, \tau_2]$ is bounded by a constant times $(\tau_2 - \tau_1)$. Condition (ref) assumes that the parameter space is compact and that the support of $Q$ is strictly larger than $\mathcal{T}$. These conditions are standard in the literature on change-point and threshold models (e.g., Seijo:Sen:11a,Seijo:Sen:11b). Condition (ref) also assumes that the conditional expectation of $\mathbb{E}[ X_{ij}^2 | Q = \cdot ]$ is bounded on $\mathcal{T}$ uniformly in $j$. Condition (ref) requires that each regressor be of the same magnitude uniformly over the threshold $\tau$. As the data-dependent weights $D_j(\tau)$ are the sample second moments of the regressors, it is not stringent to assume them to be bounded away from both zero and infinity. Condition (ref) puts some weak upper bound on $\mathbb{E}[ \left( X^{T}\delta _{0}\right) ^{2}|Q=\tau ]$ for all $\tau \in \mathcal{T}$ when $\delta_0 \neq 0$. A simple sufficient condition for condition (ref) is that the eigenvalues of $\mathbb{E}[ X_{J(\delta_0)}X_{J(\delta_0)}^T|Q=\tau ]$ are bounded uniformly in $\tau$, where $X_{J(\delta_0)}$ denotes the subvector of $X$ corresponding to the nonzero components of $\delta_0$.
Throughout the paper, we let $s \equiv |J(\alpha_0)|_0$, namely the cardinality of $J(\alpha_0)$. We allow that $s \rightarrow \infty$ as $n \rightarrow \infty$ and will give precise regularity conditions regarding its growth rates. The following theorem is concerned about the convergence of $ R(\breve{\alpha},\breve{\tau})$ with the first step estimator.
Note that Theorem (ref) holds regardless of the identifiability of $\tau_0$ (that is, whether $\delta_0 = 0$ or not). In addition, the rate $O_P(\kappa_ns)$ is achieved regardless of whether $\kappa_ns$ converges, and we have the risk consistency if $\kappa_n s \rightarrow 0$ as $n \rightarrow \infty$. The restriction on $s$ is slightly stronger than that of the standard result $s = o ( \sqrt{n/\log p} )$ in the literature for the M-estimation (see, e.g. geer and Chapter 6.6 of bulmann) since the objective function $\rho(Y, X(\tau)^T\alpha)$ is non-convex in $\tau$, due to the unknown change-point.
In Appendix (ref), we show that an improved rate of convergence, $O_P\left( \omega_{n}s\right)$, is possible for the excess risk by taking the second and third steps of estimation.
Sections (ref) and (ref) provide asymptotic properties of the proposed estimators. In Appendix (ref), we list a set of assumptions that are needed to derive these properties, in addition to Assumption (ref). We first establish the consistency of $\breve{\protect\tau}$ for $\tau_0$.
The following theorem presents the rates of convergence for the first step estimators of $\alpha_0$ and $\tau_0$. Recall that $\kappa_n$ is the first-step penalization tuning parameter that satisfies (ref).
In Theorem (ref), we have that $ R(\breve{\alpha },\breve{\tau})=O_P\left( \kappa_{n}s\right). $ The improved rate of convergence for $ R(\breve{\alpha },\breve{\tau})$ in Theorem (ref) is due to additional assumptions (in particular, compatibility conditions in Assumption (ref) among others). It is worth noting that $\breve\tau$ converges to $\tau_0$ faster than the standard parametric rate of $n^{-1/2}$, as long as $s^2(\log p)^6(\log n)^4 =o(n)$. The main reason for such super-consistency is that the objective function behaves locally linearly around $\tau_0$ with a kink at $\tau_0$, unlike in the regular estimation problem where the objective function behaves locally quadratically around the true parameter value. Moreover, the achieved convergence rate for $\breve\alpha$ is nearly minimax optimal, with an additional factor $(\log p)(\log n)$ compared to the rate of regular Lasso estimation (e.g., Bickeletal, Rasku09). This factor arises due to the unknown change-point $\tau_0.$ We will improve the rates of convergence for both $\tau_0$ and $\alpha_0$ further by taking the second and third steps of estimation.
Recall that the second-step estimator of $\tau_0$ is defined as
where $\breve{\alpha}$ is the first step estimator of $\alpha_0$ in (ref). Consider an oracle case for which $\alpha$ in $R_n(\alpha,\tau)$ is fixed at $\alpha_{0}$. Let $R_{n}^{\ast}\left( \tau\right) =R_{n} \left( \alpha_{0},\tau\right) $ and \[ \widetilde{\tau}=\operatorname*{argmin}_{\tau \in \mathcal{T}} R_{n}^{\ast}\left( \tau\right) . \]
We now give one of the main results of this paper.
The first conclusion of Theorem (ref) establishes that the second step estimator of $\tau_0$ is an oracle estimator in the sense that it is asymptotically equivalent to the infeasible, oracle estimator $\widetilde{\tau}$. As emphasized in the introduction, the oracle property is obtained without relying on the perfect model selection in the first step nor on the existence of the minimum signal condition on active covariates. The second conclusion of Theorem (ref) follows from combining well-known weak convergence results in the literature (see e.g. Pons2003, kosorok2007, Lee:Seo:08) with the argmax continuous mapping theorem by Seijo:Sen:11b.
We now consider the Step 3a estimator of $\alpha_0$ defined in (ref). Recall that $\omega_n$ is the Step 3a penalization tuning parameter that satisfies (ref).
Theorem (ref) shows that the estimator $\widehat\alpha$ defined in Step 3a achieves the optimal rate of convergence in terms of prediction and estimation. In other words, when $\omega_{n}$ is proportional to $\{ \log (p \vee n)/n \}^{1/2}$ in equation (ref) and $p$ is larger than $n$, it obtains the minimax rates as in e.g., Rasku09.
As we mentioned in Section (ref), the Step 3b estimator of $\alpha_0$ has the purpose of the variable selection. The nonzero components of $\widetilde\alpha$ are expected to identify contributing regressors. Partition $\widetilde{\alpha }=(\widetilde{\alpha }_{J},\widetilde{ \alpha }_{J^{c}})$ such that $\widetilde{\alpha }_{J}=(\widetilde{\alpha } _{j}:j\in J(\alpha_0))$ and $\widetilde{\alpha }_{J^{c}}= ( \widetilde{\alpha } _{j}:j\notin J(\alpha_0) )$. Note that $\widetilde\alpha_J$ consists of the estimators of $\beta_{0J}$ and $\delta_{0J}$, whereas $\widetilde\alpha_{J^c}$ consists of the estimators of all the zero components of $\beta_0$ and $\delta_0.$ Let $\alpha_{0J}^{(j)}$ denote the $j$-th element of $\alpha_{0J}$.
We now establish conditions under which the estimator $\widetilde\alpha$ defined in Step 3b has the change-point-oracle properties, meaning that it achieves the variable selection consistency and has the limiting distributions as though the identities of the important regressors and the location of the change point were known.
We see that ((ref)) provides a condition on the strength of the signal via $\min_{j\in J(\alpha_0)}|\alpha _{0J}^{(j)}|$, and the tuning parameter in Step 3b should satisfy $\omega_n \ll \mu_n$ and $s^2 \log s/n \ll \mu_n^2$. Hence the variable selection consistency demands a larger tuning parameter than in Step 3a.
To conduct statistical inference, we now discuss the asymptotic distribution of $\widetilde\alpha_J$. Define $\widehat{\alpha}^{\ast}_J\equiv\operatorname*{argmin}_{\alpha_J }R_{n}^{\ast}\left( \alpha_J,\tau_{0}\right)$. Note that the asymptotic distribution for $\widehat{\alpha}^{\ast}_J$ corresponds to an oracle case that we know $\tau _{0}$ as well as the true active set $J(\alpha_0)$ a priori. The limiting distribution of $\widetilde\alpha_J$ is the same as that of $\widehat{\alpha}^{\ast}_J$. Hence, we call this result the change-point-oracle property of the Step 3b estimator and the following theorem establishes this property.
Since the sparsity index ($s$) grows at a rate slower than the sample size ($n$), it is straightforward to establish the asymptotic normality of a linear transformation of $\widetilde\alpha_J,$ i.e., $\mathbf{L} \widetilde\alpha_J,$ where $\mathbf{L}:\mathbb{R}^{s}\rightarrow\mathbb{R}$ with $|\mathbf{L}|_2=1$, by combing the existing results on quantile regression with parameters of increasing dimension (see, e.g. He:Shao:00) with Theorem (ref).
In this section, we show that our estimators have desirable results even if there is no change point in the true model. The case of $\delta_0=0$ corresponds to the high-dimensional linear quantile regression model. Since $ X^T\beta_0+X^T\delta_01\{Q>\tau_0\}=X^T\beta_0, $ $\tau_0$ is non-identifiable, and there is no structural change on the coefficient. But a new analysis different from that of the standard high-dimensional model is still required because in practice we do not know whether $\delta_0=0$ or not. Thus, the proposed estimation method still estimates $\tau_0$ to account for possible structural changes. The following results show that in this case, the first step estimator of $\alpha_0$ will asymptotically behave as if $\delta_0=0$ were a priori known.
The results obtained in Theorem (ref) combined with those obtained in Theorem (ref) imply that the first step estimatior performs equally well in terms of rates of convergence for both the $\ell_1$ loss for $\breve\alpha$ and the excess risk regardless of the existence of the threshold effect. It is straightforward to obtain an improved rate result for the Step 3a estimator, equivalent to Theorem (ref) under Assumptions (ref)-(ref). We omit the details for brevity.
We now give a result that is similar to Theorem (ref) and Theorem (ref).
Theorem (ref) demonstrates that when there is in fact no change point, our estimator for $\delta_0$ is exactly zero with a high probability. Therefore, the estimator can also be used as a diagnostic tool to check whether there exists any change point. Results similar to Theorems (ref) can be established straightforwardly as well; however, their details are omitted for brevity.
In this section we provide the results of Monte Carlo experiments. The baseline model is based on the following data generating process: for $i=1,\ldots,n$,
where $U_{i}$ follows $N(0,0.5^2)$, and $Q_i$ follows the uniform distribution on the interval $[0,1]$. The $p$-dimensional covariate $X_i$ is composed of a constant and $Z_i$, i.e.\ $X:=(1,Z_i^{T})^{T}$, where $Z_i$ follows the multivariate normal distribution $N(0,\Sigma)$ with a covariance matrix $\Sigma_{ij} = ( 1 / 2 )^{\vert i-j\vert}$. Here, the variables $U_{i}, Q_i$ and $Z_i$ are independent of each other. Note that the conditional $\gamma$-quantile of $Y_i$ given $(X_i,Q_i)$ has the form:
where $\beta_{\gamma}=\beta_0+\xi_{10}\cdot Quant_{\gamma}(U)$ and $\delta_{\gamma}=\delta_0+\xi_{20} \cdot Quant_{\gamma}(U)$.
We consider three quantile regression models with $\gamma=0.25, 0.5$, and $0.75$. The $p$-dimensional parameters $\beta_0$, $\delta_0$, $\xi_{10}$, and $\xi_{20}$ are set to $\beta_0=(0,Quant_{0.75}(U)\approx0.34,0,\ldots,0)$, $\delta_0=(0,1,0,\ldots,0)$, $\xi_{10}=(0,1,0,\ldots,0)$, and $\xi_{20}=(0,0,0,\ldots,0)$, respectively. Because of the heteroskedasticity, the true parameter value $\beta_{\gamma}$ at each quantile is $\beta_{0.25}=(0,\ldots,0)$, $\beta_{0.5}=(0,0.34, \ldots,0)$, and $\beta_{0.75}=(0,0.68, \ldots,0)$. Note that nonzero coefficients are different between when $\gamma=0.25$ and when $\gamma=0.5$ or $\gamma=0.75$.
We set the change point parameter $\tau_0=0.5$ unless it is specified differently. The sample sizes are set to $n=200$ and $400$. The dimension of $X_i$ is set to $p=250$. Note that we have $500$ regressors in total. The change point $\tau$ is estimated over grid points of the sample observations $\{Q_i\}$, where the range is limited to those between the $0.15$-quantile and the $0.85$-quantile. We conduct $1{,}000$ replications of each design.
We compare estimation results of each step. To assess the performance of our estimators, we also compare the results with two “oracle estimators". Specifically, Oracle 1 knows the true active set $J(\alpha_{\gamma})$ and the change point parameter $\tau_0$, and Oracle 2 knows only $J(\alpha_{\gamma})$. The threshold parameter $\tau_0$ is re-estimated in Steps 3a and 3b using updated estimates of $\alpha_{\gamma}$.
Tables (ref)--(ref) summarize the simulation results. We abuse notation slightly and denote all estimators by $(\widehat{\alpha},\widehat{\tau})$. They would be understood as $(\breve{\alpha},\breve{\tau})$ in Step 1, $\widehat{\tau}$ in Step 2, and so on. We report the excess risk, the average number of parameters selected, $\mathbb{E}[J(\widehat{\alpha})]$, and the sum of the mean squared error of $\widehat{\alpha}$ ($\widehat{\alpha}_{J_0}$ / $\widehat{\alpha}_{J^c_0}$). For each sample, the excess risk is calculated by the simulation, $S^{-1}\sum_{s=1}^S \left[\rho(Y_s,X_s^T(\widehat{\tau})\widehat{\alpha})-\rho(Y_s,X_s^T({\tau}_0){\alpha}_0)\right]$, where $S=10{,}000$ is the number of simulations; then we report the average value of 1,000 replications. Similarly, we also calculate prediction errors by the simulation, $\left(S^{-1}\sum_{s=1}^S \left(X^T_s(\widehat{\tau})\widehat{\alpha} - X^T_s(\tau_0)\alpha_{\gamma}\right)^2\right)^{1/2}$, and report the average value.
We also report the root-mean-squared error (RMSE) and the coverage probability of the 95% confidence interval of $\widehat{\tau}$ (C.\ Prob.\ of $\widehat{\tau}$). The confidence intervals for $\tau_0$ are calculated by simulating the two-sided compound Poisson process in Theorem (ref) by adopting the approach proposed by Li:Ling:12. The details are provided in Section (ref). Li:Ling:12 showed that it is valid to simulate the compound poisson process by simulating the poisson process and the compounding factors from empirical distributions separately in the context of least squares estimation. We build on their suggestion and modify their procedure to quantile regression. We did not prove a formal justification for our procedure in this paper; however, it seems working well in simulations. It is an interesting topic for future research.
Note that the root-mean-squared error of $\widehat{\tau}$ and the coverage probability of the confidence interval at the rows of Step 3a and Step 3b in the tables are estimation results of updated $\widehat{\tau}$: we re-estimate $\tau$ as in Step 2 using $(\widehat{U}_i,\widehat{\alpha})$ and $(\widetilde{U}_i,\widetilde{\alpha})$ from Step 3a and Step 3b instead of $(\breve{U}_i,\breve{\alpha})$. Finally, we also report the oracle proportion (Oracle Prop.), namely the ratio of the correct model selection out of 1,000 replications.
Overall, the simulation results confirm the asymptotic theory developed in the previous sections. First, these results show the advantage of quantile regression models over the existing mean regression models with a change point, e.g.\ lee2012lasso. The proposed estimator (Step 3b) selects different nonzero coefficients at different quantile levels. The estimator in lee2012lasso cannot detect these heterogeneous models. In general the proposed estimators show better performance for heteroskedastic designs and for the fat-tail error distributions as will be discussed in detail below. Second, when we look at the finite sample performance of the proposed estimators in Step 3, their prediction errors are within a reasonable bound from those of Oracles 1 and 2. Recall that we estimate models with 250 times or 500 times more regressors in each design. Third, the root-mean-squared error of $\widehat{\tau}$ decreases quickly and confirms the super-consistency result of $\widehat{\tau}$. Fourth, the coverage probabilities of the confidence interval are close to 95%, especially when $n=400$. Thus, we recommend practitioners to use $\widehat{\tau}$ in Step 2 or the re-estimated version of it based on the estimates from Step 3a or Step 3b. Finally, the oracle proportion of Step 3b is quite satisfactory and confirms our results in model selection consistency.
Table (ref) compares the performance of the proposed estimator with that of the mean regression method in lee2012lasso. For the purpose of direct comparison between mean and median regression models, the tuning parameter $\lambda$ is fixed to be the same as that in Step 1 from median regression. We consider three different simulation designs at $\gamma=0.5$ with $n=200$. The first model is a homoskedastic model by setting $\xi_{01}=(1,0,\ldots,0)$ in the baseline design. The second model is the same as the heteroskedastic median regression in Table (ref). The third model is a fat-tail model, where $U_i$ follows a Cauchy distribution with a scale parameter 0.25 while keeping the heteroskedastic design as the second model. The mean regression method shows slight over-selection but its performance looks reasonable in the homoskedastic model. However, the method in lee2012lasso is not robust to the heteroskedastic errors, which we can observe in Panel B of Table (ref). Furthermore, it cannot detect different nonzero coefficients at different quantile levels while the quantile method shows such a result in Table (ref). Finally, the quantile method works well when the error distribution follows a Cauchy distribution in Panel C of Table (ref). However, the mean regression method performs poorly with a Cauchy error distribution as the conditional mean function is not well-defined in this case.
Table (ref) shows the performance of the estimator when there does not exist any change point. We use the baseline design with $\gamma=0.75$ and set $\delta=(0,\ldots,0)$. As we are interested in the performance of $\widehat{\delta}$, we report the average number of parameters selected in $\widehat{\delta}$, the MSE of $\widehat{\delta}$, and the proportion of detecting no-change point (No-change Prop.). As predicted by the theory, all measures on $\widehat{\delta}$ indicate that the estimator (Step 3b) detects no-change point models quite well. Both $\mathbb{E}[J(\widehat{\delta})]$ and MSE of $\widehat{\delta}$ are quite low and no-change proportion is high. We can also observe much improvement in these measure when the sample size increases from $n=200$ to $n=400$.
In this subsection, we consider the case when the model contains low minimal signals in $\delta$. Specifically, we consider the median regression model and set $\beta_{0.5} = (0,0.34,0,\ldots,0)$ and $\delta_{0.5} = (0,1,1/2,1/4,1/8,1/16,0\ldots,0)$. Table (ref) reports simulation results in this design. Note that the simulation design in Table (ref) is the same as that reported in Table (ref) except that $\delta_{0.5} = (0,1,0\ldots,0)$ in Table (ref). Therefore, we may view that the simulation design in Table (ref) satisfies the minimum signal condition, whereas that of this subsection does not. The simulation results in Table (ref) are consistent with asymptotic theory in Section (ref) and remarks in Section (ref) comparing estimators in step 3. The step 3b estimator performs better than the step 3a estimator in Table (ref), but it performs worse in Table (ref). Also note that the oracle proportion is zero for the step 3b estimator, which is expected given low signals in coefficients. Finally, it is important to note that the performance of the estimators of $\tau_0$ is good in terms of the MSE in the presence of low signals in $\delta$. The coverage probability of the confidence interval is much higher than the nominal level, which was not observed in previous simulations. Since the MSE and coverage probability between the infeasible oracle 2 estimator and other estimators are very similar, we interpret that the over-coverage result is not driven by high-dimensionality of regressors and variable selection. Perhaps this is due to a larger number of coefficients to estimate for the oracle 2 estimator, compared to Table (ref).
We have carried out additional Monte Carlo experiments. For the sake of brevity, we only report main findings here and show full results in the appendices. In Appendix (ref), we report simulation results when the change point $\tau_0$ and the distribution of $Q_i$ vary. In particular, we consider three different distributions of $Q_i$: Uniform$[0,1]$, $N(0,1)$, and $\chi^2(1)$. The change point parameter $\tau_0$ varies over $0.3,0.4,\ldots,0.7$ quantiles of each $Q_i$ distribution. We find that the performance of $\widehat{\tau}$ measured by the root-mean-squared error depends on the density of $Q_i$ distribution, as is expected from asymptotic theory. For instance, it is quite uniform over different $\tau_0$ when $Q_i$ follows Uniform$[0,1]$. However, when $Q_i$ follows $N(0,1)$ or $\chi^2(1)$, it performs better when $\tau_0$ is located at a point with a higher density of $Q_i$ distribution. Sensitivity analyses provided in Appendix (ref) show that the main simulation results are robust when we make changes over the range between $-15\%$ and $+15\%$ of the suggested tuning parameter values.
In summary, the proposed estimation procedure works well in finite samples and confirms the theoretical results developed earlier. The simulation studies show some advantages of the proposed estimator over the existing mean regression method, e.g.\ lee2012lasso. It also detects no-change-point models well without any pre-test. The main qualitative results are not sensitive to different simulation designs on $\tau_0$ and $Q_i$ as well as to some variation on tuning parameter values.
As an empirical illustration, we investigate the existence of tipping in the dynamics of racial segregation using the dataset constructed by card2008tipping. They show that the neighborhood's white population decreases substantially when the minority share in the area exceeds a tipping point (or threshold point), using U.S. Census tract-level data. lee2012testing develop a test for the existence of threshold effects and apply their test to this dataset. Different from these existing studies, we consider a high-dimensional setup by allowing both possibly highly nonlinear effects of the main covariate (minority share in the neighborhood) and possibly higher-order interactions between additional covariates.
We build on the specifications used in card2008tipping and lee2012testing to choose the following median regression with a constant shift due to the tipping effect:
where for census tract $i$, the dependent variable $Y_i$ is the ten-year change in the neighborhood's white population, $Q_i$ is the base-year minority share in the neighborhood, and $X_i$ is a vector of six tract-level control variables and their various interactions depending on the model specification. Both $Y_i$ and $Q_i$ are in percentage terms. The basic six variables in $X_i$ include the unemployment rate, the log of mean family income, the fractions of single-unit, vacant, and renter-occupied housing units, and the fraction of workers who use public transport to travel to work. The function $g(\cdot)$ is approximated by the cubic b-splines with 15 knots over equi-quantile locations, so the degrees of freedom are $19$ including an intercept term. In our empirical illustration, we use the census-tract-level sample of Chicago whose base year is 1980.
In the first set of models, we consider possible interactions among the six tract-level control variables up to six-way interactions. Specifically, the vector $X$ in the six-way interactions will be composed of the following 63 regressors, $$ \{X^{(1)}, \ldots,X^{(6)}, X^{(1)}X^{(2)},\ldots, X^{(5)}X^{(6)},\ldots, X^{(1)} X^{(2)}X^{(3)}X^{(4)}X^{(5)}X^{(6)}\}, $$ where $X^{(j)}$ is the $j$-th element among those tract-level control variables. Note that the lower order interaction vector (e.g.\ two-way or three-way) is nested by the higher order interaction vector (e.g.\ three-way or four-way). The total number of regressors varies from 26 (19 from b-splines, 6 from $X_i$ and $1\{Q_i>\tau\}$) when there is no interaction to 83 when there are full six-way interactions. In the next set of models, we add the square of each tract-level control variable and generate similar interactions up to six. In this case the total number of regressors varies from 32 to 2,529. For example, the number of regressors in the largest model consists of $\#(\text{b-spline basis})+\#(\text{indicator function}) + \#(\text{interactions up to six-way out of 12})=19 + 1 + \sum_{k=1}^6 \binom {12}k = 2,529$. This number is much larger than the sample size ($n=1,813$).
Table (ref) summarizes the estimation results at the 0.25, 0.5, and 0.75 quantiles, respectively. We report the total number of regressors in each model and the number of selected regressors in Step 3b. The change point $\tau$ is estimated by the grid search over 591 equi-spaced points in $[1,60]$. The lower bound value $1 \%$ corresponds to the 1.6 sample percentile of $Q_i$ and the upper bound value $60 \%$, which is about the upper sample quartile of $Q_i$, is the same as one used in card2008tipping. In this empirical example, we report the estimates of $\tau_0$ and the confidence intervals updated after Step 3b (that is, $\tau$ is re-estimated using the estimates of $\alpha_0$ in Step 3b). If this estimate is different from the previous one in Step 2, then we repeat Step 3b and Step 2 until it converges.
The estimation results suggest several interesting points. First, at each quantile, the proposed method selects sparse representations in all model specifications even when the number of regressors is relatively large. Furthermore, the number of selected regressors does not grow rapidly when we increase the number of possible covariates. It seems that the set of selected covariates overlaps across different dictionaries at each quantile. See Appendix (ref) for details on selected regressors. Second, the estimation results are different across different quantiles, indicating that there may exist heterogeneity in this application. The confidence intervals for $ \tau_0$ at the 0.25 quantile are quite tight in all cases and they provide convincing evidence of the tipping effect. If we look at the case of six-way interactions with 12 control variables, the estimated tipping point is 5.65% and the estimated jump size is $-5.50\%$. However, this strong tipping effect becomes weaker at the 0.50 and 0.75 quantiles as shown either by wider confidence intervals or by the zero jump size, i.e.\ $\widehat{\delta}=0$.
We now compare the estimation results from quantile regression with those from mean regression, which are reported in Table (ref) (full estimation results are in Appendix (ref)). We show two kinds of mean regression estimates: one with the untrimmed original data and the other with the trimmed data for which we drop top and bottom $5\%$ observations based on $\{Y_i\}$. The estimated tipping points are the same between the two datasets but the estimated jump size is much larger with the original data. Figure (ref) shows the fitted values over $Q_i$ at the sample mean of the six basic covariates. They are from the model of six-way interactions with 12 control variables and the vertical line indicates the location of a tipping point. The left panel of Figure (ref) compares the results between the mean and median regression results (without trimming the data) and the right panel shows the the interquartile range of the conditional distribution of $Y_i$ as a function of $Q_i$ given other regressors. It can be seen that the mean regression estimates are much more volatile around the tipping point than the median regression estimates, although the estimated tipping point is the same. In Figure (ref), we compare the mean regression estimates with and without trimming. Removing observations with top and bottom 5% $Y_i$'s stablize the estimates, thus demonstrating that the median regression estimates have the built-in feature that they are more stable with outliers of $Y_i$ than the mean estimates. Finally, looking at the right panel of Figure (ref), we can see that the 25 percentile of the conditional distribution drops at the tipping point of 5.65% but no such change at the 75% quantile. This shows that the quantile regression estimates can provide insights into distributional threshold effects in racial segregation.
In summary, this empirical example shows that the proposed method works well in the real empirical setup and is robust to outliers compared to the mean regression approach. The estimation results also confirm that there exists a tipping point in the racial segregation at the 0.25 quantile and that the tipping effect is heterogeneous over different quantiles.
In this paper, we have developed $\ell_1$-penalized estimators of a high-dimensional quantile regression model with an unknown change point due to a covariate threshold. We have shown among other things that our estimator of the change point achieves an oracle property without relying on a perfect covariate selection, thereby avoiding the need for the minimum level condition on the signals of active covariates. We have illustrated the usefulness of our estimation methods via Monte Carlo experiments and an application to tipping in the racial segregation.
In a recent working paper, leonardi2016 consider a high-dimensional mean regression model with multiple change points whose number may grow as the sample size increases. They have proposed a binary search algorithm to choose the number of change points. It is an important future research topic to develop a computationally efficient algorithm to detect multiple changes for high-dimensional quantile regression models.