Penalized regression with multiple loss functions and selection by vote

This article considers a linear model in a high dimensional data scenario. We propose a process which uses multiple loss functions both to select relevant predictors and to estimate parameters, and study its asymptotic properties. Variable selection is conducted by a procedure called"vote", which aggregates results from penalized loss functions. Using multiple objective functions separately simplifies algorithms and allows parallel computing, which is convenient and fast. As a special example we consider a quantile regression model, which optimally combines multiple quantile levels. We show that the resulting estimators for the parameter vector are asymptotically efficient. Simulations and a data application confirm the three main advantages of our approach: (a) reducing the false discovery rate of variable selection; (b) improving the quality of parameter estimation; (c) increasing the efficiency of computation.


Introduction
Consider a linear model where Y = (Y 1 , . . . , Y n ) T is an n-dimensional vector of responses, X = (X 1 , . . . , X n ) T is an n×p design matrix of predictors with X i = (X i1 , . . . , X ip ) T for i = 1, . . . , n, ϑ = (ϑ 1 , . . . , ϑ p ) T is a p-dimensional vector of parameters, and ε = (ε 1 , . . . , ε n ) T is an n-dimensional vector of independent and identically distributed random errors, which is independent of X. We study a scenario with high dimensional data where the dimension of the predictors p = p n is allowed to be greater than the sample size n and tend to infinity as n increases. We further assume that the model is sparse, i.e. only a fraction of the predictors significantly affects the response, while the parameters of the other predictors are zero.
One possible approach to selecting the important predictors and to estimate parameters is penalized regression. Various types of penalties, such as the Lasso (Tibshirani, 1996), the smoothly clipped absolute deviation (SCAD) penalty (Fan & Li, 2001), and the adaptive Lasso (Zou, 2006), have been applied to least squares regression. Penalized quantile regression has been considered because of its robustness (Fan et al., 2014;Wang et al., 2012;Wu & Liu, 2009). In addition to these articles working with a single loss function, Zou & Yuan (2008) introduced a composite quantile regression approach, which combines multiple quantile loss functions. They used an adaptive Lasso penalty to detect sparsity in linear models with a fixed number of parameters. Each of the above loss functions only works well for certain classes of error distributions. To obtain universal optimality, Bradic et al. (2011) developed a composite quasi-likelihood function, which approximates the log-likelihood function of the random error by a weighted linear combination of convex loss functions, and adopts a weighted Lasso penalty. The use of multiple loss functions in Zou & Yuan (2008) and Bradic et al. (2011) improves the efficiency of estimation. However, it does not necessarily reduce the false discovery rate in model selection. In fact, a loss function that yields an efficient estimator for a parameter may not be optimal for variable selection. For example, in Section 5.1 of Bradic et al. (2011), the simulation results show that penalized least squares regression can have higher false discovery rates than methods using other loss functions, even though it attains the Cramér-Rao bound because the error distribution is normal. Moreover, solving a penalized combination of multiple loss functions can be computationally complex, requiring complicated algorithms. For high dimensional data, the computation becomes increasingly intensive as the number of loss functions rises (Bradic et al., 2011, Section 4.2).
To make full use of multiple loss functions and facilitate computation, this article proposes a two-step process for the linear model (1): in the first step, the "vote procedure", the relevant predictors are determined. Then, in the second step, the parameters of the selected variables are estimated. This is different from the methods above, which generate a sparse estimator in one step by minimizing a penalized objective function. Our method uses multiple loss functions to which a weighted Lasso penalty is added. This yields different sparse (preliminary) estimators for the parameter vector ϑ. If a component of ϑ is identified as nonzero by a sufficient number of models, i.e. it has received enough "votes", the corresponding predictor is included in the final model. Our approach has smaller variance than most existing methods and excludes unimportant predictors more effectively. It requires no sophisticated algorithms and allows parallel computing to reduce processing time. A further advantage is that one can use different loss functions for variable selection and for parameter estimation, which brings more flexibility.
In the next section we introduce our approach, including the vote procedure, and study its asymptotic properties. Section 3 focuses on a special example, which optimally combines multiple quantile levels, and proves it is asymptotically efficient when the number of loss functions tends to infinity. Our method is compared with other competing methods by means of simulations in Section 4. In Section 5 we illustrate our approach with a real data analysis. Section 6 concludes the article with a brief summary and a discussion of further questions. All proofs are in the Appendix.

