# Regularized quantile regression under heterogeneous sparsity with
application to quantitative genetic traits

Qianchuan He111Equally contributed<EMAIL_ADDRESS>Public Health Sciences
Division, Fred Hutchinson Cancer Research Center, Seattle, WA 98109, U.S.A.
Linglong Kong222Equally contributed. Department of Mathematical and
Statistical Sciences, University of Alberta, Edmonton, AB Canada T6G 2G1
Yanhua Wang  School of Mathematics, Beijing Institute of Technology, Beijing,
100081, China  Sijian Wang  Department of Biostatistics & Medical Informatics,
The University of Wisconsin-Madison, Madison, WI 53706, U.S.A.  Timothy A.
Chan  Human Oncology and Pathogenesis Program, Memorial Sloan-Kettering Cancer
Center, New York, NY 10065, U.S.A.  Eric Holland  Human Biology Division, Fred
Hutchinson Cancer Research Center, Seattle, WA 98109, U.S.A.

###### Abstract

Genetic studies often involve quantitative traits. Identifying genetic
features that influence quantitative traits can help to uncover the etiology
of diseases. Quantile regression method considers the conditional quantiles of
the response variable, and is able to characterize the underlying regression
structure in a more comprehensive manner. On the other hand, genetic studies
often involve high-dimensional genomic features, and the underlying regression
structure may be heterogeneous in terms of both effect sizes and sparsity. To
account for the potential genetic heterogeneity, including the heterogeneous
sparsity, a regularized quantile regression method is introduced. The
theoretical property of the proposed method is investigated, and its
performance is examined through a series of simulation studies. A real dataset
is analyzed to demonstrate the application of the proposed method.

###### keywords:

Heterogeneous Sparsity; Quantitative Traits; Variable Selection; Quantile
Regression; Genomic Features

††journal: Computational Statistics & Data Analysis

##  1 Introduction

In many genetic studies, quantitative traits are collected for studying the
associations between the traits and certain genomic features. For example,
body mass index, lipids and blood pressure have been investigated with respect
to single nucleotide polymorphisms (SNPs) [Avery et al., 2011]. With the rapid
progress of high-throughput genome technology, new types of quantitative
traits have emerged and attracted considerable research interest, such as the
gene expression, DNA methylation, and protein quantification [Landmark-Høyvik
et al., 2013]. The analysis of these quantitative traits yields new insight
into biological processes and sheds light on the genetic basis of diseases.

Typically, quantitative genetic traits are analyzed by least-square based
methods, which seek to estimate the $E(Y|Z)$, where $Y$ is the trait and $Z$
is the set of covariates of interest. Quantile regression [Koenker and
Bassett, 1978] instead considers the conditional quantile function of $Y$
given $Z$, $Q_{\tau}(Y|Z)$, at a given $\tau\in(0,1)$. When $\tau$ is fixed at
0.5, quantile regression is simply the median regression, which is well known
to be more robust than the least-square estimation. In examining quantiles at
different $\tau$, quantile regression provides a more complete picture of the
underlying regression structure between $Y$ and $Z$.

Like least-square methods, traditional quantile regression methods only
consider a handful of covariates. With the emergence of high-dimensional data,
penalized quantile regression methods have been developed in recent years, and
can be broadly classified into two classes. The first class seeks to harness
the information shared among different quantiles to jointly estimate the
regression coefficients. Jiang et al. [2014] proposed two novel methods, Fused
Adaptive Lasso (FAL) and Fused Adaptive Sup-norm estimator (FAS), for variable
selection in interquantile regression. FAL combines the LASSO penalty
[Tibshirani, 1996] with the fused lasso penalty [Tibshirani et al., 2005],
while FAS imposes the grouped sup-norm penalty and the fused lasso penalty.
The fused lasso penalty delivers the effect of smoothing the regression slopes
of adjacent quantiles, hence FAL and FAS can be used to identify common
quantile slopes if such a smoothing property is desired. Zou and Yuan [2008]
adopted an $F_{\infty}$ penalty, which either eliminates or retains all the
regression coefficients for a covariate at multiple quantiles. The method by
Jiang et al. [2013] seeks to shrink the differences among adjacent quantiles
by resorting to the fused lasso penalty; however, this method does not perform
variable selection at the covariate level, i.e., it does not remove any
covariates from the model.

The second class of methods focuses on a single quantile at a time. Koenker
[2004] imposed a LASSO penalty on the random effects in the mixed-effect
quantile regression model; Li and Zhu [2008] adopted the LASSO penalty; Wu and
Liu [2009] explored the SCAD penalty [Fan and Li, 2001] and the adaptive LASSO
penalty [Zou, 2006], and proved the selection consistency and the normality of
the proposed estimators (for a fixed dimension of the covariates). Wang et al.
[2012] investigated several penalties for quantile regression under the
scenario of $p>n$, i.e., the dimension is larger than the sample size, and
proved the selection consistency of their proposed methods through a novel use
of the subgradient theory. A recent approach by Peng et al. [2014] shares
characteristics with both classes; its loss function targets a single $\tau$,
while its penalty borrows information across different quantiles; their
proposed penalty was shown to achieve more accurate results than the one that
uses information from only a single quantile.

If the regression coefficients associated with a covariate are treated as a
group, then some groups may be entirely zero and some other groups may be
partially zero. Thus, sparsity can occur both at the group level and within
the group level, and we call this type of sparsity as the heterogeneous
sparsity. In this paper, we propose an approach that conducts joint variable
selection and estimation for multiple quantiles under the situation that $p$
can diverge with $n$. Our proposed method is able to achieve sparsity both at
the group level and within the group level. We note that FAL can potentially
yield sparsity at the two levels, but this approach has not been evaluated
under the scenario where the dimension $p$ is high. To the best of our
knowledge, this is the first paper that explicitly investigates the
heterogeneous sparsity for quantile regression. We show that our method tends
to be more effective than the compared methods in handling heterogeneous
sparsity when the dimension is high. We also provide theoretical justification
for the proposed method. The paper is organized as the follows. In Section 2,
we describe the proposed method and the implementation details. In Section 3,
we prove the theoretical properties of the proposed method. In Section 4, we
show the results of simulation studies regarding several related methods, and
in Section 5, we present an example of real data analysis for the proposed
method.

##  2 Method

###  2.1 Data and Model

Let $Z$ represent the vector consisting of $p$ covariates, such as SNPs or
genes. Let $\gamma_{\tau}$ be the $p$-dimension coefficient vector at the
$\tau$th quantile. Let $Y$ be the random variable that denotes the phenotype
we are interested in, such as quantitative traits in genetic studies. For a
given $\tau\in(0,1)$, the linear quantile regression model is known as

$Q_{\tau}(Y|Z)=\gamma_{0}+Z^{\rm T}\gamma_{\tau},$

where $\gamma_{0}$ is the intercept, and $Q_{\tau}(Y|Z)$ is the $\tau$-th
conditional quantile of $Y$ given $Z$, that is, $P_{Y|Z}(Y\leq
Q_{\tau}(Y|Z))=\tau$.

The dimension $p$ can be potentially very high in genomic studies, but
typically it is assumed that only a limited number of genomic features
contribute to the phenotype. For this reason, one needs to find a sparse
estimation of $\gamma_{\tau}$ to identify those important genomic features. On
the other hand, we also wish to consider multiple quantile levels
simultaneously so that information shared among different quantile levels can
be utilized. To this end, we propose the following model for the joint
estimation of the regression coefficients for multiple quantiles. Given $M$
quantile levels, $0<\tau_{1}<\cdots<\tau_{M}<1$, our linear quantile
regression model is defined as, for $\tau_{m}$ $(m=1,\cdots,M)$,

$\displaystyle Q_{\tau_{m}}(Y|Z)=\gamma_{m0}+Z^{\rm T}\gamma_{\tau_{m}}.$ (1)

where $\gamma_{m0}$ is the intercept, and $\gamma_{\tau_{m}}$ is the
$p$-dimension coefficient vector. For ease of notation, we write
$\gamma_{m}=\gamma_{\tau_{m}}$. For the above model, we further define
$\gamma\equiv(\gamma_{1}^{\rm T},\cdots,\gamma_{M}^{\rm T})^{\rm T}$ with
$\gamma_{m}=(\gamma_{m1},\cdots,\gamma_{mp})^{\rm T}$, and intercept parameter
$\gamma_{0}\equiv(\gamma_{10},\cdots,\gamma_{M0})^{\rm T}$.

We now focus on the sample version of the model (1). Let $\\{(Y_{i},Z_{i}^{\rm
T})^{\rm T}\\}_{i=1}^{n}$ be an i.i.d. random sample of size $n$ from
population $(Y,Z^{\rm T})^{\rm T}$, where
$Z_{i}=(Z_{i1},Z_{i2},\cdots,Z_{ip})^{\rm T}$. The sample quantile loss
function is defined as

$\displaystyle
Q_{n}(\gamma_{0},\gamma)=\sum_{m=1}^{M}\sum_{i=1}^{n}\,\rho_{m}(Y_{i}-Z_{i}^{\rm
T}\gamma_{m}-\gamma_{m0})$

where $\rho_{m}(u)=u(\tau_{m}-I(u<0))$ is the quantile check loss function
with $I(\cdot)$ being the indicator function. To introduce sparsity to the
model, we add to the loss function a penalty function

$\displaystyle
P_{n}(\gamma)=n\lambda_{n}\sum_{j=1}^{p}\left(\sum_{m=1}^{M}\,\omega_{mj}|\gamma_{mj}|\right)^{\frac{1}{2}},$

where $\lambda_{n}$ is the tuning parameter, and
$\omega_{n}=(\omega_{mj}:1\leq m\leq M,1\leq j\leq p)$ is the weight vector
whose component $\omega_{mj}>0$ is the weight of parameter $\gamma_{mj}$. Note
that the penalty is a nonconvex function. It essentially divides the
regression coefficients into $p$ groups, and each group consists of $M$
parameters associated with the $j$th covariate. The motivation is that, while
each quantile may have its own set of regression parameters, we wish to borrow
strengths from each quantiles to select covariates that are important across
all quantiles as well as those covariates that are important to only some of
the quantiles. This type of nonconvex penalty has been considered in the Cox
regression model and other settings; see Wang et al. [2009] for an example,
but to the best of our knowledge has not been studied in the quantile
regression model. We can choose $\omega_{mj}=\tilde{\gamma}_{mj}^{-1}$, where
$\tilde{\gamma}_{mj}$ is some consistent estimate for $\gamma_{mj}$. For
example, we may use the estimates from the unpenalized quantile regression
conducted at each individual quantile level. When $p<n$ but is fixed, the
consistency of the unpenalized estimates has been proved by Koenker and
Bassett [1978]. When $p<n$ but is diverging with $n$, the estimates from
unpenalized quantile regression are consistent by adapting to Lemma A.1 of
Wang et al. [2012]. Thus, our objective function is defined as

$\sum_{m=1}^{M}\sum_{i=1}^{n}\,\rho_{m}(Y_{i}-Z_{i}^{\rm
T}\gamma_{m}-\gamma_{m0})+n\lambda_{n}\sum_{j=1}^{p}\left(\sum_{m=1}^{M}\omega_{mj}|\gamma_{mj}|\right)^{\frac{1}{2}}.$

For the sake of convenience, we define $\theta=(\theta_{1}^{\rm
T},\cdots,\theta_{M}^{\rm T})^{\rm
T}\quad\text{with}\quad\theta_{m}=(\gamma_{m0},\,\gamma_{m}^{\rm T})^{\rm
T},\quad m=1,\cdots,M,$ and the corresponding parameter space by
$\Theta_{n}\subset\mathbb{R}^{M(p+1)}$. Further define $U_{i}=(1,Z_{i}^{\rm
T})^{\rm T}$, $i=1,\cdots,n$. Then, $Q_{n}(\gamma_{0},\gamma)$ can be written
as $Q_{n}(\theta)$. Emphasizing that $\gamma$ is a subset of $\theta$, we can
write the objective function as

$\displaystyle L_{n}(\theta)$ $\displaystyle\equiv
Q_{n}(\theta)+P_{n}(\gamma)=\sum_{m=1}^{M}\sum_{i=1}^{n}\,\rho_{m}(Y_{i}-U_{i}^{\rm
T}\theta_{m})+n\lambda_{n}\sum_{j=1}^{p}\left(\sum_{m=1}^{M}\omega_{mj}|\gamma_{mj}|\right)^{\frac{1}{2}}.$
(2)

Let $\hat{\theta}$ be a local minimizer of $L_{n}(\theta)$ in (2) for
$\theta\in\Theta_{n}$. Because the heterogeneity of sparsity is explicitly
taken into account in this model, we name our proposed method as Heterogeneous
Quantile Regression (Het-QR). Our model can be modified to accommodate
different weights for the losses at different quantiles. That is,
$Q_{n}(\theta)$ may take the form
$\sum_{m=1}^{M}\pi_{m}\sum_{i=1}^{n}\,\rho_{m}(Y_{i}-U_{i}^{\rm
T}\theta_{m})$, where $\pi_{m}$ is the weight for the $m$th quantile. Some
examples on the choice of weight $\pi_{m}$ can be found in Koenker [2004] and
Zhao and Xiao [2014].

###  2.2 Implementation

We design the following algorithm to implement the proposed method. First, we
show in the Appendix that the objective function can be transformed into

${\rm
argmin}_{\theta,\xi}\sum_{m=1}^{M}\sum_{i=1}^{n}\,\rho_{m}(Y_{i}-U_{i}^{\rm
T}\theta_{m})+\lambda_{1}\sum_{j=1}^{p}\xi_{j}+\sum_{j=1}^{p}\xi_{j}^{-1}\left(\sum_{m=1}^{M}\omega_{mj}|\gamma_{mj}|\right),$

where $\xi=(\xi_{1},\ldots,\xi_{p})$ are newly introduced nonnegative
parameters. Then, the new objective function can be solved by the following
iterative algorithm:
Step 1: We first fix $\theta$ to solve $\xi_{j},j=1,...,p$. To this end,
$\xi_{j}$ has a closed-form solution. That is,
$\hat{\xi}_{j}=\left(\sum_{m=1}^{M}\omega_{mj}|\gamma_{mj}|\right)^{1/2}\lambda_{1}^{-1/2},j=1,...,p$.
Step 2: We fix $\xi_{j},j=1,...,p$, to solve $\theta$. That is, we aim to
solve

${\rm argmin}_{\theta}\sum_{m=1}^{M}\sum_{i=1}^{n}\,\rho_{m}(Y_{i}-U_{i}^{\rm
T}\theta_{m})+\sum_{j=1}^{p}\sum_{m=1}^{M}\xi_{j}^{-1}\omega_{mj}|\gamma_{mj}|.$

We can formulate this objective function as a linear program and derive its
dual form (see Appendix for details), then the optimization can be conducted
by recoursing to existing linear programming packages; we utilize the Quantreg
R package [Koenker, 2015].
Step 3: Iterate step 1 and step 2 until convergence. Due to the nonconvexity
of the penalty function, the estimate is a local minimizer.

##  3 Theoretical Properties

Now we investigate the asymptotic properties of the proposed method. FAL and
FAS considered $p$ to be fixed. We study the situation where $p$ can diverge
with $n$. Let the true value of $\theta$ be $\theta^{*}$, where the
corresponding true values of $\gamma_{m0},\gamma_{mj},\gamma$ are
$\gamma_{m0}^{*},\gamma_{mj}^{*},\gamma^{*}$, respectively. Let the number of
nonzero elements in $\gamma^{*}$ be $s$. To emphasize that $s$ and $p$ can go
to infinity, we use $s_{n}$ and $p_{n}$ when necessary. For Theorems 1 and 2
(to be shown), $p_{n}$ is at the order lower than $O(n^{1/2})$; for Theorem 3,
$p_{n}$ is at the order lower than $O(n^{1/6})$.

We define some index sets to be used in our theorems. Let
$\mathcal{N}=\\{(m,j):1\leq m\leq M,1\leq j\leq p_{n}\\}$. For the true
parameters, define the oracle index set
$I=\\{(m,j)\in\mathcal{N}:\gamma_{mj}^{*}\neq 0\\}$ and its complementary set
$II=\\{(m,j)\in\mathcal{N}:\gamma_{mj}^{*}=0\\}$. Assume that $I$ has
cardinality $|I|=s_{n}$.

We define some notations used for our theorems. Define $d_{nI}=\max_{(m,j)\in
I}\omega_{mj}$ and $d_{nII}=\min_{(m,j)\in
II}\omega_{mj}\bigl{(}\max_{(m,j)\in\mathcal{N}}\omega_{mj}\bigr{)}^{-\frac{1}{2}}$.
Define $\theta_{I}^{*}$ and $\hat{\theta}_{I}$ as the subvectors of the
vectors $\theta^{*}$ and $\hat{\theta}$ corresponding to the oracle index set
$I$, respectively. For every fixed $1\leq m\leq M$, define the index set
$I_{m}=\\{1\leq j\leq p_{n}:\gamma_{mj}^{*}\neq 0\\}$. Let

$\displaystyle\mathbf{\Sigma}_{n}$
$\displaystyle=(\mathbf{\Sigma}_{lm})_{M\times
M}\quad\text{with}\quad\mathbf{\Sigma}_{lm}=(\min(\tau_{m},\tau_{l})-\tau_{m}\tau_{l})\,E(U_{ilI}U_{imI}^{\rm
T}),$ (3) $\displaystyle\mathbf{B}_{nI}$
$\displaystyle=Diag(B_{1},B_{2},\cdots,B_{M})\quad\text{with}\quad
B_{m}=\sum_{i=1}^{n}\,f(U_{i}^{\rm
T}\theta_{m}^{*}|\,Z_{i})\,U_{imI}U_{imI}^{\rm T},$ (4)

where $U_{imI}=(1,Z_{imI}^{\rm T})^{\rm T}$ with $Z_{imI}$ being the subvector
of $Z_{i}$ corresponding to index set $I_{m}$.

Let $F(y|z)$ and $f(y|z)$ be the conditional distribution function and the
conditional density function of $Y$ given $Z=z$, respectively. For any proper
square matrix $A$, let $\lambda_{min}(A)$ and $\lambda_{max}(A)$ denote the
minimum and maximum eigenvalue of $A$, respectively.

Before stating the main theorems, we need the following regularity conditions
labeled by $\mathcal{L}$:

  1. (L1)

The conditional density $f(y|z)$ has first order derivative $f^{\prime}(y|z)$
with respect to $y$; And $f(y|z)$ and $f^{\prime}(y|z)$ are uniformly bounded
away from $0$ and $\infty$ on the support set of $Y$ and the support set of
$Z$;

  2. (L2)

For randon sample $Z_{i}=(Z_{i1},Z_{i2},\cdots,Z_{ip})^{\rm T}$, $1\leq i\leq
n$, there exists a postive constant $C_{1}$ such that $\max_{1\leq i\leq
n,1\leq j\leq p}|Z_{ij}|\leq C_{1}$;

  3. (L3)

For $U_{i}=(1,Z_{i}^{\rm T})^{\rm T}$, $i=1,\cdots,n$, let
$\mathbf{S}_{n}=\sum_{i=1}^{n}\,U_{i}U_{i}^{\rm T}$. There exist positive
constants $C_{2}<C_{3}$ such that
$C_{2}\leq\lambda_{min}(n^{-1}\mathbf{S}_{n})\leq\lambda_{max}(n^{-1}\mathbf{S}_{n})\leq
C_{3}$;

  4. (L4)

The dimension $s_{n}$ satisfies that $s_{n}=a_{0}n^{\alpha_{0}}$, and the
dimension $p_{n}$ satisfies that $p_{n}=a_{1}n^{\alpha_{1}}$, where
$0<\alpha_{0}<\alpha_{1}<\frac{1}{2}$, and $a_{0}$ and $a_{1}$ are two
positive constants;

  5. (L5)

The matrix $\mathbf{B}_{n}$ given in $\eqref{0-anbn}$ (see Appendix) satisfies
that
$C_{4}\leq\lambda_{min}(n^{-1}\mathbf{B}_{n})\leq\lambda_{max}(n^{-1}\mathbf{B}_{n})\leq
C_{5}$, where $C_{4}$ and $C_{5}$ are positive constants.

  6. (L6)

The matrix $\mathbf{\Sigma}_{n}$ satisfies that
$\lambda_{min}(\mathbf{\Sigma}_{n})\geq C_{6}$, where $C_{6}$ is a positive
constant.

Conditions (L1)-(L3) and (L5)-(L6) are seen in typical theoretical
investigation of quantile regression. Condition (L4) specifies the magnitude
of $s_{n}$ and $p_{n}$ with respect to the sample size. Under the
aforementioned regularity conditions, we present the following three theorems.
The proof is relegated to the Appendix. Define $\theta_{II}^{*}$ and
$\hat{\theta}_{II}$ as the subvectors of the vectors $\theta^{*}$ and
$\hat{\theta}$ corresponding to the index set $II$, respectively. Clearly,
$\theta_{II}^{*}=0$. Due to the nonconvexity of the penalty function, all the
following theorems and their proofs (in the Appendix) are regarding to a local
minimizer of the objective function.

Theorem 1. Under conditions (L1)-(L5), if
$\lambda_{n}d_{nI}^{\frac{1}{2}}=o(s_{n}^{-\frac{3}{4}}p_{n}^{\frac{3}{4}}n^{-\frac{3}{4}})$,
then the estimator $\hat{\theta}$ of $\theta^{*}$ exists, is a local
minimizer, and satisfies the estimation consistency that
$\|\hat{\theta}-\theta^{*}\|_{2}=O_{p}(n^{-\frac{1}{2}}p_{n}^{\frac{1}{2}})$.

Theorem 1 shows that the proposed method is consistent in parameter
estimation. The convergence rate $O_{p}(n^{-\frac{1}{2}}p_{n}^{\frac{1}{2}})$
is typical for the settings where $p$ diverges with $n$.

Theorem 2. Under conditions (L1)-(L5), if
$\lambda_{n}d_{nI}^{\frac{1}{2}}=o(n^{-\frac{3}{4}}s_{n}^{-\frac{3}{4}}p_{n}^{\frac{3}{4}})$
and $n^{-\frac{1}{2}}p_{n}^{\frac{1}{2}}=o(\lambda_{n}\,d_{nII})$, then
$P(\hat{\theta}_{II}=0)\to 1$.

Theorem 2 indicates that our method can distinguish the truly zero
coefficients from the nonzero coefficients with probability tending to 1. It
can be seen that the penalty weight $d_{nII}$ plays a critical role in the
property of selection consistency.

Theorem 3. Under conditions (L1)-(L3) and (L5)-(L6), if
$\lambda_{n}d_{nI}^{\frac{1}{2}}=O(n^{-1}p_{n}^{\frac{1}{2}})$, and
$n^{-\frac{1}{2}}p_{n}^{\frac{1}{2}}=o(\lambda_{n}\,d_{nII})$, and the powers
of $s_{n}$ and $p_{n}$ in condition $(L4)$ satisfy
$0<\alpha_{0}<\alpha_{1}<\frac{1}{6}$, then for any unit vector
$b\in\mathbb{R}^{M+s_{n}}$ we have

$\displaystyle(nb^{\rm T}\mathbf{\Sigma}_{n}b)^{-\frac{1}{2}}\,b^{\rm
T}\mathbf{B}_{nI}\bigl{(}\hat{\theta}_{I}-\theta_{I}^{*}\bigr{)}\to N(0,1).$

Theorem 3 suggests that the estimated nonzero coefficients have the asymptotic
normality. Heuristically, for given $n,\lambda_{n}$ and $\omega_{mj}$, the
considered penalty has its slope tending to infinity when $\gamma_{mj}$ goes
to 0, thus the penalty tends to dominate small $\gamma_{mj}$. On the other
hand, when $\lambda_{n}$ is sufficiently small, the penalty has little impact
on the estimation of relatively large $\gamma_{mj}$. These properties, in
combination with proper choice of the tuning parameter, play major roles in
the oracle property of the proposed estimator. The oracle property for
coefficients within a group is mainly due to the penalty weights, which put
large penalty on small coefficients (and small penalty on large coefficients).

##  4 Simulation Studies

We conduct simulation studies to evaluate the proposed method along with the
following methods: the QR method, which applies quantile regression to each
individual quantile level without any variable selection; the QR-LASSO method,
which adopts the $L_{1}$-penalized quantile regression to each quantile level;
the QR-aLASSO method, which imposes the adaptive LASSO penalty to each
quantile level [Wu and Liu, 2009]; the FAL and the FAS method [Jiang et al.,
2014]. Both FAL and FAS contain a fused-LASSO type of penalty, which
encourages the equality of the regression coefficients among different
quantiles. FAL allows within-group sparsity, while FAS generates sparsity only
at the group level. For Het-QR, we set the penalty weight $\omega_{mj}$ to be
the inverse of the estimate from the unpenalized quantile regression (unless
specified otherwise).

We first consider a model where important covariates have nonzero regression
coefficients across all (or almost all) quantiles. We simulate 6 independent
covariates, each of which follows the uniform(0,1) distribution. Then we
simulate the trait as

$Y=1.0+\beta_{1}Z_{1}+\beta_{2}Z_{2}+\beta_{6}Z_{6}+\kappa Z_{6}\epsilon,$

where $\beta_{1}=1,\beta_{2}=1,\beta_{6}=2$, $\kappa=2$ and $\epsilon\sim
N(0,1)$. Under this set up, $Z_{1}$ and $Z_{2}$ have constant regression
coefficients across all quantiles, while $Z_{6}$’s regression coefficient is
determined by $2+2\times\Phi^{-1}(\tau)$, which varies among different
quantiles. That is, the $\tau$th quantile of $Y$ given $Z_{1},Z_{2}$ and
$Z_{6}$ is

$Q_{\tau}\left(Y|Z_{1},Z_{2},Z_{6}\right)=1.0+Z_{1}+Z_{2}+\left(2+2\times\Phi^{-1}(\tau)\right)Z_{6}.$

This model is in line with the model considered by Jiang et al. [2014]. All
the other 3 covariates, $Z_{3},Z_{4}$ and $Z_{5}$, have no contribution to
$Y$. The sample size $n$ is set to 500. To select the tuning parameter, we
follow the lines of Mazumder et al. [2011] and Wang et al. [2012] to generate
another dataset with sample size of $10n$, and then pick the tuning parameter
at which the check loss function is minimized. The total number of simulations
for each experiment is $100$.

We consider various criteria to evaluate the performance of the compared
methods, such as the model size and the parameter estimation error (PEE). The
model size refers to the number of estimated non-zero coefficients among the
$M$ quantile levels. The PEE is calculated by
$\sum_{m=1}^{M}\sum_{j=1}^{p}|\hat{\gamma}_{mj}-\gamma^{*}_{mj}|/M$. To
evaluate the prediction error, we simulate an independent dataset,
$(Y_{pred},Z_{pred})$, with sample size of 100n, and then calculate the
F-measure (FM) [Gasso et al., 2009], the quantile prediction error (QPE) and
the prediction error (PE). The FM is equal to $2\times S_{a}/M_{a}$, where
$S_{a}$ is the number of truly nonzero slopes being captured, and $M_{a}$ is
the sum of the estimated model size and the true model size. The QPE is
defined as the sample version of the
$\sum_{m=1}^{M}(Q_{\tau_{m}}(Y_{pred}|Z_{pred})-Z_{pred}^{\rm
T}\hat{\gamma}_{\tau_{m}}-\hat{\gamma}_{m0})^{2}/M$, averaged across all the
subjects. The PE is defined as $Q_{n}(\hat{\theta})/n$, i.e., the check loss
averaged across all considered quantiles for all the test samples, evaluated
on $(Y_{pred},Z_{pred})$.

For the purpose of illustration, we consider three quantiles,
$\tau=0.25,0.5,0.75$. The results are shown in the upper panel of Table 1. It
can be seen that when $p$ is 6, FAL has the lowest parameter estimation error
and FAS has the lowest PE, though the difference between these two methods and
the other compared methods is generally quite small. Next, we increase $p$ to
100 to evaluate the the methods under higher dimension. As shown in the lower
panel of Table 1, when $p$ is equal to 100, both FAL and FAS have deteriorated
performance; for instance, their model sizes tend to be twice (or more) as the
true model size and their PEE and PE are higher than Het-QR. This experiment
shows that the performance of FAL and FAS is suboptimal when the
dimensionality grows large; one potential explanation is that the penalties of
FAL and FAS may overemphasize the interquantile shrinkage, which make them
less efficient when many noise covariates are present. Further research is
merited. As to computation, we did not observe non-convergence for Het-QR in
our experiments.

Table 1. Comparison of Het-QR and other methods in the absence of within-group
sparsity (standard error of the sample mean shown in the parenthesis )

Method | Model-size | FM (%) | $PEE\times 100$ | $QPE\times 10^{3}$ | $PE\times 10^{3}$ |
---|---|---|---|---|---|---
| $p=6$ |
QR | 18 | - | 53.3(1.5) | 10.2(0.6) | 1041.2(0.6) |
QR-LASSO | 15.0(0.2) | 75(0.7) | 36.1(1.2) | 6.7(0.5) | 1038.4(0.6) |
QR-aLASSO | 10.9(0.2) | 91(0.7) | 25.6(0.9) | 5.8(0.4) | 1037.2(0.5) |
FAL | 11.1(0.2) | 91(0.9) | 25.3(1.0) | 6.2(0.5) | 1037.2(0.5) |
FAS | 12.1(0.2) | 86(0.9) | 26.4(1.1) | 5.9(0.4) | 1037.1(0.5) |
Het-QR | 9.6(0.1) | 97(0.6) | 26.0(0.9) | 6.5(0.5) | 1037.5(0.5) |
| $p=100$ |
QR | 300 | - | 1556.4(12.7) | 325.0(5.1) | 1242.1(2.7) |
QR-LASSO | 47.3(1.2) | 33(0.7) | 120.5(3.5) | 23.4(1.2) | 1052.4(0.9) |
QR-aLASSO | 16.8(0.4) | 72(1.1) | 46.9(1.9) | 10.5(0.8) | 1042.1(0.7) |
FAL | 17.7(0.6) | 70(1.4) | 41.9(1.8) | 10.2(0.8) | 1041.0(0.6) |
FAS | 23.7(0.7) | 58(1.2) | 58.0(2.4) | 13.9(0.9) | 1043.4(0.7) |
Het-QR | 9.3(0.1) | 99(0.4) | 29.5(1.2) | 8.4(0.6) | 1039.9(0.6) |

Next, we systematically evaluate the situation where within-group sparsity
exists. To introduce correlations into covariates, we simulate 20 blocks of
covariates, each block containing 5 correlated covariates. For each block, we
first simulate a multivariate normal distribution with mean being the unit
vector and covariance matrix following either the compound symmetry or the
auto-regressive correlation structure with correlation coefficient $\rho=0.5$;
next, we take the absolute value of the simulated random normal variables as
the covariates $Z$. The total number of covariates is 100. We specify the
conditional quantile regression coefficient function $\gamma(\tau)$ as
follows. For $\tau\in(0,0.3]$, the first 8 regression slopes for $Z$ are
$(0.5,0,0,0,0,0.6,0,0)$; for $\tau\in(0.3,0.7]$, the first 8 regression slopes
are $(0.5,0,0,0,0,0.6,0,0.7)$; for $\tau\in(0.7,1.0)$, the corresponding
slopes are $(0.6,0,0,0,0,0.7,0,0.7)$. All other regression slopes are 0. Thus,
the first and the sixth covariates are active among all quantiles, while the
eighth covariate is active only for the last two quantile levels. To generate
$Y$, we first simulate a random number $\tau\in\mbox{Uniform}(0,1)$, and then
determine the $\gamma(\tau)$ based on $\tau$; subsequently, we obtain

$Y=1.0+Z^{\rm T}\gamma(\tau)+F^{-1}(\tau),$

where $F^{-1}$ is the inverse cumulative function of some distribution $F$.
That is, the $\tau$th quantile of $Y$ given $Z$ is

$Q_{\tau}\left(Y|Z\right)=\begin{cases}1.0+F^{-1}(\tau)+0.5Z_{1}+0.6Z_{6}&\text{if
}0<\tau\leq 0.3\\\ 1.0+F^{-1}(\tau)+0.5Z_{1}+0.6Z_{6}+0.7Z_{8}&\text{if
}0.3<\tau\leq 0.7\\\ 1.0+F^{-1}(\tau)+0.6Z_{1}+0.7Z_{6}+0.7Z_{8}&\text{if
}0.7<\tau<1\end{cases}.$

We explore different distributions for $F$: the standard normal distribution,
the $T$-distribution with degrees of freedom equal to 3 ($T_{3}$), and the
exponential distribution with shape parameter equal to 1.

We first consider the normal distribution for $F$. The results are shown in
Table 2. Because no variable selection is conducted, QR has much larger PEE,
QPE, and PE than the other methods; for example, the PEE and QPE of QR are
more than 10 times higher than the compared methods. QR-LASSO, QR-aLASSO, FAL,
and FAS have more tamed model sizes, but still contain a number of noise
features. Het-QR yields a model that is closer to the true model, in which the
three considered quantiles contains 2, 2, and 3 nonzero slopes, respectively.
Het-QR also appears to have the highest FM, and lowest errors for parameter
estimation and prediction. Next, we consider the distribution to be the
$T_{3}$ (Table 3) and the exponential distribution (Table 4), and the results
show a similar pattern. These experiments indicate that Het-QR can handle
higher dimension as well as the heterogeneous sparsity better than the other
methods.