Variable selection and parameter estimation with multiple loss functions
In the following we will, for convenience of notation, let the lower case letter c represent a generic constant. The lower case letter r denotes a vector of a proper length whose components all equal one and the lower case letter v represents a unit vector of a proper length.
For a matrix B, let λ max (B) and λ min (B) denote its maximum and minimum eigenvalues.
, B ∞ equal to the largest absolute value of the entries in B and Assume that in the model (1), the parameter vector ϑ is sparse, i.e. there exists a set Q n ⊂ {1, 2, . . . , p n } such that ϑ j = 0 for any j ∈ Q n and ϑ j = 0 for any j ∈ Q c n = {1, 2, . . . , p n }\Q n . Without loss of generality, let Q n = {1, 2, . . . , q n } for some positive sequence q n < p n . To identify the set Q n we consider K different loss functions ℓ 1 (·), . . . , ℓ K (·). A weighted Lasso penalty with weights d k1 , . . . , d kpn and a tuning parameter λ n,k , is applied to the k th loss function. Then we obtain preliminary estimators for k = 1, . . . , K, where θ = (θ 1 , . . . , θ pn ) T . For some positive integer α ≤ K, the set Q n is estimated by the vote procedure, i.e., where I(·) is the indicator function and α is a threshold. This means that the j th component ϑ j of ϑ is included in the model if it receives α or more votes, i.e. at least α of the K estimates ϑ 1j , . . . , ϑ Kj are nonzero. In Theorem 1, we will see that the choice of α does not affect the result of variable selection asymptotically. When the sample size is finite we recommend cross validation to determine α for simulations and data analysis. A similar vote procedure was employed by Meinshausen & Bühlmann (2010) to develop stability selection and by Chen & Xie (2014) to apply the split and conquer strategy to penalized regression.
Unlike these two articles, which calculate multiple estimators on subsets of data, our method obtains each preliminary estimator (2) on the whole data set. This is desirable when the sample size is is not too large. In addition, the same penalty tuning parameter is used for all estimators in Meinshausen & Bühlmann (2010), while we choose a tuning parameter λ n,k for each estimator independently.
In the above process multiple loss functions are used not only to increase the efficiency of parameter estimation by the weighted combination (5) but also to improve the result of variable selection by the vote procedure (3). For the calculation of the estimators in (2) and (4), we only need algorithms for the K single objective functions. The multiple minimization procedures can then be conducted in parallel; see Section 4 for details. The two estimation steps may use different sets of loss functions. For example, if the error distribution is thought to be normal, we can use multiple loss functions for (2) and only the quadratic loss function for (4). For simplicity of notation, in the rest of this section and the next section we assume . Except for minor differences in notation, the conclusions are exactly the same if we drop these assumptions.
Assumptions 1 and 2 regulate the loss functions in (2) and (4). Common loss functions, such as the square function and the check function, satisfy these conditions. Assumption 3 is a standard condition on the growth rate of the model size for linear models with a diverging number of parameters, which can also be found in Bradic et al. (2011) and Wang et al. (2012), among others. Assumptions 4 and 5 guarantee good behaviour of the design matrix. Assumptions 6 and 7 are imposed on the weighted Lasso penalty to ensure important predictors can be detected. Assumption 8 excludes situations where the values of nonzero parameters decay too fast. Conditions similar to Assumptions 1-7 were required in Bradic et al. (2011) for penalized regression with a weighted linear combination of loss functions, while Assumption 8 is necessary for the nonzero parameters to be identified with probability approaching one by the preliminary estimators (2) and vote procedure (3).
In Theorem 1, we first state that, with probability close to one, the preliminary estimator which is the minimizer of the penalized objective function in the set Θ n = {θ ∈ R pn : θ j = 0 for j / ∈ Q n }, for k = 1, . . . , K. This indicates that ϑ k can exclude the unimportant variables. Then we show that by aggregating multiple such preliminary estimators, i.e. ϑ 1 , . . . , ϑ K , the vote procedure (3) owns model selection consistency. This means, with probability tending to one, the procedure can recover the index set of the nonzero parameters Q n .
We now specify the optimal weights, which minimizes the asymptotic variance in Theorem 2, for the estimator in (5) and show that the limiting distribution is not changed if we replace the optimal weights by a consistent estimator.
Corollary 1 Suppose the conditions in Theorem 2 hold. Then the optimal weight vector for the weighted estimator ϑ Qn (w) in Theorem 2 is w * = (r T H −1 r) −1 (H −1 r). In addition, if w * converges to w * in probability as n → ∞, the estimator ϑ Qn ( w * ) converges in distribution as Theorems 1 and 2 establish the model selection consistency and the asymptotic normality of the estimator with arbitrary weights when multiple loss functions are used for both variable selection and parameter estimation. Corollary 1 gives the optimal weights and the asymptotic distribution of the estimator (5)