We finally consider the situation where $p>n$. While theoretical development
is still needed for this setting, our experiment is to evaluate the practical
performance of the proposed approach. We let $n=500$ and $p=600$. For
$\tau\in(0,0.3]$, the first 8 regression slopes for $Z$ are
$(0.6,0,0,0,0,0.6,0,0)$; for $\tau\in(0.3,0.7]$, the first 8 regression slopes
are $(0.6,0,0.8,0,0,0.7,0,0.8)$; for $\tau\in(0.7,1.0)$, the corresponding
slopes are $(0.8,0,0.8,0,0,0.8,0,1.0)$. In this scenario, $Z_{3}$ and $Z_{8}$
have zero coefficients for the first quantile, but nonzero coefficients for
the other two quantiles. That is, the $\tau$th quantile of $Y$ given $Z$ is

$Q_{\tau}\left(Y|Z\right)=\begin{cases}1.0+F^{-1}(\tau)+0.6Z_{1}+0.6Z_{6}&\text{if
}0<\tau\leq 0.3\\\
1.0+F^{-1}(\tau)+0.6Z_{1}+0.8Z_{3}+0.7Z_{6}+0.8Z_{8}&\text{if }0.3<\tau\leq
0.7\\\ 1.0+F^{-1}(\tau)+0.8Z_{1}+0.8Z_{3}+0.8Z_{6}+1.0Z_{8}&\text{if
}0.7<\tau<1\end{cases}.$

We omit QR, FAL and FAS because they are not designed to handle the setting of
‘$p>n$’. QR-LASSO can be directly applied to the data that have higher
dimension than sample size. For QR-aLASSO, we derive the penalty weights using
the estimates from QR-LASSO. For Het-QR, we first run Het-QR with the penalty
weights equal to 1 to obtain the initial estimators for
$\gamma_{m}^{*},m=1,\ldots,M$, and then use the inverse of the initial
estimators as the penalty weights; finally, we run Het-QR to obtain the
$\hat{\theta}$. The results are shown in Table 5. Het-QR tends to yield a
smaller model than the compared methods and have better performance in
estimating the regression coefficients as well as in prediction.

Table 2. Comparison of Het-QR and other methods for $p=100$ under the normal
distribution (standard error of the sample mean shown in the parenthesis )

Method | Model-size | FM (%) | $PEE\times 100$ | $QPE\times 10^{3}$ | $PE\times 10^{3}$ |
---|---|---|---|---|---|---
| Correlation structure: Auto-regressive |
QR | 300 | - | 1189.3(7.2) | 791.6(9.1) | 1606.5(3.1) |
QR-LASSO | 41.1(1.3) | 34(0.9) | 108.7(3.1) | 78.5(3.1) | 1353.1(1.2) |
QR-aLASSO | 17.8(0.5) | 63(1.2) | 60.4(2.5) | 49.2(3.2) | 1341.5(1.2) |
FAL | 20.3(0.8) | 60(1.6) | 59.5(2.6) | 49.9(3.5) | 1340.4(1.2) |
FAS | 24.8(1.1) | 53(1.5) | 85.1(2.6) | 86.4(3.1) | 1351.2(1.1) |
Het-QR | 8.6(0.1) | 96(0.7) | 34.7(1.6) | 32.6(2.8) | 1334.9(1.1) |
| Correlation structure: Compound symmetry |
QR | 300 | - | 1218.1(8.4) | 792.7(9.5) | 1606.1(3.4) |
QR-LASSO | 40.1(1.1) | 35.0(0.9) | 105.9(3.1) | 74.8(3.0) | 1351.9(1.3) |
QR-aLASSO | 18.9(0.6) | 61(1.3) | 61.2(2.7) | 48.0(3.1) | 1341.6(1.3) |
FAL | 19.9(0.8) | 61(1.4) | 61.0(2.9) | 54.4(4.0) | 1341.2(1.3) |
FAS | 24.1(1.1) | 55(1.6) | 83.8(2.8) | 89.3(3.2) | 1351.3(1.2) |
Het-QR | 8.9(0.2) | 94(0.8) | 34.6(1.9) | 31.4(2.9) | 1334.7(1.2) |

Table 3. Comparison of the Het-QR and other methods for $p=100$ under the
$T_{3}$ distribution (standard error of the sample mean shown in the
parenthesis )

Method | Model-size | FM (%) | $PEE\times 100$ | $QPE\times 10^{3}$ | $PE\times 10^{3}$ |
---|---|---|---|---|---|---
| Correlation structure: Auto-regressive |
QR | 300 | - | 1452.1(10.2) | 1149.1(15.2) | 2114.3(4.5) |
QR-LASSO | 39.4(1.3) | 35(0.9) | 123.6(3.3) | 103.7(3.7) | 1796.8(1.4) |
QR-aLASSO | 19.7(0.6) | 58(1.3) | 81.0(3.0) | 77.3(4.2) | 1787.7(1.6) |
FAL | 20.6(0.7) | 58(1.3) | 76.4(3.1) | 74.8(4.6) | 1786.2(1.6) |
FAS | 24.2(0.9) | 53(1.3) | 99.6(3.4) | 107.8(4.2) | 1795.6(1.6) |
Het-QR | 8.9(0.2) | 92(1.0) | 47.3(2.3) | 53.5(4.0) | 1779.4(1.5) |
| Correlation structure: Compound symmetry |
QR | 300 | - | 1487.4(11.2) | 1156.3(15.6) | 2115.2(4.7) |
QR-LASSO | 38.7(1.0) | 35(1.2) | 120.5(3.2) | 99.4(3.4) | 1795.5(1.6) |
QR-aLASSO | 20.3(0.7) | 56(1.3) | 81.6(3.3) | 76.0(4.1) | 1787.6(1.8) |
FAL | 22.5(0.9) | 56(1.5) | 82.5(3.9) | 78.1(5.1) | 1786.6(1.8) |
FAS | 25.6(1.0) | 51(1.3) | 99.1(3.5) | 107.0(4.1) | 1794.7(1.7) |
Het-QR | 9.1(0.3) | 91(1.1) | 48.6(2.9) | 55.4(4.7) | 1780.1(1.8) |

Table 4. Comparison of Het-QR and other methods for $p=100$ under the
exponential distribution (standard error of the sample mean shown in the
parenthesis )

Method | Model-size | FM (%) | $PEE\times 100$ | $QPE\times 10^{3}$ | $PE\times 10^{3}$ |
---|---|---|---|---|---|---
| Correlation structure: Auto-regressive |
QR | 300 | - | 1024.2(7.1) | 618.4(8.8) | 1446.4(3.1) |
QR-LASSO | 39.3(1.1) | 36(0.9) | 88.7(2.5) | 60.6(2.2) | 1220.9(1.0) |
QR-aLASSO | 17.1(0.5) | 66(1.2) | 47.9(1.9) | 37.0(2.3) | 1210.2(0.9) |
FAL | 18.9(0.7) | 63(1.5) | 43.0(1.9) | 33.3(3.0) | 1209.5(1.0) |
FAS | 26.9(1.2) | 51(1.7) | 74.5(2.3) | 79.1(3.5) | 1223.8(1.1) |
Het-QR | 8.7(0.1) | 96(0.6) | 28.8(1.3) | 24.2(1.8) | 1205.4(0.8) |
| Correlation structure: Compound symmetry |
QR | 300 | - | 1041.6(7.5) | 610.9(8.3) | 1444.4(2.9) |
QR-LASSO | 39.6(1.2) | 36(0.9) | 86.1(2.6) | 58.0(2.4) | 1220.2(1.1) |
QR-aLASSO | 17.0(0.5) | 66(1.3) | 46.6(2.0) | 35.7(2.3) | 1210.1(1.1) |
FAL | 18.7(0.7) | 63(1.4) | 41.8(2.0) | 31.0(3.0) | 1208.8(1.1) |
FAS | 27.6(1.4) | 51(1.8) | 74.3(2.5) | 79.6(3.3) | 1223.5(1.1) |
Het-QR | 8.7(0.1) | 96(0.7) | 27.2(1.3) | 22.3(1.9) | 1205.1(0.9) |

Table 5. Comparison of Het-QR and other methods for $p>n$ (standard error of
the sample mean shown in the parenthesis )

Method | Model-size | FM (%) | $PEE\times 100$ | $QPE\times 10^{3}$ | $PE\times 10^{3}$ |
---|---|---|---|---|---|---
| Correlation structure: Auto-regressive |
QR-LASSO | 67.6(1.8) | 25(0.6) | 199.0(4.2) | 215.5(6.0) | 1801.4(1.7) |
QR-aLASSO | 15.2(0.3) | 73(1.1) | 81.3(3.0) | 113.3(6.2) | 1768.8(1.6) |
Het-QR | 10.6(0.1) | 96(0.6) | 44.6(2.5) | 61.0(7.0) | 1753.9(1.6) |
| Correlation structure: Compound symmetry |
QR-LASSO | 62.9(1.8) | 27(0.7) | 182.4(4.2) | 199.8(5.9) | 1796.4(1.8) |
QR-aLASSO | 15.2(0.3) | 74(1.0) | 76.7(2.9) | 103.5(5.9) | 1766.6(1.6) |
Het-QR | 10.6(0.1) | 96(0.7) | 46.4(3.0) | 63.1(7.8) | 1754.2(1.7) |

##  5 Real Data Analysis

We collect 206 brain tumor patients each with 91 genes expression levels. All
patients were de-identified. All patients were diagnosed to have glioma, one
of the deadliest cancers among all cancer types. Indeed, many patients died
within 1 year after the diagnosis. Glioma is associated with a number of
genes. We focus on the PDGFRA gene, which encodes the alpha-type platelet-
derived growth factor receptor and has been shown to be an important gene for
brain tumors [Holland, 2000, Puputti et al., 2006]. We use this data set to
investigate how the expression of PDGFRA is influenced by other genes.

For demonstration, we set $\tau$ to 0.25, 0.5, and 0.75. For QR-LASSO, QR-
aLASSO, and Het-QR, we use cross-validation to ascertain the tuning parameter.
That is, (1) we divide the data into 3 folds; (2) we use 2 folds to build the
model and 1 fold to calculate the prediction errors, and this is done three
times; (3) we choose the $\lambda_{n}$ that minimizes the prediction error as
the best tuning parameter (we were not able to obtain an independent dataset
with sample size $10n$ to determine the tuning parameter; it would be
meaningful to compare the two procedures when such a dataset becomes available
in future). For FAL and FAS, we follow [Jiang et al., 2014] to use BIC and AIC
for determining the tuning parameter, and the corresponding methods are named
as FAL-BIC, FAL-AIC, FAS-BIC, FAS-AIC. Hence, in total 7 approaches are
compared. We examine the final model sizes after applying the 7 approaches.
Taken all the three quantiles together, the number of nonzero covariates of
the seven models are 89 (QR-LASSO), 47 (QR-aLASSO), 25 (Het-QR), 49 (FAL-BIC),
182 (FAL-AIC), 93 (FAS-BIC), 167 (FAS-AIC). For illustration, we list some of
the estimated regression coefficients in Table 6. It can be seen that the
coefficients for a given gene often differ among different quantiles. For a
better view of the regression coefficients among different quantiles, we plot
the estimated coefficients for the first 30 covariates (Figure 1). Table 6 and
Figure 1 show that most models (except FAS-BIC and FAS-AIC) demonstrate
heterogeneous sparsity, i.e., some covariates have nonzero effects in only one
or two of the three quantiles. FAS-BIC and FAS-AIC do not show this type of
sparsity due to the sup-norm penalty they adopt, as this penalty either
selects or removes a covariate for all the quantiles. FAL-AIC and FAS-AIC
models contain more nonzero estimates than FAL-BIC and FAS-BIC, consistent
with the fact that BIC favors smaller models than AIC. Compared to other
methods, Het-QR yields a smaller model which may be easier to interpret and
prioritize candidate genes for further functional study.

The covariates selected by Het-QR are shown in Table 7. Consistent with the
model assumptions, the estimated regression coefficients show heterogeneity
among quantiles. For example, the CDKN2C gene has zero coefficient at
$\tau=0.25$, and nonzero coefficients at $\tau=0.5$ and 0.75. In contrast,
some other genes, such as BMP2 and SLC4A4, have nonzero coefficients across
all the considered quantiles. This suggests that the expression of PDGFRA is
influenced by other genes in a delicate manner that may not be fully
characterized by least square methods or quantile regression methods that fail
to account for the genetic heterogeneity. CDKN2C encodes a cyclin-dependent
kinase, and BMP2 and SLC4A4 encode a bone morphogenetic protein and a sodium
bicarbonate cotransporter, respectively. This indicates that PDGFRA’s
expression is associated with genes with a wide spectrum of cellular
functions. The gene EGFR has non-positive regression coefficients, suggesting
that there may be some negative control between PDGFRA and EGFR. Future
biological studies may provide new insight on the gene regulation of PDGFRA.

One main purpose of variable selection is to apply the selected variables from
one dataset to other datasets to guide statistical analysis. Along this line,
we further collect the brain tumor data from the cancer genome atlas (TCGA)
project, which contains 567 subjects. We apply the models selected by
different methods from the training data to the TCGA data to assess the
prediction accuracy of the different models. We randomly split the TCGA data
into two halves, and use one half to estimate the regression coefficients and
the other half to calculate the prediction error; the prediction error is then
averaged across the two halves. We repeat the random-splitting 400 times, and
calculate the average of the prediction errors. Het-QR appears to have a
slightly lower prediction error than the other compared ones, but the
difference among the seven methods is generally small; in detail, the observed
prediction errors are 1.349 (QR-LASSO), 1.351 (QR-aLASSO), 1.345 (Het-QR),
1.362 (FAL-BIC), 1.513 (FAL-AIC), 1.355 (FAS-BIC), and 1.430 (FAS-AIC).

Figure 1: Graphic view of the first 30 regression coefficients estimated by
different methods. Estimates are thresholded at 0.4 and -0.4, and only nonzero
estimates are shown.

Table 6. A snapshot of the estimated regression coefficients (only 5
covariates are shown)

Gene | $\tau$ | QR-LASSO | QR-aLASSO | Het-QR | FAL-BIC | FAL-AIC | FAS-BIC | FAS-AIC |
---|---|---|---|---|---|---|---|---|---
| 0.25 |  |  |  |  | 0.10 |  | 0.06 |
POLR2A | 0.5 |  |  |  |  | 0.20 |  | 0.06 |
| 0.75 |  |  |  |  | 0.20 |  | 0.06 |
| 0.25 |  |  |  | 0.13 | 0.30 | 0.09 | 0.19 |
SDHA | 0.5 |  |  |  | 0.13 | 0.26 | 0.09 | 0.19 |
| 0.75 | 0.03 | 0.32 | 0.2 | 0.13 | 0.58 | 0.09 | 0.19 |
| 0.25 |  |  |  |  | 0.09 | 0.01 | 0.04 |
CDKN2A | 0.5 | 0.03 |  |  |  | 0.03 | 0.01 | 0.04 |
| 0.75 | 0.02 |  |  |  |  | 0.01 | 0.04 |
| 0.25 |  |  |  |  | 0.10 | 0.10 | 0.19 |
CDKN2C | 0.5 | 0.02 | 0.07 | 0.13 | 0.19 | 0.26 | 0.10 | 0.19 |
| 0.75 | 0.05 | 0.33 | 0.25 | 0.19 | 0.34 | 0.10 | 0.19 |
| 0.25 | 0.07 |  |  |  |  | 0.03 | 0.05 |
DLL3 | 0.5 | 0.06 |  |  |  |  | 0.03 | 0.05 |
| 0.75 | 0.06 | 0.12 | 0.07 | 0.01 | 0.16 | 0.03 | 0.05 |

Note: Zero estimates are left blank.

Table 7. The model selected by the Het-QR method

|  | Estimated Regression Coefficients |
---|---|---|---
Gene |  |  | $\tau=0.25$ | $\tau=0.5$ | $\tau=0.75$ |  |
SDHA |  |  |  |  | 0.2 |  |
BMP2 |  |  | 0.35 | 0.31 | 0.34 |  |
CDKN2C |  |  |  | 0.13 | 0.25 |  |
DLL3 |  |  |  |  | 0.07 |  |
EGFR |  |  |  | -0.08 | -0.29 |  |
GRIA2 |  |  | 0.28 | 0.22 | 0.18 |  |
LTF |  |  |  |  | 0.07 |  |
OLIG2 |  |  | 0.14 | 0.30 | 0.38 |  |
PLAT |  |  | 0.20 | 0.21 | 0.25 |  |
SLC4A4 |  |  | -0.21 | -0.25 | -0.24 |  |
TAGLN |  |  |  |  | -0.20 |  |
TMEM100 |  |  | 0.20 | 0.17 |  |  |