An example
A special example of the weighted estimator (5) is the estimator from a quantile regression model optimally combining multiple quantile levels, which was considered in Section 3 of Zhao & Xiao (2014) for a fixed number of parameters without sparsity. These authors used and β k is the τ k quantile of the random error ε 1 , which can be estimated along with the slope vector ϑ as an additional parameter. In this scenario, the (i, j) th component of the matrix H in Corollary 1 is where f (·) is the density function of the random error ε 1 . We show the asymptotic efficiency of ϑ under the following assumption, which regulates the error distribution.

Assumption 9
The error density f (·) is positive, twice differentiable and bounded over the is the distribution function of ε 1 . For some constant The above assumption is satisfied by most of common continuous distributions. The value of ν 3 is related to the support of the error density. For example, the constant ν 3 equals zero if the support is (−∞, ∞). When the support has the form [s 1 , s 2 ], (−∞, s 2 ] or [s 1 , ∞) for some constants s 1 < s 2 , it is easy to see that ν 3 is infinite.
The following theorem gives the limit of the asymptotic variance (r T H −1 r) −1 of the optimally weighted estimator ϑ Qn (w * ) from Corollary 1 when the number of quantiles K tends to infinity.
Theorem 3 Consider the asymptotic variance in Corollary 1 and suppose Assumption 9 holds true, then the reciprocal of that variance satisfies where ν 3 is the constant in Assumption 9.
In the above conclusion, the integral on the right-hand side is the Fisher information. Consequently, we see that (r T H −1 r) coincides with the Fisher information, i.e. we have asymptotic efficiency, when the number of quantiles tends to infinity and ν 3 equals zero. This is the case, for example, if the error distribution is normal. For error distributions with ν 3 greater than zero, the limit of (r T H −1 r) is larger than the Fisher information, i.e. the asymptotic variance becomes even smaller. This holds for irregular cases such as the uniform distribution on [−1, 1]. The above shows that the estimator ϑ is close to being asymptotically efficient if the number of quantiles is large.