##  6 Discussion

In this article, we have proposed a variable selection method that is able to
conduct joint variable selection and estimation for multiple quantiles
simultaneously. The joint selection/estimation allows one to harness the
strength shared among multiple quantiles and to achieve a model that is closer
to the truth. In particular, our approach is able to handle the heterogeneous
sparsity, under which a covariate contributes to some (but not all) of the
quantiles. By considering the heterogeneous sparsity, one can better dissect
the regression structure of the trait over the covariates, which in turn leads
to more accurate characterization of the underlying biological mechanism.

We have conducted a series of simulation studies to evaluate the performance
of our proposed approach as well as other approaches. Our simulation studies
show that the proposed method has superior performance to its peer methods. In
real data analysis, our method tends to yield a sparser model than the
compared methods. The benefit of achieving a sparse model is of great
importance to biological studies, because it helps biological investigators to
narrow down important candidate covariates (such as genes or proteins), so
that research efforts can be leveraged more efficiently. Our analysis
indicates that the regression coefficients at different quantiles can be quite
heterogeneous. We suggest that the interpretation of the results be guided by
biological knowledge and scientific insight, and that the variability be
examined by experimental studies. FAL and FAS were mainly designed to generate
interquantile shrinkage for quantile regression; when a smooth $\gamma(\tau)$
(with respect to $\tau$) is desired, these two methods are highly suitable and
are indeed the only available methods to achieve such a goal.

We have also provided theoretical proof for the proposed method under the
situation that $p$ can grow to infinity. Our exploratory experiments suggest
that Het-QR can be potentially applied to ‘$p>n$’, although theoretical work
is still needed to guide future experiments in this direction. Wang et al.
(2012) proposed a novel approach for studying asymptotics under the ‘$p>n$’
situation, and they focused on the penalties that can be written as the
difference of two convex functions. The group penalty considered herein does
not seem to fall into their framework. Further theoretical development is
merited. While we have imposed equal weights for multiple quantiles in this
paper, our method can be easily extended to accommodate different weights for
different quantiles. Properly chosen weights may lead to improved efficiency
of the estimated parameters (Zhao and Xiao, 2013).

## Acknowledgments

The authors are grateful to the AE and two reviewers for many helpful and
constructive comments. Dr. He’s research is supported by the Institutional
Support from the Fred Hutchinson Cancer Research Center. Dr. Kong’s research
is supported by Natural Sciences and Engineering Research Council of Canada.
The authors thank Dr. Huixia Judy Wang for generously providing the code for
FAL and FAS, which greatly expands the breadth of this manuscript. The authors
thank Dr. Li Hsu for helpful discussions.

The results shown here are in part based upon data generated by the TCGA
Research Network: http://cancergenome.nih.gov/.

Conflict of Interest: None declared.

## Appendix

### Transformation of the objective function

Our proof is in vein with the proof of Proposition 1 in Huang et al. [2009].
Consider the transformed objective function

$\displaystyle{\min}_{\theta,\xi}\sum_{m=1}^{M}\sum_{i=1}^{n}\,\rho_{m}(Y_{i}-U_{i}^{\rm
T}\theta_{m})+\lambda_{1}\sum_{j=1}^{p}\xi_{j}+\sum_{j=1}^{p}\xi_{j}^{-1}\left(\sum_{m=1}^{M}\omega_{mj}|\gamma_{mj}|\right)$
(5) $\displaystyle=$
$\displaystyle{\min}_{\theta}\left\\{\sum_{m=1}^{M}\sum_{i=1}^{n}\,\rho_{m}(Y_{i}-U_{i}^{\rm
T}\theta_{m})+\min_{\xi}\left\\{\lambda_{1}\sum_{j=1}^{p}\xi_{j}+\sum_{j=1}^{p}\xi_{j}^{-1}\left(\sum_{m=1}^{M}\omega_{mj}|\gamma_{mj}|\right)\right\\}\right\\}$

By Cauchy-Schwartz inequality, we have

$\lambda_{1}\sum_{j=1}^{p}\xi_{j}+\sum_{j=1}^{p}\xi_{j}^{-1}\left(\sum_{m=1}^{M}\omega_{mj}|\gamma_{mj}|\right)\geq
2\sqrt{\lambda}_{1}\sum_{j=1}^{p}\left(\sum_{m=1}^{M}\omega_{mj}|\gamma_{mj}|\right)^{\frac{1}{2}},$

then it follows that (5) is equivalent to

$\displaystyle\min_{\theta}\sum_{m=1}^{M}\sum_{i=1}^{n}\,\rho_{m}(Y_{i}-U_{i}^{\rm
T}\theta_{m})+2\sqrt{\lambda}_{1}\sum_{j=1}^{p}\left(\sum_{m=1}^{M}\omega_{mj}|\gamma_{mj}|\right)^{\frac{1}{2}}.$
(6)

Now, let $2\sqrt{\lambda}_{1}=n\lambda_{n}$, then (6) is identical to the
original objective function (2). $\Box$

### Derivation of the primal and dual problem

Let $\lambda_{mj}=\xi_{j}^{-1}\omega_{mj}$, then in step 2 of Section 2.2, we
aim to solve

${\min}_{\gamma,\gamma_{0}}\left\\{\sum_{m=1}^{M}\sum_{i=1}^{n}\,\rho_{m}(Y_{i}-Z_{i}^{\rm
T}\gamma_{m}-\gamma_{m0})+\sum_{j=1}^{p}\sum_{m=1}^{M}\lambda_{mj}|\gamma_{mj}|\right\\}.$

Let $e_{n}$ denote the unit vector of length $n$, and
$\boldsymbol{\lambda}_{m}$ the vector of $\lambda_{mj}(j=1,\ldots,p)$. With
slight abuse of notation, let $Y$ be the $n\times 1$ vector consisting of
$Y_{i}$, and $Z$ the $n\times p$ matrix consisting of $Z_{i}$. The above
objective function is equivalent to

$\min_{u_{m},v_{m},\gamma_{m},s_{m},t_{m}}\sum_{m=1}^{M}\tau_{m}e_{n}^{T}u_{m}+(1-\tau_{m})e_{n}^{T}v_{m}+\boldsymbol{\lambda}_{m}^{T}s_{m}+\boldsymbol{\lambda}_{m}^{T}t_{m},$

subject to $u_{m}-v_{m}=Y-Z\gamma_{m}-\gamma_{m0}e_{n}$,
$s_{m}-t_{m}=\gamma_{m}$, $u_{m}\geq 0,v_{m}\geq 0,s_{m}\geq 0,$ and
$t_{m}\geq 0$.

Let ${\bf 0}_{r}$ be the zero vector of length $r$,
$\lambda^{*}=(\boldsymbol{\lambda}_{1}^{T},\ldots,\boldsymbol{\lambda}_{M}^{T})^{T}$,
$s^{*}=(s_{1}^{T},\ldots,s_{M}^{T})^{T}$,
$t^{*}=(t_{1}^{T},\ldots,t_{M}^{T})$,
$u^{*}=(u_{1}^{T},\ldots,u_{M}^{T})^{T}$,
$v^{*}=(v_{1}^{T},\ldots,v_{M}^{T})^{T}$. Let $Y_{(M)}$ denote the vector in
which $Y$ is stacked by $M$ times. Let

$c=({\bf 0}_{Mp}^{T},{\bf
0}_{M}^{T},{\lambda^{*}}^{T},{\lambda^{*}}^{T},\tau_{1}e_{n}^{T},\ldots,\tau_{M}e_{n}^{T},(1-\tau_{1})e_{n}^{T},\ldots,(1-\tau_{M})e_{n}^{T})^{T},$

$x=(\gamma^{T},\gamma_{0}^{T},{s^{*}}^{T},{t^{*}}^{T},{u^{*}}^{T},{v^{*}}^{T})^{T},$
and $b=(Y_{(M)}^{T},{\bf 0}_{Mp}^{T})^{T}$. Then, the linear program primal of
the above objective function can be written as

$\min_{x}c^{T}x,$

subject to $Ax=b$ and
$({s^{*}}^{T},{t^{*}}^{T},{u^{*}}^{T},{v^{*}}^{T})^{T}\geq 0$, where $A$ is
defined as follows. $A$ is a matrix consisting of two rows of blocks. The
first row of $A$ consists of 6 blocks, $A_{11}=I_{M}\otimes Z$,
$A_{12}=I_{M}\otimes e_{n}$, $A_{13}=[{\bf 0}]_{nM\times Mp}$, $A_{14}=[{\bf
0}]_{nM\times Mp}$, $A_{15}=I_{nM}$, and $A_{16}=-I_{nM}$. The second row of
$A$ consists of the following 6 blocks, $A_{21}=I_{Mp}$, $A_{22}=[0]_{Mp\times
M}$, $A_{23}=-I_{Mp}$, $A_{24}=I_{Mp}$, $A_{25}=[{\bf 0}]_{nM\times nM}$, and
$A_{26}=[{\bf 0}]_{nM\times nM}$.

Then using standard linear program arguments, we obtain the dual as

$\max_{\tilde{d}}b^{T}\tilde{d},$

subject to $\tilde{Z}^{T}\tilde{d}=S_{1}+S_{2}$ and
$\tilde{d}\in[0,1]^{nM+Mp},$ where

$S_{1}=\left((1-\tau_{1})e_{n}^{T}Z,\ldots,(1-\tau_{M})e_{n}^{T}Z,n(1-\tau_{1}),\ldots,n(1-\tau_{M})\right)^{T},$

$S_{2}=1/2\times(R,[{\bf 0}]_{Mp\times M})^{T}e_{Mp}$,
$R=diag((2\boldsymbol{\lambda}^{T},\ldots,2\boldsymbol{\lambda}^{T}_{M})^{T})$,
and $\tilde{Z}$ is defined as follows. $\tilde{Z}$ consists of two rows of
blocks. The first row of $\tilde{Z}$ includes two blocks,
$\tilde{Z}_{11}=I_{M}\otimes Z$ and $\tilde{Z}_{12}=I_{M}\otimes e_{n}$. The
second row includes two blocks, $\tilde{Z}_{21}=R$ and $\tilde{Z}_{22}=[{\bf
0}]_{Mp\times M}$.

### Computation time

For $p=100$ under the auto-regressive structure, i.e., Table 2, we calculate
the summary statistics of the CPU time (seconds). The average time (and the
standard error of the sample mean) for QR, QR-LASSO, QR-aLASSO, FAL, FAS, Het-
QR is 1.6(0.003), 10.2(0.017), 10.0(0.030), 867.0(8.059), 455.3(2.641),
90.3(0.391), respectively. For $p=600$ under the auto-regressive structure,
the CPU time for QR-LASSO is 372.6(1.453); the time for QR-aLASSO and Het-QR
is 362.0(1.729) and 3717.5(36.401), respectively (excluding the time for
calculating penalty weights).

### Proof of the theorems

We now give the proof of the theorems. We note that throughout the proof, the
upper letter $C$ in different formulas stands for different constants.

Recall that the definitions of $\mathcal{N}$, $I$, $II$, $s_{n}$ and $I_{m}$
have been given in the main text. Define the index set $J=\\{1\leq j\leq
p_{n}:\text{there exists}\,1\leq m\leq M\,\text{such
that}\,\gamma_{mj}^{*}\neq 0\\}$ with the cardinality $|J|=k_{n}$. For every
fixed $j\in J$, define the index set $M_{j}=\\{1\leq m\leq
M:\gamma_{mj}^{*}\neq 0\\}$. Clearly, the oracle index set $I=\\{(m,j):m\in
M_{j},\,j\in J\\}$. For every fixed $1\leq m\leq M$, define $II_{m}=\\{1\leq
j\leq p_{n}:\gamma_{mj}^{*}=0\\}$.

We need to define some notations for our proof. Let the vectors $\gamma_{mI}$,
$\gamma_{mI}^{*}$ and $\hat{\gamma}_{mI}$ be the subvectors of $\gamma_{m}$,
$\gamma_{m}^{*}$ and $\hat{\gamma}_{m}$ corresponding to the index set
$I_{m}$, respectively. Define the subvectors of $\gamma$, $\gamma^{*}$ and
$\hat{\gamma}$ corresponding to the oracle index set $I$ as
$\gamma_{I}=(\gamma_{1I}^{\rm T},\cdots,\gamma_{MI}^{\rm T})^{\rm T}$,
$\gamma_{I}^{*}=(\gamma_{1I}^{*\rm T},\cdots,\gamma_{MI}^{*\rm T})^{\rm T}$
and $\hat{\gamma}_{I}=(\hat{\gamma}_{1I}^{\rm T},\cdots,\hat{\gamma}_{MI}^{\rm
T})^{\rm T}$. Let $\theta_{mI}=(\gamma_{m0},\gamma_{mI}^{\rm T})^{\rm T}$,
$\theta_{mI}^{*}=(\gamma_{m0}^{*},\gamma_{mI}^{*\rm T})^{\rm T}$ and
$\hat{\theta}_{mI}=(\hat{\gamma}_{m0},\hat{\gamma}_{mI}^{\rm T})^{\rm T}$.
Define the vector $\theta_{I}$ as the subvector of the parameter vector
$\theta$ corresponding to the oracle index set $I$. Recall that the vectors
$\theta_{I}^{*}$ and $\hat{\theta}_{I}$ are the subvectors of the vectors
$\theta^{*}$ and $\hat{\theta}$ corresponding to the oracle index set $I$.
Clearly, $\theta_{I}=(\theta_{1I}^{\rm T},\cdots,\theta_{MI}^{\rm T})^{\rm
T}$, $\theta_{I}^{*}=(\theta_{1I}^{*\rm T},\cdots,\theta_{MI}^{*\rm T})^{\rm
T}$ and $\hat{\theta}_{I}=(\hat{\theta}_{1I}^{\rm
T},\cdots,\hat{\theta}_{MI}^{\rm T})^{\rm T}$.

Similarly, let the vectors $\gamma_{mII}$, $\gamma_{mII}^{*}$ and
$\hat{\gamma}_{mII}$ be the subvectors of $\gamma_{m}$, $\gamma_{m}^{*}$ and
$\hat{\gamma}_{m}$ corresponding to the index set $II_{m}$, respectively.
Define the subvectors of $\gamma$, $\gamma^{*}$ and $\hat{\gamma}$
corresponding to the index set $II$ as $\gamma_{II}=(\gamma_{1II}^{\rm
T},\cdots,\gamma_{MII}^{\rm T})^{\rm T}$, $\gamma_{II}^{*}=(\gamma_{1II}^{*\rm
T},\cdots,\gamma_{MII}^{*\rm T})^{\rm T}$ and
$\hat{\gamma}_{II}=(\hat{\gamma}_{1II}^{\rm T},\cdots,\hat{\gamma}_{MII}^{\rm
T})^{\rm T}$. Define the vector $\theta_{II}$ as the subvector of the
parameter vector $\theta$ corresponding to the index set $II$. Recall that the
vectors $\theta_{II}^{*}$ and $\hat{\theta}_{II}$ are the subvectors of the
vectors $\theta^{*}$ and $\hat{\theta}$ corresponding to the index set $II$,
respectively. Clearly, $\theta_{II}=\gamma_{II}$,
$\theta_{II}^{*}=\gamma_{II}^{*}=0$ and $\hat{\theta}_{II}=\hat{\gamma}_{II}$.

For convenience, write $\gamma=(\gamma_{I}^{\rm T},\gamma_{II}^{\rm T})^{\rm
T}$, $\gamma^{*}=(\gamma_{I}^{*\rm T},\gamma_{II}^{*\rm T})^{\rm T}$ and
$\hat{\gamma}=(\hat{\gamma}_{I}^{\rm T},\hat{\gamma}_{II}^{\rm T})^{\rm T}$;
and $\theta=(\theta_{I}^{\rm T},\theta_{II}^{\rm T})^{\rm T}$,
$\theta^{*}=(\theta_{I}^{*\rm T},\theta_{II}^{*\rm T})^{\rm T}$ and
$\hat{\theta}=(\hat{\theta}_{I}^{\rm T},\hat{\theta}_{II}^{\rm T})^{\rm T}$.

We first give a lemma related to the loss function $Q_{n}(\theta)$. The lemma
plays an important role in the proof of our theorems.

Lemma 1. Under conditions (L1)-(L4), we have

$\displaystyle Q_{n}(\theta)=Q_{n}(\theta^{*})+\mathbf{A}_{n}^{\rm
T}(\theta-\theta^{*})+\frac{1}{2}(\theta-\theta^{*})^{\rm
T}\mathbf{B}_{n}(\theta-\theta^{*})+R_{n}(\theta),$

where
$\sup_{1\leq\,m\leq\,M}\,\sup_{\|\theta_{m}-\theta_{m}^{*}\|_{2}\,\leq\,\eta_{n}}\,|R_{n}(\theta)|=O_{p}(n^{\frac{1}{2}}p^{\frac{3}{4}}\eta_{n}^{\frac{3}{2}})+O_{p}(np^{\frac{1}{2}}\eta_{n}^{3})$,
and

$\displaystyle\mathbf{A}_{n}$ $\displaystyle=(A_{1}^{\rm T},\cdots,A_{M}^{\rm
T})^{\rm T}\quad\text{with}\quad
A_{m}=\sum_{i=1}^{n}\,\bigl{(}I(Y_{i}<U_{i}^{\rm
T}\theta_{m}^{*})-\tau_{m}\bigr{)}\,U_{i},$ $\displaystyle\mathbf{B}_{n}$
$\displaystyle=Diag(B_{11},B_{22},\cdots,B_{MM})\quad\text{with}\quad
B_{mm}=\sum_{i=1}^{n}\,f(U_{i}^{\rm T}\theta_{m}^{*}|\,Z_{i})\,U_{i}U_{i}^{\rm
T}.$ (7)

Proof of Lemma 1. Let $\psi_{m}(u)$ be a sub-derivative of the quantile
function $\rho_{m}(u)$, then $\psi_{m}(u)=\tau_{m}-I(u<0)+l_{m}\,I(u=0)$ with
$l_{m}\in[-1,0]$. Let $T_{n}=Q_{n}(\theta)-Q_{n}(\theta^{*})$. Then, there
exists an $l_{mi}\in[-1,0]$ for every $1\leq m\leq M$ and $1\leq i\leq n$ such
that

$\displaystyle T_{n}$
$\displaystyle=-\sum_{m=1}^{M}\sum_{i=1}^{n}\,\psi_{m}(Y_{i}-U_{i}^{\rm
T}\,\bar{\theta}_{m})\,U_{i}^{\rm T}(\theta_{m}-\theta_{m}^{*})$
$\displaystyle=-\sum_{m=1}^{M}\sum_{i=1}^{n}\,\bigl{(}\tau_{m}-I(Y_{i}<U_{i}^{\rm
T}\,\bar{\theta}_{m})+l_{mi}\,I(Y_{i}=U_{i}^{\rm
T}\,\bar{\theta}_{m})\bigr{)}\,U_{i}^{\rm T}(\theta_{m}-\theta_{m}^{*})$
$\displaystyle=\sum_{m=1}^{M}\sum_{i=1}^{n}\,\bigl{(}I(Y_{i}<U_{i}^{\rm
T}\,\theta_{m}^{*})-\tau_{m}\bigr{)}\,U_{i}^{\rm
T}(\theta_{m}-\theta_{m}^{*})$
$\displaystyle\quad-\sum_{m=1}^{M}\sum_{i=1}^{n}\,\bigl{(}I(Y_{i}<U_{i}^{\rm
T}\,\theta_{m}^{*})-I(Y_{i}<U_{i}^{\rm T}\,\bar{\theta}_{m})\bigr{)}U_{i}^{\rm
T}(\theta_{m}-\theta_{m}^{*})$
$\displaystyle\quad-\sum_{m=1}^{M}\sum_{i=1}^{n}\,l_{mi}\,I(Y_{i}=U_{i}^{\rm
T}\,\bar{\theta}_{m})U_{i}^{\rm T}(\theta_{m}-\theta_{m}^{*})$
$\displaystyle\equiv T_{n1}-T_{n2}-T_{n3}.$ (8)

where $\bar{\theta}_{m}$ is on the linear segment between $\theta_{m}$ and
$\theta_{m}^{*}$, and may be written as
$\bar{\theta}_{m}=\theta_{m}^{*}+\eta_{m}\,(\theta_{m}-\theta_{m}^{*})$ with
$\eta_{m}\in(0,1)$. For $T_{n3}$, note that $Y_{i}$ has a continuous
conditional distribution given $Z_{i}$, hence almost surely
$I(Y_{i}=U_{i}^{\rm T}\,\bar{\theta}_{m})=0$ for all $i=1,\cdots,n$ and
$m=1,\cdots,M$, thus $T_{n3}=0$ almost surely. Subsequently, we can write

$T_{n}=\bigl{(}T_{n1}-\mathbf{E}(T_{n1})\bigr{)}-\bigl{(}T_{n2}-\mathbf{E}(T_{n2})\bigr{)}+\mathbf{E}(T_{n}),$

where $\mathbf{E}(T_{n})=\mathbf{E}(T_{n1})+\mathbf{E}(T_{n2})$, and
$\mathbf{E}$ denotes the conditional expectation given $Z$. Note that
$\mathbf{E}(T_{n1})=0$ because $\mathbf{E}(I(Y_{i}<U_{i}^{\rm
T}\,\theta_{m}^{*})-\tau_{m}|Z_{i})=F(U_{i}^{\rm
T}\,\theta_{m}^{*}|Z_{i})-\tau_{m}=0$ from $\eqref{0-mod-1}$. Rename
$\bigl{(}T_{n2}-\mathbf{E}(T_{n2})\bigr{)}$ as $R_{n2}$ and
$\mathbf{E}(T_{n})$ as $T_{n4}$, then we have

$\displaystyle T_{n}$ $\displaystyle=T_{n1}-R_{n2}+T_{n4}.$ (9)

For $R_{n2}$ (recall $T_{n2}$ in $\eqref{1-tn}$), let
$\zeta_{mi}(t)=I(Y_{i}-U_{i}^{\rm T}\,\theta_{m}^{*}<0)-I(Y_{i}-U_{i}^{\rm
T}\,\theta_{m}^{*}<U_{i}^{\rm T}\,t)$, then

$\displaystyle R_{n2}$
$\displaystyle=\sum_{m=1}^{M}\biggl{(}\sum_{i=1}^{n}\,\\{\zeta_{mi}(\eta_{m}(\theta_{m}-\theta_{m}^{*}))-\mathbf{E}(\zeta_{mi}(\eta_{m}(\theta_{m}-\theta_{m}^{*})))\\}\,U_{i}\biggr{)}^{\rm
T}(\theta_{m}-\theta_{m}^{*})$ $\displaystyle=\sum_{m=1}^{M}\phi_{n}(m)^{\rm
T}(\theta_{m}-\theta_{m}^{*})$ (10)

where
$\phi_{n}(m)=\sum_{i=1}^{n}\,\\{\zeta_{mi}(\eta_{m}(\theta_{m}-\theta_{m}^{*}))-\mathbf{E}(\zeta_{mi}(\eta_{m}(\theta_{m}-\theta_{m}^{*})))\\}\,U_{i}.$

Note that $|\zeta_{mi}(t)|=|I(Y_{i}-U_{i}^{\rm
T}\,\theta_{m}^{*}<0)-I(Y_{i}-U_{i}^{\rm T}\,\theta_{m}^{*}<U_{i}^{\rm
T}\,t)|\leq I(|Y_{i}-U_{i}^{\rm T}\,\theta_{m}^{*}|\leq|U_{i}^{\rm T}\,t|)$,
and that $f(t|Z_{i})$ is bounded under condition (L1) and $\|U_{i}\|_{2}\leq
Cp^{\frac{1}{2}}$ under condition (L2). Hence, making use of independence,
under conditions (L1)and (L2), for all $1\leq m\leq M$ and
$\|\theta_{m}-\theta_{m}^{*}\|_{2}\leq\eta_{n}$, we can see

$\displaystyle\mathbf{E}\|\phi_{n}(m)\|_{2}^{2}$
$\displaystyle=\mathbf{E}\biggl{\|}\sum_{i=1}^{n}\,\\{\zeta_{mi}(\eta_{m}(\theta_{m}-\theta_{m}^{*}))-\mathbf{E}(\zeta_{mi}(\eta_{m}(\theta_{m}-\theta_{m}^{*})))\\}\,U_{i}\biggr{\|}_{2}^{2}$
$\displaystyle=\sum_{i=1}^{n}\,\mathbf{E}\\{\zeta_{mi}(\eta_{m}(\theta_{m}-\theta_{m}^{*}))-\mathbf{E}(\zeta_{mi}(\eta_{m}(\theta_{m}-\theta_{m}^{*})))\\}^{2}\|U_{i}\|_{2}^{2}$
$\displaystyle\leq
Cnp\mathbf{E}\zeta_{mi}(\eta_{m}(\theta_{m}-\theta_{m}^{*}))^{2}\leq
Cnp\,P\bigl{(}|Y_{i}-U_{i}^{\rm T}\,\theta_{m}^{*}|\leq|U_{i}^{\rm
T}(\theta_{m}-\theta_{m}^{*})|\bigr{)}$ $\displaystyle\leq
CnpP\bigl{(}|Y_{i}-U_{i}^{\rm T}\,\theta_{m}^{*}|\leq
Cp^{\frac{1}{2}}\|\theta_{m}-\theta_{m}^{*}\|_{2})$

Note that

$\displaystyle P\bigl{(}|Y_{i}-U_{i}^{\rm T}\,\theta_{m}^{*}|\leq
Cp^{\frac{1}{2}}\|\theta_{m}-\theta_{m}^{*}\|_{2})$
$\displaystyle=F(U_{i}^{\rm
T}\,\theta_{m}^{*}+Cp^{\frac{1}{2}}\|\theta_{m}-\theta_{m}^{*}\|_{2})-F(U_{i}^{\rm
T}\,\theta_{m}^{*}-Cp^{\frac{1}{2}}\|\theta_{m}-\theta_{m}^{*}\|_{2})$
$\displaystyle\leq\,|f(\xi_{mi}|Z_{i})|p^{\frac{1}{2}}\|\theta_{m}-\theta_{m}^{*}\|_{2}\leq
Cp^{\frac{1}{2}}\|\theta_{m}-\theta_{m}^{*}\|_{2},$

where $\xi_{mi}$ is between $U_{i}^{\rm
T}\,\theta_{m}^{*}+Cp^{\frac{1}{2}}\|\theta_{m}-\theta_{m}^{*}\|_{2}$ and
$U_{i}^{\rm
T}\,\theta_{m}^{*}-Cp^{\frac{1}{2}}\|\theta_{m}-\theta_{m}^{*}\|_{2}$. Thus,
$\mathbf{E}\|\phi_{n}(m)\|_{2}^{2}\leq Cnp^{\frac{3}{2}}\eta_{n}.$ By
Chebyshev’s inequality, we get
$\sup_{1\leq\,m\leq\,M}\sup_{\|\theta_{m}-\theta_{m}^{*}\|_{2}\,\leq\,\eta_{n}}\,\|\phi_{n}(m)\|_{2}=O_{p}(n^{\frac{1}{2}}p^{\frac{3}{4}}\eta_{n}^{\frac{1}{2}})$.
Together with $\eqref{1-rn2-0}$, by Cauchy-Schwartz’s inequality, we get

$\displaystyle\sup_{1\leq\,m\leq\,M}\,\sup_{\|\theta_{m}-\theta_{m}^{*}\|_{2}\,\leq\,\eta_{n}}\,|R_{n2}|=O_{p}(n^{\frac{1}{2}}p^{\frac{3}{4}}{\eta_{n}}^{\frac{1}{2}})O(\eta_{n})=O_{p}(n^{\frac{1}{2}}p^{\frac{3}{4}}{\eta_{n}}^{\frac{3}{2}}).$
(11)

For $T_{n4}$, the third term in $\eqref{1-tn-1}$, write
$\mathbf{E}(T_{n})=e_{n}(\theta)-e_{n}(\theta^{*})$, where

$\displaystyle
e_{n}(\theta)\equiv\sum_{m=1}^{M}\sum_{i=1}^{n}\,\mathbf{E}\,\rho_{m}\bigl{(}Y_{i}-U_{i}^{\rm
T}\theta_{m}\bigr{)}.$

$\mathbf{E}\,\rho_{m}\bigl{(}Y_{i}-U_{i}^{\rm T}\theta_{m}\bigr{)}$ is second
order differentiable with respect to $\theta_{m}$ under condition (L1), with
gradient $G_{mi}(\theta)=-\mathbf{E}\,\bigl{\\{}(\tau_{m}-I(Y_{i}<U_{i}^{\rm
T}\theta_{m}))U_{i}\bigr{\\}}=\bigl{(}F(U_{i}^{\rm
T}\theta_{m}\bigr{|}Z_{i})-\tau_{m}\bigr{)}\,U_{i}$, and Hessian matrix
$H_{mi}(\theta)=f(\bigl{(}U_{i}^{\rm
T}\theta_{m})\bigr{|}Z_{i}\bigr{)}\,U_{i}U_{i}^{\rm T}$.

Let $G(\theta)$ and $H(\theta)$ be gradient and Hessian matrix of
$e_{n}(\theta)$, then $G(\theta)=\sum_{m=1}^{M}\sum_{i=1}^{n}\,G_{mi}(\theta)$
and
$H(\theta)=\sum_{m=1}^{M}\sum_{i=1}^{n}\,H_{mi}(\theta)=\sum_{m=1}^{M}\sum_{i=1}^{n}\,f(\bigl{(}U_{i}^{\rm
T}\theta_{m})\bigr{|}Z_{i}\bigr{)}\,U_{i}U_{i}^{\rm T}$. It is easy to see
that $G(\theta^{*})=0$ by $F(U_{i}^{\rm T}\theta_{m}^{*}|Z_{i})=\tau_{m}$ in
$\eqref{0-mod-1}$. By Taylor expansion of
$T_{n4}=\mathbf{E}(T_{n})=e_{n}(\theta)-e_{n}(\theta^{*})$ at $\theta^{*}$, we
have

$\displaystyle T_{n4}$ $\displaystyle=\frac{1}{2}(\theta-\theta^{*})^{\rm
T}\,H\bigl{(}\theta^{*}+\xi(\theta-\theta^{*})\bigr{)}(\theta-\theta^{*})$
$\displaystyle=\frac{1}{2}\sum_{m=1}^{M}\sum_{i=1}^{n}\,f(\zeta_{mi}|Z_{i})\,(\theta_{m}-\theta_{m}^{*})^{\rm
T}\,U_{i}U_{i}^{\rm T}(\theta_{m}-\theta_{m}^{*}),$
$\displaystyle=\frac{1}{2}\sum_{m=1}^{M}\sum_{i=1}^{n}\,\bigl{\\{}f(\zeta_{mi}|Z_{i})-f(U_{i}^{\rm
T}\theta_{m}^{*}|Z_{i})\bigr{\\}}\,(\theta_{m}-\theta_{m}^{*})^{\rm
T}\,U_{i}U_{i}^{\rm T}(\theta_{m}-\theta_{m}^{*})$
$\displaystyle\quad+\frac{1}{2}\sum_{m=1}^{M}\sum_{i=1}^{n}\,f(U_{i}^{\rm
T}\theta_{m}^{*}|Z_{i})(\theta_{m}-\theta_{m}^{*})^{\rm T}U_{i}U_{i}^{\rm
T}(\theta_{m}-\theta_{m}^{*})$ $\displaystyle\equiv R_{n4}+T_{n41}.$ (12)

where $\xi\in(0,1)$ and $\zeta_{mi}$ is between $U_{i}^{\rm T}\theta_{m}$ and
$U_{i}^{\rm T}\theta_{m}^{*}$. Trivially,

$\displaystyle|R_{n4}|\leq C\sum_{m=1}^{M}\sup_{1\leq i\leq
n}\bigl{|}f(\zeta_{mi}|Z_{i})-f(U_{i}^{\rm
T}\theta_{m}^{*}|Z_{i})\bigr{|}(\theta_{m}-\theta_{m}^{*})^{\rm
T}\sum_{i=1}^{n}\,\,U_{i}U_{i}^{\rm T}(\theta_{m}-\theta_{m}^{*}).$ (13)