Simulations
In this section we study the numerical performance of our method. We consider samples of size n = 200 throughout. We draw random vectors X 1 , . . . , X n independently from a p-dimensional multivariate normal distribution with mean zero and a variance-covariance matrix whose (i, j) th component is 0.5 |i−j| . The full model size is p = 12 or p = 400 and the nonzero parameters are (ϑ 1 , ϑ 2 , ϑ 5 ) = (3.0, 1.5, 2.0). The response vector is Y = X T ϑ + ε.
In the objective functions (2) and (4) for the preliminary and the final estimator we set K = K ′ = 9 and use the same loss function, namely the check function ℓ k (x) = ℓ ′ k (x) = (x−β k ){τ k −I(x < β k )} with τ k = k/10, k = 1, . . . , 9. We call this method "weighted quantile regression through vote" (WQR-vote). To guarantee the best performance we use an iterative scheme for the preliminary estimator ϑ k . The initial value is ϑ |θ j |} and the updates in the t th iteration are Here λ k d k (x) = λ k I(|x| ≤ λ k ) + (b − 1) −1 (bλ k − |x|) + I(|x| > λ k ) is the derivative of the SCAD penalty with b being a constant that is usally set to 3.7 (Fan & Li, 2001). We repeat Table 1: Mean numbers of correctly selected nonzero parameters (MNC) and mean numbers of incorrectly selected zero parameters (MNI) of the weighted quantile regression through vote method (WQR-vote), least absolute deviation regression (LADR), least squares regression (LSR) and composite quantile regression (CQR), and the relative efficiency (RE) of the WQR-Vote to the three competing methods for various error distributions. Higher values of relative efficiency indicate better performance of the WQR-vote in estimation. The full model size is p = 12; The abbreviation SMN stands for a scale mixture of normal distributions 0.5N(0, 6)+0.5N(0, 6×0.5 6 ); LMN is a location mixture of normal distributions 0.5N(−2, 1)+0.5N(2, 1) and DE denotes the double exponential distribution with mean 0 and variance 2.  (Zou & Li, 2008).
To solve (8) we use the quick iterative coordinate descent algorithm (Peng & Wang, 2015). We minimize the multiple objective functions in (2) and (4) in parallel by using the foreach function from the R package foreach. The tuning parameters in the penalty term and the threshold α for the vote procedure (3) are determined by minimizing prediction errors that are calculated over an independent validation set {( X i , Y i ) : i = 1, . . . , n ′ } with n ′ = 2, 000, which is generated from the distribution of (X 1 , Y 1 ). We choose the tuning parameter λ k in (8) by setting λ k = arg min a∈L where L is a fine grid and ϑ a,k is the outcome of the iterative process (8) with λ k = a. The candidate set for α is A = {⌈K/2⌉, ⌈K/2⌉ + 1, . . . , K − 1} and the criterion is α = arg min a∈A Here the parameter estimator is ϑ a,k = arg min θ∈ Θn(a) n i=1 ℓ k (Y i − X T i θ) with Θ n (a) = {θ ∈ R pn : θ j = 0 for j / ∈ Q n (a)}. The vector ξ = ( ξ 1 , . . . , ξ K ) T is a consistent estimator of the optimal weights ξ = R −1 ψ, where R is a K × K matrix with (i, j) th entry Bradic et al., 2011, Section 4.2). Recall that β k is the τ k quantile of the error distribution for k = 1, . . . , K. After obtaining the set Q(α), we construct the estimator ϑ in (5) as described in Remark 1 in Section 3. Remark 1 also explains how to estimate f (β k ) for k = 1, . . . , K.
For comparison we consider least absolute deviation regression, least squares regression and composite quantile regression with the same nine quantile levels 1/10, . . . , 9/10. The iterative penalty in (8) is also applied to the three competing methods. For each of the four estimators, three indices from 200 simulated data sets are recorded in the following tables: 1. mean number of correctly selected nonzero parameters; 2. mean number of incorrectly selected zero parameters; where ϑ oracle is the oracle version of the four respective estimators and ϑ WQR−vote is the estimator from the WQR-vote. Here "oracle" means knowing the index set Q = {1, 2, 5} of the nonzero parameters before estimating and applying no penalty.
In Table 1 (Bradic et al., 2011, Sections 3.4 and 4.2), which takes extra time. In Table 2 we present simulation results as in Table 1, but now for the high dimensional scenario with p = 400 predictors. When dealing with high dimensional data, the WQR-vote still yields the smallest numbers of incorrect selections for all of the cases and the lowest mean squared errors of estimation for most of the cases. In the 200 simulations with normal errors the average time required for the WQR-vote and the composite quantile regression is 27.12 seconds and 106.41 seconds, respectively, when both methods use the nine quantile levels. The advantage of adopting parallel computing in the WQR-vote procedure is obvious.