Note that $f^{\prime}(t|Z_{i})$ is bounded under condition (L1),
$\|U_{i}\|_{2}\leq Cp^{\frac{1}{2}}$ under condition (L2), and
$\lambda_{max}\bigl{(}n^{-1}\sum_{i=1}^{n}\,U_{i}U_{i}^{\rm T}\bigr{)}\leq C$
under condition (L3). Hence, for all $1\leq i\leq n$, $1\leq m\leq M$, and
$\|\theta_{m}-\theta_{m}^{*}\|_{2}\leq\eta_{n}$, we have
$|f(\zeta_{mi}|Z_{i})-f(U_{i}^{\rm T}\theta_{m}^{*}|Z_{i})|\leq C|U_{i}^{\rm
T}(\theta_{m}-\theta_{m}^{*})|\leq\,C\|U_{i}\|_{2}\|\theta_{m}-\theta_{m}^{*}\|_{2}\leq\,Cp^{\frac{1}{2}}\eta_{n}$,
and $(\theta_{m}-\theta_{m}^{*})^{\rm T}\sum_{i=1}^{n}\,\,U_{i}U_{i}^{\rm
T}(\theta_{m}-\theta_{m}^{*})\leq
n\lambda_{max}\bigl{(}n^{-1}\sum_{i=1}^{n}\,U_{i}U_{i}^{\rm
T}\bigr{)}\|\theta_{m}-\theta_{m}^{*}\|_{2}^{2}\leq Cn\eta_{n}^{2}$. Hence,
from $\eqref{1-rn4}$, we get
$\sup_{1\leq\,m\leq\,M}\,\sup_{\|\theta_{m}-\theta_{m}^{*}\|_{2}\,\leq\,\eta_{n}}\,|R_{n4}|\leq
Cnp^{\frac{1}{2}}\eta_{n}^{3}$. Together with $\eqref{1-tn-1}$,
$\eqref{1-rn2}$ and $\eqref{1-tn4}$, we obtain
$T_{n}=T_{n1}+T_{n41}+R_{n}(\theta)$, where
$\sup_{1\leq\,m\leq\,M}\,\sup_{\|\theta_{m}-\theta_{m}^{*}\|_{2}\,\leq\,\eta_{n}}\,|R_{n}(\theta)|=O_{p}(n^{\frac{1}{2}}p^{\frac{3}{4}}{\eta_{n}}^{\frac{3}{2}})+O_{p}(np^{\frac{1}{2}}\eta_{n}^{3})$,
and

$\displaystyle T_{n1}$
$\displaystyle=\sum_{m=1}^{M}\sum_{i=1}^{n}\,\bigl{(}I(Y_{i}<U_{i}^{\rm
T}\theta_{m}^{*})-\tau_{m}\bigr{)}\,U_{i}^{\rm
T}(\theta_{m}-\theta_{m}^{*})=\mathbf{A}_{n}^{\rm T}(\theta-\theta^{*}),$
$\displaystyle T_{n41}$
$\displaystyle=\frac{1}{2}\sum_{m=1}^{M}\sum_{i=1}^{n}\,f(U_{i}^{\rm
T}\theta_{m}^{*}|Z_{i})(\theta_{m}-\theta_{m}^{*})^{\rm T}U_{i}U_{i}^{\rm
T}(\theta_{m}-\theta_{m}^{*})=\frac{1}{2}(\theta-\theta^{*})^{\rm
T}\mathbf{B}_{n}(\theta-\theta^{*}).$

This completes the proof of the lemma. $\Box$

Proof of Theorem 1. Recall the definition of
$L_{n}(\theta)=Q_{n}(\theta)+P_{n}(\gamma)$ in $\eqref{0-ln}$. Let
$\theta-\theta^{*}=\nu_{n}\,u$, where $\nu_{n}>0$, $u\in\mathbb{R}^{M(p+1)}$
and $\|u\|_{2}=1$. It is easy to see that $\|\theta-\theta^{*}\|_{2}=\nu_{n}$.
Based on the continuity of $L_{n}$, if we can prove that in probability

$\displaystyle\inf_{\|u\|_{2}=1}\,L_{n}(\theta^{*}+\nu_{n}u)>L_{n}(\theta^{*}),$
(14)

then the minimal value point of $L_{n}(\theta^{*}+\nu_{n}u)$ on
$\\{u:\|u\|_{2}\leq 1\\}$ exists and lies in the unit ball $\\{u:\|u\|_{2}\leq
1\\}$ in probability. We will prove that $\eqref{2-1}$ holds.

For $Q_{n}(\theta)$, because
$\|\theta_{m}-\theta_{m}^{*}\|_{2}\,\leq\,\|\theta-\theta^{*}\|_{2}$ for all
$m=1,2,\cdots,M$, by Lemma 1 with $\eta_{n}=\nu_{n}$ under conditions (L1) to
(L4), we get

$\displaystyle q_{n}(\theta)\equiv
Q_{n}(\theta)-Q_{n}(\theta^{*})=\mathbf{A}_{n}^{\rm
T}(\theta-\theta^{*})+\frac{1}{2}(\theta-\theta^{*})^{\rm
T}\mathbf{B}_{n}(\theta-\theta^{*})+R_{n}(\theta),$ (15)

where
$\sup_{\|\theta-\theta^{*}\|_{2}\leq\nu_{n}}\,|R_{n}(\theta)|=O_{p}(n^{\frac{1}{2}}p^{\frac{3}{4}}{\nu_{n}}^{\frac{3}{2}})+O_{p}(np^{\frac{1}{2}}\nu_{n}^{3})$.

For $P_{n}(\gamma)$, let $p_{n}(\gamma)=P_{n}(\gamma)-P_{n}(\gamma^{*})$.
Define $p_{1n}(\gamma)=p_{n}(\gamma_{I},0)$, that is,

$\displaystyle p_{1n}(\gamma)=n\lambda_{n}\,\sum_{j\in J}\,\bigl{(}\sum_{m\in
M_{j}}\,\omega_{mj}|\gamma_{mj}|\bigr{)}^{\frac{1}{2}}-n\lambda_{n}\,\sum_{j\in
J}\,\bigl{(}\sum_{m\in
M_{j}}\,\omega_{mj}|\gamma_{mj}^{*}|\bigr{)}^{\frac{1}{2}}.$ (16)

Clearly, $p_{1n}(\gamma)\leq p_{n}(\gamma)$ and
$p_{1n}(\gamma^{*})=p_{n}(\gamma^{*})=0$.

Define $l_{n}(\theta)=q_{n}(\theta)+p_{n}(\gamma)$ and
$l_{1n}(\theta)=q_{n}(\theta)+p_{1n}(\gamma)$, both of which are continuous.
Clearly, $l_{1n}(\theta)\leq l_{n}(\theta)$ and
$l_{1n}(\theta^{*})=l_{n}(\theta^{*})=0$. Note that $\eqref{2-1}$ is
equivalent to that in probability

$\displaystyle\inf_{\|u\|_{2}=1}\,l_{1n}(\theta^{*}+\nu_{n}u)>0.$ (17)

Note that

$\displaystyle|p_{1n}(\gamma)|$
$\displaystyle=n\lambda_{n}\,\bigl{|}\sum_{j\in J}\,\bigl{(}\sum_{m\in
M_{j}}\,\omega_{mj}|\gamma_{mj}|\bigr{)}^{\frac{1}{2}}-\sum_{j\in
J}\,\bigl{(}\sum_{m\in
M_{j}}\,\omega_{mj}|\gamma_{mj}^{*}|\bigr{)}^{\frac{1}{2}}\bigr{|}$
$\displaystyle\leq n\lambda_{n}\,\sum_{j\in J}\,\bigl{|}\bigl{(}\sum_{m\in
M_{j}}\,\omega_{mj}|\gamma_{mj}|\bigr{)}^{\frac{1}{2}}-\bigl{(}\sum_{m\in
M_{j}}\,\omega_{mj}|\gamma_{mj}^{*}|\bigr{)}^{\frac{1}{2}}\bigr{|}\leq
n\lambda_{n}\,\sum_{j\in J}\,\bigl{(}\sum_{m\in
M_{j}}\,\omega_{mj}|\gamma_{mj}-\gamma_{mj}^{*}|\bigr{)}^{\frac{1}{2}},$

where the last inequality follows from the fact that
$|\sqrt{|x|}-\sqrt{|y|}|\leq\sqrt{|x-y|}$. Note that
$\gamma_{mj}-\gamma_{mj}^{*}=\nu_{n}\,u_{mj}$ where $u_{mj}$ is a component of
$u$. Hence,

$\displaystyle|p_{1n}(\gamma)|$ $\displaystyle\leq
n\lambda_{n}\,d_{nI}^{\frac{1}{2}}\nu_{n}^{\frac{1}{2}}\,\sum_{j\in
J}\,\bigl{(}\sum_{m\in M_{j}}\,|u_{mj}|\bigr{)}^{\frac{1}{2}}\leq
n\lambda_{n}\,d_{nI}^{\frac{1}{2}}\nu_{n}^{\frac{1}{2}}k_{n}^{\frac{1}{2}}\bigl{(}\sum_{(m,j)\in
I}\,|u_{mj}|\bigr{)}^{\frac{1}{2}}$ $\displaystyle\leq
n\lambda_{n}d_{nI}^{\frac{1}{2}}\nu_{n}^{\frac{1}{2}}k_{n}^{\frac{1}{2}}s_{n}^{\frac{1}{4}}\|u\|_{2}^{\frac{1}{2}}=n\lambda_{n}d_{nI}^{\frac{1}{2}}\nu_{n}^{\frac{1}{2}}s_{n}^{\frac{3}{4}}\|u\|_{2}^{\frac{1}{2}},$
(18)

where $k_{n}=|J|\leq s_{n}$. By $\eqref{2-qn}$ and the above, we have

$\displaystyle
l_{1n}(\theta^{*}+\nu_{n}u)=q_{n}(\theta^{*}+\nu_{n}u)+p_{1n}(\theta^{*}+\nu_{n}u)=\nu_{n}\mathbf{A}_{n}^{\rm
T}u+\frac{1}{2}\nu_{n}^{2}u^{\rm T}\mathbf{B}_{n}u+R_{L}(u),$ (19)

where
$\sup_{\|u\|_{2}\,=\,1}\,|R_{L}(u)|=O_{p}(n^{\frac{1}{2}}p^{\frac{3}{4}}{\nu_{n}}^{\frac{3}{2}})+O_{p}(np^{\frac{1}{2}}\nu_{n}^{3})+O(n\lambda_{n}\,d_{nI}^{\frac{1}{2}}\nu_{n}^{\frac{1}{2}}s_{n}^{\frac{3}{4}})$.

For the quadratic term in $\eqref{2-ln-1}$, from condition (L5),

$\displaystyle\frac{1}{2}\nu_{n}^{2}u^{\rm
T}\mathbf{B}_{n}u\geq\frac{C_{4}}{2}n\nu_{n}^{2}\|u\|_{2}^{2}.$ (20)

For the linear term $\mathbf{A}_{n}^{\rm T}u$ in $\eqref{2-ln-1}$, by the
independence of $(Z_{i},Y_{i})$ and $(Z_{j},Y_{j})$ for all $i\neq j$ and the
fact that $E(I(Y_{i}<U_{i}^{\rm T}\theta_{m}^{*})-\tau_{m}|Z_{i})=0$, we get

$\displaystyle E(A_{n}^{\rm T}A_{n})$
$\displaystyle=E\sum_{m=1}^{M}\,\sum_{i=1}^{n}\,\bigl{(}I(Y_{i}<U_{i}^{\rm
T}\theta_{m}^{*})-\tau_{m}\bigr{)}\,U_{i}^{\rm
T}\sum_{j=1}^{n}\,\bigl{(}I(Y_{j}<U_{j}^{\rm
T}\theta_{m}^{*})-\tau_{m}\bigr{)}\,U_{j}$
$\displaystyle=\sum_{m=1}^{M}\,\sum_{i=1}^{n}\,E\bigl{(}\bigl{(}I(Y_{i}<U_{i}^{\rm
T}\theta_{m}^{*})-\tau_{m}\bigr{)}^{2}\|U_{i}\|_{2}^{2}\bigl{)}\leq Cnp.$

Then, it follows that

$\displaystyle\|\mathbf{A}_{n}\|_{2}=O_{p}((np)^{\frac{1}{2}}),$ (21)

which implies that $\sup_{\|u\|_{2}\leq 1}|\mathbf{A}_{n}^{\rm
T}u|=O_{p}((np)^{\frac{1}{2}})$. Together with $\eqref{2-ln-1}$ and
$\eqref{2-quad-1}$, in probability,

$\displaystyle\inf_{\|u\|_{2}=1}\,l_{1n}(\theta^{*}+\nu_{n}u)$
$\displaystyle\geq\frac{C_{4}}{2}n\nu_{n}^{2}-C(np)^{\frac{1}{2}}\nu_{n}-Cn^{\frac{1}{2}}p^{\frac{3}{4}}{\nu_{n}}^{\frac{3}{2}}-Cnp^{\frac{1}{2}}\nu_{n}^{3}-Cn\lambda_{n}\,d_{nI}^{\frac{1}{2}}\nu_{n}^{\frac{1}{2}}s_{n}^{\frac{3}{4}}$
$\displaystyle\geq\frac{C_{4}}{2}n\nu_{n}\bigl{\\{}\nu_{n}-Cn^{-\frac{1}{2}}p^{\frac{1}{2}}-Cn^{-\frac{1}{2}}p^{\frac{3}{4}}{\nu_{n}}^{\frac{1}{2}}-Cp^{\frac{1}{2}}\nu_{n}^{2}-C\lambda_{n}d_{nI}^{\frac{1}{2}}s_{n}^{\frac{3}{4}}\nu_{n}^{-\frac{1}{2}}\bigr{\\}}.$
(22)

Now take $\nu_{n}=C_{0}(n^{-{\frac{1}{2}}}p^{\frac{1}{2}})$ where $C_{0}$ is a
sufficiently large constant. Under condition (L4), i.e.,
$\lambda_{n}d_{nI}^{\frac{1}{2}}=o(s_{n}^{-\frac{3}{4}}p_{n}^{\frac{3}{4}}n^{-\frac{3}{4}})$,
for the last three terms in $\eqref{2-order-2-2}$, we can check that
$n^{-\frac{1}{2}}p^{\frac{3}{4}}{\nu_{n}}^{\frac{1}{2}}\leq
Cn^{-\frac{1}{4}+\frac{1}{2}\alpha_{1}}\nu_{n}=o(\nu_{n})$,
$p^{\frac{1}{2}}\nu_{n}^{2}\leq
Cn^{-\frac{1}{2}+\alpha_{1}}\nu_{n}=o(\nu_{n})$, and
$\lambda_{n}\,d_{nI}^{\frac{1}{2}}s_{n}^{\frac{3}{4}}\nu_{n}^{-\frac{1}{2}}\leq
C\lambda_{n}\,d_{nI}^{\frac{1}{2}}s_{n}^{\frac{3}{4}}p_{n}^{-\frac{3}{4}}n^{\frac{3}{4}}\nu_{n}=o(\nu_{n})$.
Hence,

$\displaystyle\inf_{\|u\|_{2}=1}\,l_{1n}(\theta^{*}+\nu_{n}u)\geq
Cn\nu_{n}^{2}\to\infty\quad\text{in probability}.$ (23)

Therefore, in probability there exists a local minimizer $\hat{\theta}$ of
$L_{n}(\theta)$ such that $\|\hat{\theta}-\theta^{*}\|_{2}<\nu_{n}$. This
completes the proof of the theorem. $\Box$

Proof of Theorem 2. For the quantile function $Q_{n}(\theta)$, because
$\|\theta_{m}-\theta_{m}^{*}\|_{2}\,\leq\,\|\theta-\theta^{*}\|_{2}$ for all
$1\leq\,m\leq\,M$, by Lemma 1 with $\eta_{n}=\nu_{n}$, under conditions
(L1)-(L4), we have

$\displaystyle Q_{n}(\theta)-Q_{n}(\theta^{*})=\mathbf{A}_{n}^{\rm
T}(\theta-\theta^{*})+\frac{1}{2}(\theta-\theta^{*})^{\rm
T}\mathbf{B}_{n}(\theta-\theta^{*})+R_{n}(\theta),$ (24)

where
$\sup_{\|\theta-\theta^{*}\|_{2}\,\leq\,\nu_{n}}\,|R_{n}(\theta)|=O_{p}(n^{\frac{1}{2}}p^{\frac{3}{4}}{\nu_{n}}^{\frac{3}{2}})+O_{p}(np^{\frac{1}{2}}\nu_{n}^{3})$.

Let $\theta-\theta^{*}=\nu_{n}u$ where $\nu_{n}>0$ and
$u\in\mathbb{R}^{M(p+1)}$. Then $\|\theta-\theta^{*}\|_{2}\leq\nu_{n}$ if and
only if $\|u\|_{2}\leq 1$. Let $u=(u_{I}^{\rm T},u_{II}^{\rm T})^{\rm T}$,
where $u_{I}$ and $u_{II}$ are are subvectors of $u$ corresponding to the
index sets $I$ and $II$, respectively. Clearly,
$\|u\|_{2}^{2}=\|u_{I}\|_{2}^{2}+\|u_{II}\|_{2}^{2}$. Note that
$\theta_{I}=\theta_{I}^{*}+\nu_{n}u_{I}$ and
$\theta_{II}=\theta_{II}^{*}+\nu_{n}u_{II}=\nu_{n}u_{II}$.

Define the ball
$\tilde{\Theta}_{n}=\\{\theta=\theta^{*}+\nu_{n}u\in\Theta_{n}:\,\|u\|_{2}\leq
1\\}$ with $\nu_{n}=C_{0}n^{-\frac{1}{2}}p^{\frac{1}{2}}$. For any
$\theta=(\theta_{I}^{\rm T},\theta_{II}^{\rm T})^{\rm
T}\in\tilde{\Theta}_{n}$, we can see $\|(\theta_{I}^{\rm T},0^{\rm T})^{\rm
T}-\theta^{*}\|_{2}=\|\theta_{I}-\theta_{I}^{*}\|_{2}\leq\nu_{n}$, and
$\|\theta_{II}\|_{2}=\nu_{n}\|u_{II}\|_{2}$, where $\|u_{II}\|_{2}\leq 1$.

Consider that
$Q_{n}(\theta_{I},\theta_{II})-Q_{n}(\theta_{I},0)=Q_{n}(\theta_{I},\theta_{II})-Q_{n}(\theta^{*})-\bigl{(}Q_{n}(\theta_{I},0)-Q_{n}(\theta^{*})\bigr{)}$

$\displaystyle=\frac{1}{2}(0^{\rm T},\theta_{II}^{\rm T})\mathbf{B}_{n}(0^{\rm
T},\theta_{II}^{\rm T})^{\rm T}+\mathbf{A}_{n}^{\rm T}(0^{\rm
T},\theta_{II}^{\rm T})^{\rm T}-(\theta_{I}^{\rm T}-\theta_{I}^{*\rm T},0^{\rm
T})\mathbf{B}_{n}(0^{\rm T},\theta_{II}^{\rm T})^{\rm
T}+\bigl{(}R_{n}(\theta)-R_{n}(\theta_{I},0)\bigr{)}$ $\displaystyle\equiv
I_{n1}+I_{n2}+I_{n3}+r_{1n}(\theta)$

From $\eqref{2-an-2}$, we can see
$\sup_{\theta\in\tilde{\Theta}_{n}}\,|I_{n2}|\leq\|\mathbf{A}_{n}\|_{2}\|\theta_{II}\|_{2}=O_{p}((np)^{\frac{1}{2}})\nu_{n}\|u_{II}\|_{2}=O_{p}(p\|u_{II}\|_{2})$.
Under condition (L5), $\|\mathbf{B}_{n}(0^{\rm T},\theta_{II}^{\rm T})^{\rm
T}\|_{2}^{2}=(0^{\rm T},\theta_{II}^{\rm T})\mathbf{B}_{n}^{2}(0^{\rm
T},\theta_{II}^{\rm T})^{\rm T}\leq
n^{2}\lambda_{\max}^{2}(n^{-1}B_{n})\|\theta_{II}\|_{2}^{2}\leq
Cn^{2}\nu_{n}^{2}\|u_{II}\|_{2}^{2}$, we have
$\sup_{\theta\in\tilde{\Theta}_{n}}\,|I_{n3}|\leq\|\theta_{I}-\theta_{I}^{*}\|_{2}\|\mathbf{B}_{n}(0^{\rm
T},\theta_{II}^{\rm T})^{\rm T}\|_{2}\leq Cn\nu_{n}^{2}\leq Cp$. From
$\eqref{3-qn-0}$, we have
$\sup_{\theta\in\tilde{\Theta}_{n}}\,|r_{1n}(\theta)|=O_{p}(n^{-\frac{1}{4}}p^{\frac{3}{2}})$.
Hence,

$\displaystyle
Q_{n}(\theta_{I},\theta_{II})-Q_{n}(\theta_{I},0)=I_{n1}+r_{2n}(\theta),$ (25)

where
$\sup_{\theta\in\tilde{\Theta}_{n}}\,|r_{2n}(\theta)|=O_{p}(p\|u_{II}\|_{2})+O_{p}(n^{-\frac{1}{4}}p^{\frac{3}{2}})$.
Recall
$L_{n}(\theta_{I},\theta_{II})-L_{n}(\theta_{I},0)=Q_{n}(\theta_{I},\theta_{II})-Q_{n}(\theta_{I},0)+P_{n}(\gamma_{I},\gamma_{II})-P_{n}(\gamma_{I},0)$.
From $\eqref{3-qn}$, we get

$\displaystyle L_{n}(\theta_{I},\theta_{II})-L_{n}(\theta_{I},0)\geq
P_{n}(\gamma_{I},\gamma_{II})-P_{n}(\gamma_{I},0)+r_{2n}(\theta).$ (26)

Note that
$(n\lambda_{n})^{-1}\bigl{(}P(\gamma_{I},\gamma_{II})-P_{n}(\gamma_{I},0)\bigr{)}$

$\displaystyle=\sum_{j=1}^{p}\,\bigl{(}\sum_{m=1}^{M}\,\omega_{mj}|\gamma_{mj}|\bigr{)}^{\frac{1}{2}}-\sum_{j\in
J}\,\bigl{(}\sum_{m\in M_{j}}\,\omega_{mj}|\gamma_{mj}|\bigr{)}^{\frac{1}{2}}$
$\displaystyle=\sum_{j\in
J^{c}}\,\bigl{(}\sum_{m=1}^{M}\,\omega_{mj}|\gamma_{mj}|\bigr{)}^{\frac{1}{2}}+\sum_{j\in
J}\,\bigl{(}\sum_{m=1}^{M}\,\omega_{mj}|\gamma_{mj}|\bigr{)}^{\frac{1}{2}}-\sum_{j\in
J}\,\bigl{(}\sum_{m\in M_{j}}\,\omega_{mj}|\gamma_{mj}|\bigr{)}^{\frac{1}{2}}$
$\displaystyle\geq\sum_{j\in
J^{c}}\frac{\sum_{m=1}^{M}\,\omega_{mj}|\gamma_{mj}|}{2\bigl{(}\sum_{m=1}^{M}\,\omega_{mj}|\gamma_{mj}|\bigr{)}^{\frac{1}{2}}}+\sum_{j\in
J}\frac{\sum_{m\in
M_{j}^{c}}\omega_{mj}|\gamma_{mj}|}{2\bigl{(}\sum_{m=1}^{M}\omega_{mj}|\gamma_{mj}|\bigr{)}^{\frac{1}{2}}}.$
(27)

For all $\theta\in\tilde{\Theta}_{n}$, we have
$|\gamma_{mj}|\leq|\gamma_{mj}^{*}|+1\leq C$, which implies that
$\sum_{m=1}^{M}\omega_{mj}|\gamma_{mj}|\leq C\|\omega_{n}\|_{\infty}$ for all
$j=1,2,\cdots,p$, where $\|\omega_{n}\|_{\infty}=\max_{1\leq m\leq M,\,1\leq
j\leq p}|\omega_{mj}|$. Recall that $d_{nII}=\min_{(m,j)\in
II}\\{\omega_{mj}\\}\|\omega_{n}\|_{\infty}^{-\frac{1}{2}}$. From
$\eqref{3-pn-1}$, it follows that for all $\theta\in\tilde{\Theta}_{n}$,

$\displaystyle P_{n}(\gamma_{I},\gamma_{II})-P_{n}(\gamma_{I},0)$
$\displaystyle\geq
Cn\lambda_{n}\,\|\omega_{n}\|_{\infty}^{-\frac{1}{2}}\bigl{(}\sum_{j\in
J^{c}}\sum_{m=1}^{M}\,\omega_{mj}|\gamma_{mj}|+\sum_{j\in J}\sum_{m\in
M_{j}^{c}}\omega_{mj}|\gamma_{mj}|\bigr{)}$
$\displaystyle=Cn\lambda_{n}\,\|\omega_{n}\|_{\infty}^{-\frac{1}{2}}\sum_{(m,j)\in
II}\,\omega_{mj}|\gamma_{mj}|\geq Cn\lambda_{n}\,d_{nII}\sum_{(m,j)\in
II}|\gamma_{mj}|$ $\displaystyle\geq
Cn\lambda_{n}\,d_{nII}\|\gamma_{II}\|_{2}=C\lambda_{n}\,d_{nII}(np)^{\frac{1}{2}}\|u_{II}\|_{2}.$
(28)

Define
$\Omega_{n}=\\{\theta=\theta^{*}+\nu_{n}u\in\tilde{\Theta}_{n}:\|u_{II}\|_{2}>0\\}$
and
$\Omega_{n}^{c}=\\{\theta=\theta^{*}+\nu_{n}u\in\tilde{\Theta}_{n}:u_{II}=0\\}$.
Clearly, $\tilde{\Theta}_{n}=\Omega_{n}\cup\Omega_{n}^{c}$. From
$\eqref{3-ln-0}$ and $\eqref{3-pn-0}$, we obtain in probability

$\displaystyle\inf_{\theta\in\Omega_{n}}\,\bigl{(}L_{n}(\theta_{I},\theta_{II})-L_{n}(\theta_{I},0)\bigr{)}$
$\displaystyle\geq\inf_{\theta\in\Omega_{n}}\,\bigl{(}P_{n}(\gamma_{I},\gamma_{II})-P_{n}(\gamma_{I},0)\bigr{)}-\sup_{\theta\in\tilde{\Theta}_{n}}\,|r_{2n}(\theta)|$
$\displaystyle\geq\tilde{C}_{1}\lambda_{n}\,d_{nII}(np)^{\frac{1}{2}}\|u_{II}\|_{2}-\tilde{C}_{2}p\|u_{II}\|_{2}-\tilde{C}_{2}n^{-\frac{1}{4}}p^{\frac{3}{2}}$
$\displaystyle\geq
p\biggl{(}\|u_{II}\|_{2}\bigl{(}\tilde{C}_{1}\lambda_{n}\,d_{nII}n^{\frac{1}{2}}p^{-\frac{1}{2}}-\tilde{C}_{2}\bigr{)}-\tilde{C}_{2}n^{-\frac{1}{4}}p^{\frac{1}{2}}\biggr{)}.$

where $\tilde{C}_{1}$ and $\tilde{C}_{2}$ are positive constants. Under the
given conditions,
$\lambda_{n}\,d_{nII}n^{\frac{1}{2}}p^{-\frac{1}{2}}\to\infty$ and
$n^{-\frac{1}{4}}p^{\frac{1}{2}}\to 0$ as $n\to\infty$. Hence,
$\inf_{\theta\in\Omega_{n}}\,\bigl{(}L_{n}(\theta_{I},\theta_{II})-L_{n}(\theta_{I},0)\bigr{)}>0$
in probability. Thus,
$\inf_{\theta\in\Omega_{n}}\,L_{n}(\theta)\geq\inf_{\theta\in\Omega_{n}}\,\bigl{(}L_{n}(\theta_{I},\theta_{II})-L_{n}(\theta_{I},0)\bigr{)}+\inf_{\theta\in\Omega_{n}}L_{n}(\theta_{I},0)>\inf_{\theta\in\Omega_{n}}L_{n}(\theta_{I},0)=\inf_{\theta\in\tilde{\Theta}_{n}}L_{n}(\theta_{I},0)\geq\inf_{\theta\in\tilde{\Theta}_{n}}L_{n}(\theta)$,
which implies that
$\inf_{\theta\in\tilde{\Theta}_{n}}\,L_{n}(\theta)=\inf_{\theta\in\Omega_{n}^{c}}\,L_{n}(\theta)$.
Therefore, the minimal value point of $L_{n}(\theta)$ on $\tilde{\Theta}_{n}$
only lies in its subset $\Omega_{n}^{c}$.

From Theorem 1, we know that in probability
$\hat{\theta}\in\tilde{\Theta}_{n}$ and that $\hat{\theta}$ is a local
minimizer of $L_{n}(\theta)$. Hence, $\hat{\theta}\in\Omega_{n}^{c}$ in
probability, which implies that $\hat{\gamma}_{II}=0$. $\Box$

Proof of Theorem 3. Let $\theta_{I}-\theta_{I}^{*}=\nu_{n}u$ where $\nu_{n}>0$
and $u\in\mathbb{R}^{M+s_{n}}$. Because
$\|\theta_{mI}-\theta_{mI}^{*}\|_{2}\,\leq\,\|\theta_{I}-\theta_{I}^{*}\|_{2}=\nu_{n}\|u\|_{2}$
for all $1\leq\,m\leq\,M$, and due to conditions (L1)-(L3) and that
$0<\alpha_{0}<\alpha_{1}<\frac{1}{6}$, Lemma 1 implies that

$\displaystyle
Q_{n}(\theta_{I},0)-Q_{n}(\theta_{I}^{*},0)=Q_{n}(\theta_{I}^{*}+\nu_{n}u,0)-Q_{n}(\theta_{I}^{*},0)=\nu_{n}\mathbf{A}_{nI}^{\rm
T}u+\frac{1}{2}\nu_{n}^{2}u^{\rm T}\mathbf{B}_{nI}u+r_{n}(u),$ (29)

where
$\sup_{\|u\|_{2}\,\leq\,1}\,|r_{n}(u)|=O_{p}(n^{\frac{1}{2}}s_{n}^{\frac{3}{4}}\nu_{n}^{\frac{3}{2}})+O_{p}(ns_{n}^{\frac{1}{2}}\nu_{n}^{3})$,
$\mathbf{B}_{nI}$ is given in $\eqref{4-bn}$,

$\displaystyle\mathbf{A}_{nI}$ $\displaystyle=(A_{1}^{\rm T},\cdots,A_{M}^{\rm
T})^{\rm T}\quad\text{with}\quad
A_{m}=\sum_{i=1}^{n}\,\bigl{(}I(Y_{i}<U_{i}^{\rm
T}\theta_{m}^{*})-\tau_{m}\bigr{)}\,U_{imI},$ (30)

and $U_{imI}$ is given in Section 3. Note that $A_{nI}$ and $B_{nI}$ are the
sub-vector(matrix) of $\mathbf{A}_{n}$ and $\mathbf{B}_{n}$ in
$\eqref{0-anbn}$ corresponding to the index set $I$, respectively.

For $P_{n}(\gamma)$, define $p_{n}(\gamma)=P_{n}(\gamma)-P_{n}(\gamma^{*})$.
Let $u_{I}$ be the subvector of $u$ corresponding to the subvector
$\gamma_{I}$ in $\theta_{I}$. Then, we have

$\displaystyle p_{n}(\gamma_{I},0)=p_{n}(\gamma_{I}^{*}+\nu_{n}u_{I},0)$
$\displaystyle=n\lambda_{n}\,\sum_{j\in J}\,\bigl{(}\sum_{m\in
M_{j}}\,\omega_{mj}|\gamma_{mj}|\bigr{)}^{\frac{1}{2}}-n\lambda_{n}\,\sum_{j\in
J}\,\bigl{(}\sum_{m\in
M_{j}}\,\omega_{mj}|\gamma_{mj}^{*}|\bigr{)}^{\frac{1}{2}},$

which is $p_{1n}(\gamma)$ given in $\eqref{2-p1n}$. By $\eqref{2-tpn}$ in the
proof of Theorem 1, we have
$|p_{n}(\gamma_{I},0)|=|p_{n}(\gamma_{I}^{*}+\nu_{n}u_{I},0)|\leq
n\lambda_{n}d_{nI}^{\frac{1}{2}}\nu_{n}^{\frac{1}{2}}s_{n}^{\frac{3}{4}}\|u\|_{2}^{\frac{1}{2}}$.
This, combined with $\eqref{4-qn}$, implies that

$\displaystyle
L_{n}(\theta_{I}^{*}+\nu_{n}u,0)-L_{n}(\theta_{I}^{*},0)=\nu_{n}\mathbf{A}_{nI}^{\rm
T}u+\frac{1}{2}\nu_{n}^{2}u^{\rm T}\mathbf{B}_{nI}u+r_{l}(u),$ (31)

where
$\sup_{\|u\|_{2}\,\leq\,1}\,|r_{l}(u)|=O_{p}(n^{\frac{1}{2}}s_{n}^{\frac{3}{4}}{\nu_{n}}^{\frac{3}{2}})+O_{p}(ns_{n}^{\frac{1}{2}}\nu_{n}^{3})+O(n\lambda_{n}\,d_{nI}^{\frac{1}{2}}\nu_{n}^{\frac{1}{2}}s_{n}^{\frac{3}{4}})$.