A data application
In this section we analyze a microarray dataset from Scheetz et al. (2006), which is available at https://www.ncbi.nlm.nih.gov/geo/geo2r/?acc=GSE5680. The dataset contains the gene expression values of 31,042 probes on 120 rats. We are interested in how the expression of TRIM32, which is related to human hereditary diseases of the retina and corresponding to probe 1389163 at, depends on that of other genes. To exclude genes without sufficient variability, we first remove probes whose maximum among the 120 rats is less than the 25 th percentile of all the expression values or whose range among the 120 rats is less than 1. We then sort the remaining 18,984 probes by the absolute values of their correlation coefficients with the response, 1389163 at. The top 300 probes are used as predictors in the anlaysis.
To identify important predictors by the vote procedure (3), we set K = 9 and ℓ k (x) = x{k/10 − I(x < β k )} in the objective function (2), where β k is the k/10 quantile of the error distribution and can be estimated as an extra parameter. Then three different estimators are computed based on one loss function, i.e. K ′ = 1 in (4): the square function, the absolute value function and the composite check function ℓ ′ 1 (x) = 9 m=1 x{m/10 − I(x < β m )}. For comparison, the data are also analyzed by penalized regression using each of the above three loss functions. The iterative penalty (8) is applied to all methods.
We first use the entire dataset to fit models. Then the dataset is randomly divided into a training set of 80 observations and a validation set of 40 observations. Models are fitted with the training set and prediction errors are calculated by the loss functions on the validation set. Based on five-fold cross validation, we choose the tuning parameters in the penalty and the threshold α in (3) by the criteria stated in Section 4. This is repeated 50 times. Table 3 gives the sizes of the models that are selected by the vote procedure and by the three competing penalized regression approaches. The left panel provides the sizes that are obtained using all data. The first column of the right panel lists the means and, in parentheses, the standard deviations of the model sizes that are obtained using the 50 randomly generated partitions. The top half of the second column provides the prediction errors of the parameter estimates that are obtained in the second step of our method, i.e. after model selection by the vote procedure, using least absolute deviation regression, least Table 3: Sizes of selected subsets based on all data and on the 50 randomly generated partitions. We consider least absolute deviation regression (LADR), least squares regression (LSR) and composite quantile regression (CQR). The term LADR-Vote denotes the LADR appoach (without penalty) after variable selection through the vote procedure; LSR-Vote and CQR-Vote are defined analogously. Among all methods, the vote method generates the smallest models in the situations of both all data and random partition. In addition to achieving more sparsity, the submodel selected by the vote method yields a smaller mean prediction error compared with that selected by the penalized regression approach, when the same loss function is used to estimate the parameters. This confirms the superiority of the vote method in variable selection. Furthermore, this analysis illustrates the flexibility of our approach from employing different loss functions for selection and estimation. One may also apply the method optimally combining multiple quantile levels described in Section 3 to the data. We do not consider it here because the competing methods use different loss functions for estimation. The prediction errors are therefore not directly comparable.

Conclusion and discussion
We have developed a process that uses multiple loss functions separately to analyze linear models in the high dimensional data situation. Our approach has three notable advantages: (a) it lowers the false discovery rate of variable selection by using our newly-developed vote procedure; (b) it improves the quality of parameter estimation by combining results from multiple loss functions; (c) it increases the speed of computation by adopting parallel computing.
A specific instance of our approach, which optimally combines multiple quantile levels, achieves asymptotic efficiency of parameter estimation under mild conditions.
In practical applications some loss functions may be more relevant for model selection than others. Therefore a weighted version of the variable selector (3), namely with a nonnegative weight vector w = (w 1 , . . . , w K ) T , may improve our vote procedure, which uses uniform weights w k = 1 for k = 1, . . . , K. Specifying and estimating the weight vector w may be data-driven and vary from case to case. However, using uniform weights makes the vote procedure easier to implement. According to the results of Sections 4 and 5, this suffices to improve on competing methods with respect to variable selection.
The equations (12) and (13) gurantee the conditions (27) and (28) of Lemma 1 in Bradic et al. (2011) are satisfied, which implies that ϑ o k is the unique global minimizer of the objective function in (2) on Γ n,k . This combined with (10), the definition of ϑ k in (2) and Assumption 6 implies This gives the first conclusion of the theorem. The equation (14) and the definition of ϑ o k in (6) further yield Moreover, we know that, with probability tending to one, In the above the first two steps use (15) and (14) and the last step follows from (11). Then we have with ϑ kQn = ( ϑ k1 , . . . , ϑ kqn ) T , where the convergence follows from (16) and Assumption 8.
Let w * be a consistent estimator of w * . Then The last step uses the fact that |v T (X T X) 1/2 ( ϑ kQn − ϑ Qn )| = O p (1) for 1 ≤ k ≤ K, which holds by Theorem 2, and the consistency of w * . Combining (23) and (24) yields Proof of Theorem 3: The i th diagonal element of the K × K matrix H −1 is 2(K + 1){f (β i )} 2 for i = 1, . . . , K. The (i, i + 1) th and (i + 1, i) th entries are −(K + 1)f (β i )f (β i+1 ) for i = 1, . . . , K − 1, and the other entries are zero. Hence we have where The sixth step in the above uses the Cauchy-Schwarz inequality and the last step follows from Assumption 9. Then, using (25), we obtain The last step follows from Assumption 9, the equation (26) and the fact that τ 1 = (K + 1) −1 → 0 as K → ∞. This completes the proof.