Define the ball $\Theta_{nI}=\\{(\theta_{I}^{\rm T},0^{\rm T})^{\rm
T}\,:\,\theta_{I}-\theta_{I}^{*}=\nu_{n}\,u\\}$ with
$\nu_{n}=C_{0}(n^{-1}p_{n})^{\frac{1}{2}}$, where $u\in\mathbb{R}^{M+s_{n}}$
and $C_{0}$ is a positive constant. Given that
$\nu_{n}=C_{0}(n^{-1}p_{n})^{\frac{1}{2}}$,
$\lambda_{n}d_{nI}^{\frac{1}{2}}=O(n^{-1}p_{n}^{\frac{1}{2}})$, and
$0<\alpha_{0}<\alpha_{1}<\frac{1}{6}$ we see that
$n^{\frac{1}{2}}s_{n}^{\frac{3}{4}}{\nu_{n}}^{\frac{3}{2}}\leq
Cn^{-\frac{1}{4}}s_{n}^{\frac{3}{4}}p_{n}^{\frac{3}{4}}=o_{p}(1)$,
$n\sqrt{s_{n}}\nu_{n}^{3}\leq
Cn^{-\frac{1}{2}}s_{n}^{\frac{1}{2}}p_{n}^{\frac{3}{2}}=o(n^{-\frac{1}{4}}s_{n}^{\frac{3}{4}}p_{n}^{\frac{3}{4}})=o_{p}(1)$
and
$n\lambda_{n}\,d_{nI}^{\frac{1}{2}}\nu_{n}^{\frac{1}{2}}s_{n}^{\frac{3}{4}}\leq\lambda_{n}\,d_{nI}^{\frac{1}{2}}\,n^{\frac{3}{4}}s_{n}^{\frac{3}{4}}p_{n}^{\frac{1}{4}}=O(n^{-\frac{1}{4}}s_{n}^{\frac{3}{4}}p_{n}^{\frac{3}{4}})=o_{p}(1)$.
Hence, considering $\eqref{4-ln-1}$, for all $(\theta_{I}^{\rm T},0^{\rm
T})^{\rm T}\in\Theta_{nI}$ we have
$L_{n}(\theta_{I},0)-L_{n}(\theta_{I}^{*},0)=\mathbf{A}_{nI}^{\rm
T}(\theta_{I}-\theta_{I}^{*})+\frac{1}{2}\,(\theta_{I}-\theta_{I}^{*})^{\rm
T}\,\mathbf{B}_{nI}\,(\theta_{I}-\theta_{I}^{*})+o_{p}(1)$, which implies that
for all $(\theta_{I}^{\rm T},0^{\rm T})^{\rm T}\in\Theta_{nI}$,

$\displaystyle
L_{n}(\theta_{I},0)-L_{n}(\theta_{I}^{*},0)+\frac{1}{2}\,\mathbf{A}_{nI}^{\rm
T}\mathbf{B}_{nI}^{-1}\mathbf{A}_{nI}$ (32) $\displaystyle=$
$\displaystyle\frac{1}{2}\bigl{(}\mathbf{B}_{nI}^{-\frac{1}{2}}\,\mathbf{A}_{nI}+\mathbf{B}_{nI}^{\frac{1}{2}}(\theta_{I}-\theta_{I}^{*})\bigr{)}^{\rm
T}\bigl{(}\mathbf{B}_{nI}^{-\frac{1}{2}}\,\mathbf{A}_{nI}+\mathbf{B}_{nI}^{\frac{1}{2}}(\theta_{I}-\theta_{I}^{*})\bigr{)}+o_{p}(1).$

By Theorem 1 and Theorem 2, we know that in probability a local minimizer
$\hat{\theta}$ of $L_{n}(\theta)-L_{n}(\theta^{*})$ lies in the ball
$\Theta_{nI}$, which implies that $\hat{\theta}=(\hat{\theta}_{I}^{\rm
T},0^{\rm T})^{\rm T}\in\Theta_{nI}$ in probability. Hence, from $\eqref{4-ln-
eq}$, we have

$\displaystyle\mathbf{B}_{nI}^{\frac{1}{2}}(\hat{\theta}_{I}-\theta_{I}^{*})=-\mathbf{B}_{nI}^{-\frac{1}{2}}\,\mathbf{A}_{nI}+o_{p}(1)t,$
(33)

where $t\in\mathbb{R}^{M+s_{n}}$ is an unit vector. Since $t^{\rm
T}\mathbf{B}_{nI}t\leq\lambda_{\max}(\mathbf{B}_{nI})\leq\lambda_{\max}(\mathbf{B}_{n})\leq
Cn$ by condition (L5), we have
$\mathbf{B}_{nI}^{\frac{1}{2}}t=O_{p}(n^{\frac{1}{2}})t$. Multiplying both
sides of $\eqref{4-11}$ by a vector $b^{\rm T}\mathbf{B}_{nI}^{\frac{1}{2}}$,
where $b\in\mathbb{R}^{M+s_{n}}$ is any unit vector, we obtain

$\displaystyle b^{\rm
T}\mathbf{B}_{nI}(\hat{\theta}_{I}-\theta_{I}^{*})=-b^{\rm
T}\mathbf{A}_{nI}+o_{p}(n^{1/2}).$ (34)

Let $\xi_{n}=b^{\rm T}\mathbf{A}_{nI}$, and write $b=(b_{1}^{\rm
T},\cdots,\,b_{M}^{\rm T})^{\rm T}$ where $b_{m}$ is the subvector of $b$
corresponding to the subvector $\theta_{mI}^{*}$ of $\theta_{I}^{*}$. By the
definition of $\mathbf{A}_{nI}$ in $\eqref{4-an}$, we see that

$\displaystyle\xi_{n}$
$\displaystyle=\sum_{m=1}^{M}\,\sum_{i=1}^{n}\,\bigl{(}\psi_{mi}^{*}\,U_{imI}^{\rm
T}b_{m}\bigr{)}=\sum_{i=1}^{n}\,\zeta_{i},$

where $\zeta_{i}=\sum_{m=1}^{M}\,\psi_{mi}^{*}\,U_{imI}^{\rm T}b_{m}$ with
$\psi_{mi}^{*}=I(Y_{i}<U_{imI}^{\rm T}\theta_{mI}^{*})-\tau_{m}$. Clearly,
$\\{\zeta_{i},\,i=1,\cdots,n\\}$ is an independent sequence. Next, we will
verify that $\xi_{n}$ satisfies the Lindeberg’s condition

$\displaystyle\sigma_{n}^{-2}\,\sum_{i=1}^{n}\,E\bigl{(}\zeta_{i}^{2}I(|\zeta_{i}|\geq\sigma_{n})\bigr{)}\to
0,$ (35)

where $\sigma_{n}^{2}=Var(\xi_{n})$.

For $\zeta_{i}=\sum_{m=1}^{M}\,\psi_{mi}^{*}\,U_{imI}^{\rm T}b_{m}$, it is
easy to see that
$E(\zeta_{i})=\sum_{m=1}^{M}\,E\bigl{(}E(\psi_{mi}^{*}|Z_{i})\,U_{imI}^{\rm
T}b_{m}\bigr{)}=0$ from the fact $E(\psi_{mi}^{*}|Z_{i})=F(U_{imI}^{\rm
T}\theta_{mI}^{*}|Z_{i})-\tau_{m}=F(U_{i}^{\rm
T}\theta_{m}^{*}|Z_{i})-\tau_{m}=0$, and that

$\displaystyle E(\zeta_{i}^{2})$
$\displaystyle=E\bigl{(}\sum_{m=1}^{M}\,\psi_{mi}^{*}\,U_{imI}^{\rm
T}b_{m}\bigr{)}^{2}=\sum_{m=1}^{M}\,\sum_{l=1}^{M}\,E\bigl{(}E(\psi_{mi}^{*}\psi_{li}^{*}|Z_{i})\,b_{m}^{\rm
T}(U_{imI}U_{ilI}^{\rm T})b_{l}\bigr{)}$
$\displaystyle=\sum_{m=1}^{M}\,\sum_{l=1}^{M}\,E\bigl{(}(\min(\tau_{m},\tau_{l})-\tau_{m}\tau_{l})\,b_{m}^{\rm
T}(U_{imI}U_{ilI}^{\rm T})b_{l}\bigr{)}=b^{\rm T}\mathbf{\Sigma}_{n}b,$

where $\mathbf{\Sigma}_{n}$ is given in $\eqref{4-sigma}$. Hence, we obtain
that $E(\xi_{n})=\sum_{i=1}^{n}E(\zeta_{i})=0$ and that, by independence
$Var(\xi_{n})=\sum_{i=1}^{n}E(\zeta_{i}^{2})=nb^{\rm T}\mathbf{\Sigma}_{n}b$.
Under condition (L6), we obtain

$\displaystyle\sigma_{n}^{2}=nb^{\rm T}\mathbf{\Sigma}_{n}b\geq
n\lambda_{\min}(\mathbf{\Sigma}_{n})b^{\rm T}b\geq Cn.$ (36)

Note that $(U_{imI}^{\rm
T}b_{m})^{2}\leq\|U_{imI}\|_{2}^{2}\|b_{m}\|_{2}^{2}\leq
Cs_{n}\|b\|_{2}^{2}\leq Cs_{n}$ under condition (L2). By Cauchy-Schwartz’s
inequality, $\zeta_{i}^{2}=(\sum_{m=1}^{M}\,\psi_{mi}^{*}\,U_{imI}^{\rm
T}b_{m})^{2}\leq M\sum_{m=1}^{M}\,(\psi_{mi}^{*}\,U_{imI}^{\rm
T}b_{m})^{2}\leq Cs_{n}\sum_{m=1}^{M}\,(\psi_{mi}^{*})^{2}$, which implies
that $\zeta_{i}^{2}\leq Cs_{n}$ from the fact $|\psi_{mi}^{*}|\leq
1+\tau_{m}\leq 2$. Hence, we have

$\displaystyle\sum_{i=1}^{n}\,E\bigl{(}\zeta_{i}^{2}I(|\zeta_{i}|\geq\sigma_{n})\bigr{)}$
$\displaystyle\leq
Cs_{n}\sum_{i=1}^{n}\,E\bigl{(}I(\zeta_{i}^{2}\geq\sigma_{n}^{2})\bigr{)}\leq
Cs_{n}\sum_{i=1}^{n}\,P\bigl{(}Cs_{n}\sum_{m=1}^{M}\,(\psi_{mi}^{*})^{2}\geq\sigma_{n}^{2}\bigr{)}$
$\displaystyle\leq
Cs_{n}\sum_{i=1}^{n}\,E\bigl{(}Cs_{n}\sum_{m=1}^{M}\,(\psi_{mi}^{*})^{2}\bigr{)}\,\sigma_{n}^{-2}\leq
Cns_{n}^{2}\sigma_{n}^{-2}.$

Together with $\eqref{4-sigma2}$, this implies that, as $n\to\infty$,

$\displaystyle\sigma_{n}^{-2}\,\sum_{i=1}^{n}\,E\bigl{(}\zeta_{i}^{2}I(|\zeta_{i}|\geq\sigma_{n})\bigr{)}\leq
Cns_{n}^{2}\sigma_{n}^{-4}\leq C\frac{s_{n}^{2}}{n}\leq Cn^{-1+2\alpha_{0}}\to
0,$

which shows that Lindeberg’s condition $\eqref{4-ldbg}$ holds. Hence,

$\displaystyle(nb^{\rm T}\mathbf{\Sigma}_{n}b)^{-\frac{1}{2}}\,b^{\rm
T}\mathbf{A}_{nI}=\frac{b^{\rm
T}\mathbf{A}_{nI}}{\sigma_{n}}=\frac{\xi_{n}-E(\xi_{n})}{\sigma_{n}}\to
N(0,1).$

Together with $\eqref{4-13}$ and $\eqref{4-sigma2}$, this implies that

$\displaystyle(nb^{\rm T}\mathbf{\Sigma}_{n}b)^{-\frac{1}{2}}\,b^{\rm
T}\mathbf{B}_{nI}\bigl{(}\hat{\theta}_{I}-\theta_{I}^{*}\bigr{)}=-(nb^{\rm
T}\mathbf{\Sigma}_{n}b)^{-\frac{1}{2}}\,b^{\rm T}\mathbf{A}_{nI}+o_{p}(1)\to
N(0,1).$

This completes the proof of the theorem. $\Box$

## References

## References

  * Avery et al. [2011] Avery, C.L., He, Q., North, K.E., Ambite, J.L., Boerwinkle, E., Fornage, M., Hindorff, L.A., Kooperberg, C., Meigs, J.B., Pankow, J.S., et al., 2011\.  A phenomics-based strategy identifies loci on APOC1, BRAP, and PLCG1 associated with metabolic syndrome phenotype domains.  PLoS Genet 7, e1002322.
  * Fan and Li [2001] Fan, J., Li, R., 2001.  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American Statistical Association 96, 1348–1360.
  * Gasso et al. [2009] Gasso, G., Rakotomamonjy, A., Canu, S., 2009.  Recovering sparse signals with a certain family of nonconvex penalties and DC programming.  IEEE Transactions on Signal Processing 57, 4686–4698.
  * Holland [2000] Holland, E.C., 2000.  Glioblastoma multiforme: the terminator.  Proceedings of the National Academy of Sciences 97, 6242–6244.
  * Huang et al. [2009] Huang, J., Ma, S., Xie, H., Zhang, C.H., 2009.  A group bridge approach for variable selection.  Biometrika 96, 339–355.
  * Jiang et al. [2014] Jiang, L., Bondell, H.D., Wang, H.J., 2014.  Interquantile shrinkage and variable selection in quantile regression.  Computational Statistics & Data Analysis 69, 208–219.
  * Jiang et al. [2013] Jiang, L., Wang, H.J., Bondell, H.D., 2013.  Interquantile shrinkage in regression models.  Journal of Computational and Graphical Statistics 22, 970–986.
  * Koenker [2004] Koenker, R., 2004.  Quantile regression for longitudinal data.  Journal of Multivariate Analysis 91, 74–89.
  * Koenker [2015] Koenker, R., 2015.  quantreg: Quantile Regression.  URL: http://CRAN.R-project.org/package=quantreg. R package version 5.11.
  * Koenker and Bassett [1978] Koenker, R., Bassett, J.G., 1978\.  Regression quantiles.  Econometrica 46, 33–50.
  * Landmark-Høyvik et al. [2013] Landmark-Høyvik, H., Dumeaux, V., Nebdal, D., Lund, E., Tost, J., Kamatani, Y., Renault, V., Børresen-Dale, A.L., Kristensen, V., Edvardsen, H., 2013\.  Genome-wide association study in breast cancer survivors reveals SNPs associated with gene expression of genes belonging to MHC class I and II.  Genomics 102, 278–287.
  * Li and Zhu [2008] Li, Y., Zhu, J., 2008.  L1-norm quantile regression.  Journal of Computational and Graphical Statistics 17, 163–185.
  * Mazumder et al. [2011] Mazumder, R., Friedman, J.H., Hastie, T., 2011.  Sparsenet: Coordinate descent with nonconvex penalties.  Journal of the American Statistical Association 106, 1125–1138.
  * Peng et al. [2014] Peng, L., Xu, J., Kutner, N., 2014.  Shrinkage estimation of varying covariate effects based on quantile regression.  Statistics and Computing 24, 853–869.
  * Puputti et al. [2006] Puputti, M., Tynninen, O., Sihto, H., Blom, T., Mäenpää, H., Isola, J., Paetau, A., Joensuu, H., Nupponen, N.N., 2006.  Amplification of KIT, PDGFRA, VEGFR2, and EGFR in gliomas.  Molecular Cancer Research 4, 927–934.
  * Tibshirani [1996] Tibshirani, R., 1996.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B 58, 267–288.
  * Tibshirani et al. [2005] Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., Knight, K., 2005.  Sparsity and smoothness via the fused lasso.  Journal of the Royal Statistical Society: Series B 67, 91–108.
  * Wang et al. [2012] Wang, L., Wu, Y., Li, R., 2012\.  Quantile regression for analyzing heterogeneity in ultra-high dimension.  Journal of the American Statistical Association 107, 214–222.
  * Wang et al. [2009] Wang, S., Nan, B., Zhu, N., Zhu, J., 2009.  Hierarchically penalized Cox regression with grouped variables.  Biometrika 96, 307–322.
  * Wu and Liu [2009] Wu, Y., Liu, Y., 2009.  Variable selection in quantile regression.  Statistica Sinica 19, 801–817.
  * Zhao and Xiao [2014] Zhao, Z., Xiao, Z., 2014.  Efficient regressions via optimally combining quantile information.  Econometric Theory 30, 1272–1314.
  * Zou [2006] Zou, H., 2006.  The adaptive lasso and its oracle properties.  Journal of the American Statistical Association 101, 1418–1429.
  * Zou and Yuan [2008] Zou, H., Yuan, M., 2008.  Regularized simultaneous model selection in multiple quantiles regression.  Computational Statistics & Data Analysis 52, 5296–5304.