M-estimation in high-dimensional linear model

We mainly study the M-estimation method for the high-dimensional linear regression model and discuss the properties of the M-estimator when the penalty term is a local linear approximation. In fact, the M-estimation method is a framework which covers the methods of the least absolute deviation, the quantile regression, the least squares regression and the Huber regression. We show that the proposed estimator possesses the good properties by applying certain assumptions. In the part of the numerical simulation, we select the appropriate algorithm to show the good robustness of this method.


Introduction
For the classical linear regression model Y = Xβ + ε, we are interested in the problem of variable selection and estimation, where Y = (y 1 , y 2 , . . . , y n ) T is the response vector, X = (X 1 , X 2 , . . . , X p n ) = (x 1 , x 2 , . . . , x n ) T = (x ij ) n×p n is an n × p n design matrix, and ε = (ε 1 , ε 2 , . . . , ε n ) T is a random vector. The main topic is how to estimate the coefficients vector β ∈ R p n when p n increases with sample size n and many elements of β equal zero. We can transfer this problem into a minimization of a penalized least squares objective functionβ where · is the l 2 norm of the vector, λ n is a tuning parameter, and p λ n (|t|) a penalty term. It is well known that the least squares estimation is not robust, especially when in the data there exist abnormal values or the error term has a heavy-tailed distribution.
In this paper we consider the loss function to be the least absolute deviation, i.e., we minimize the following objective function: where the loss function is the least absolute deviation (LAD for short) that does not need the noise to obey a gaussian distribution and be more robust than a least squares estimation. In fact, the LAD estimation is a special case of the M-estimation, which was mentioned by Huber (1964Huber ( , 1973Huber ( , 1981 [1][2][3] firstly and which can be obtained by minimizing the objective function where the function ρ can be selected. For example, if we choose ρ(x) = 1 2 x 2 1 |x|≤c + (c|x|c 2 /2)1 |x|>c , where c > 0, the Huber estimator can be obtained; if we choose ρ(x) = |x| q , where 1 ≤ q ≤ 2, L q estimator will be obtained, with two special cases: LAD estimator for q = 1 and OLS estimator for q = 2. If we choose ρ(x) = αx + + (1α)(-x) + , where 0 < α < 1, x + = max(x, 0), we call it a quantile regression, and we can also get the LAD estimator for α = 1/2 especially.
When p n approaches infinity as n tends to infinity, we assume that the function ρ is convex and not monotone, and the monotone function ϕ is the derivative of ρ. By imposing the appropriate regularity conditions, Huber (1973), Portnoy (1984) [4], Welsh (1989) [5] and Mammen (1989) [6] have proved that the M-estimator enjoyed the properties of consistency and asymptotic normality, where Welsh (1989) gave the weaker condition imposed on ϕ and the stronger condition on p n /n. Bai and Wu [7] further pointed out that the condition on p n could be a part of the integrable condition imposed on the design matrix. Moreover, He and Shao (2000) [8] studied the asymptotic properties of the M-estimator in the case of a generalized model setting and the dimension p n getting bigger and bigger. Li (2011) [9] obtained the Oracle property of the non-concave penalized M-estimator in high-dimensional model with the condition of p n log n/n → 0, p 2 n /n → 0, and proposed RSIS to make a variable selection by applying a rank sure independence screening method in the ultra high-dimensional model. Zou and Li (2008) [10] combined a penalized function and a local linear approximation method (LLA) to prove that the obtained estimator enjoyed good asymptotic properties, and they demonstrated that this method improved the computational efficiency of a local quadratic approximation (LQA) in a simulation.
Inspired by this, in this paper we consider the following problem: where p λ n (·) is the derivative of the penalized function, andβ n = (β n1 ,β n2 , . . . ,β np n ) T is the non-penalized estimator.
In this paper, we assume that the function ρ is convex, hence the objective function is still convex and the obtained local minimizer is a global minimizer.
Next, we state some assumptions which will be needed in the following results.
(A 1 ) The function ρ is convex on R, and its left derivative and right derivative ϕ The error term ε is i.i.d, and the distribution function F of ε i satisfies F(S) = 0, where S is the set of discontinuous points of ϕ.
Let z i be the transpose of the ith row vector of X (1) , such that lim n→∞ n -1 2 × max 1≤i≤n z T i z i = 0. It is worth mentioning that conditions (A 1 ) and (A 2 ) are classical assumptions for an M-estimation in a linear model, which can be found in many references, for example Bai, Rao and Wu (1992) [11] and Wu (2007) [12]. The condition (A 3 ) is frequently used for a sparse model in the linear model regression theory, which requires that the eigenvalues of the matrices D and D 11 are bounded. The condition (A 4 ) is weaker than that in previous references. Under the condition (A 4 ) we broaden the order of p n to n 1/2 , but Huber (1973) and Li, Peng and Zhu (2011) [9] required that p 2 n /n → 0, Portnoy (1984) required p n log p n /n → 0, and Mammen (1989) required p 3/2 n log p n /n → 0. Compared with these results, it is obvious that our sparse condition is much weaker. The condition (A 5 ) is the same as that in Huang, Horowitz and Ma (2008) [13], which is used to prove the asymptotic properties of the nonzero part of M-estimation.  Remark 2.2 By Theorem 2.2, we see that under the suitable conditions the global Mestimation of the zero-coefficient variables goes to zero with a high probability when n is large enough. This also shows that the model is sparse.

Theorem 2.3 (Oracle property) If the conditions (A 1 )-(A 5 ) hold and λ min (D) >
with probability converging to one, the non-concave penalized M-estimation β n = (β T n (1) ,β T n (2) ) T has the following properties: (1) (The consistency of the model selection)β n(2) = 0; (2) (Asymptotic normality) where s 2 n = σ 2 γ -1 u T D -1 11 u, and u is any k n dimensional vector such that u ≤ 1. Meanwhile, z i is the transpose of the ith row vector of a k n × k n matrix X (1) . Remark 2.3 From Theorem 2.3, the M-estimation enjoys the Oracle property, that is, the M-estimator can correctly select covariates with nonzero coefficients with probability converging to one and the estimators of the nonzero coefficients has the same asymptotic distribution that they would have if the zero coefficients were known in advance. (2004) [14], the authors showed that the non-concave penalized M-estimation has the property of consistency with the condition p 4 n /n → 0, and enjoyed the property of asymptotic normality with the condition p 5 n /n → 0. By Theorems 2.1-2.3, we can see that the corresponding conditions we impose are quite weak.

Proofs of main results
The proof of Theorem 2.1 Let α n = (p n /n) 1/2 + p 1/2 n c n . For any p n -dimensional vector u with u = C, we only need to prove that there exists a great enough positive constant C such that for any ε > 0, that is, there at least exists a local minimizerβ n such that β n - Firstly, by the triangle inequality we get Combining with the Von-Bahr-Esseen inequality and the fact that We can easily obtain E(T 111 ) = 0. From the Von-Bahr-Esseen inequality, the Schwarz inequality and the condition (B 3 ), it follows that so together with the Markov inequality this yields As for T 112 , Finally, considering T 2 , we can easily obtain This together with (3.1)-(3.5) shows that we can choose a great enough constant C such that T 111 and T 2 are controlled by T 112 , from which it follows that there at least exists a local minimizerβ n such that β nβ 0 = O P (α n ) in the closed ball {β 0 + α n u : u ≤ C}.
The proof of Theorem 2.3 It is obvious that the conclusion (1) can be obtained instantly by Theorem 2.2, so we only need to prove the conclusion (2). It follows from Theorem 2.1 that β n is consistent with β 0 andβ n(2) = 0 with probability converging to one from Theorem 2.2. Therefore forβ n(1) In the following part we give the Taylor expansion of upper left first term: Noticing that which shows that Then, as long as u ≤ 1, holds, for any k n -dimensional vector u. For the upper right third term, we can obtain Now let us deal with the upper right second term. Theorem 2.1 and the condition (A 3 ) yield where the upper third inequality sign holds because of Lemma 3 of Mammen (1989). Combining (3.9)-(3.10), we have that is, Denote s 2 n = σ 2 γ -1 u T D -1 11 u, F in = n -1/2 s -1 n γ -1 u T D -1 11 z T i , where z i is a k n × k n matrix and the transpose of the ith row vector of X (1) , then n 1/2 u T (β n(1)β 0(1) ) = n i=1 F in ϕ(ε i ) + o P (1). It follows from (A 5 ) that Applying the Slutsky theorem, we see that

Simulation results
In this section we evaluate the performance of the M-estimator proposed in (1.1) by simulation studies.
Then for the loss function. In this section we can choose some special loss functions, such as the LAD loss function, the OLS loss function and the Huber loss function. In this paper we choose the LAD loss function and the Huber loss function.
About stimulation algorithm. For the proposed LLA method, we connect the penalty function with independent variables and an independent variable, respectively, then we write a program by using the quantile package in R. For the Lasso method, we use the Lars package to simulate. Now we address the selection of the tuning parameter. We apply the BIC criterion to select the tuning parameter. The criterion is where DF λ n is the generalized degree of freedom used by Fan and Li (2001). About the selection of the evaluation index. In order to evaluate the performance of the estimators, we select four measures called EE, PE, C, IC and CP, which are obtained by 500 replicates. EE is the median of ββ 0 2 to evaluate the estimation accuracy, and PE is the prediction error defined by the median of n -1 Y -Xβ 2 . The other three measures are to qualify the performance of model consistency, where C and IC refer to the average number of correctly selected zero covariates and the average number of incorrectly selected zero covariates, and CP is the proportion of the number of the correct selection of zero variables to the total number of zero variables.
In the following we will compare the performances of the method LLA we proposed, the Lasso method and the Oracle estimation. Set n = 200, 500, 700, respectively, and p = [2 √ n]. From Table 1, we notice that the indices EE, C, IC, CP of our proposed LLA method perform better when ε ∼ N(0, 1). In particular, for the index CP, LLA outperforms Lasso. The reason for this may be that we impose different penalties for important and unimportant variables, while Lasso imposes the same penalties for all variables. Moreover, with the increase of the sample size, the ability of LLA method to correctly identify unimportant variables is also increasing. When the sample size is 700 and the number of explanatory variables is 53, an average of 48.9617 unimportant variables-zero variables are estimated to be zero on average, with an average accuracy of 99.92%.
An interesting fact can be found from Table 2, that is, when the error term is chosen as t 5 , the accuracy of the method LLA proposed to correctly exclude incorrect variables   is slightly higher than that of the case where the error term is a standardized normal distribution. The reason is that when the error term is heavy tailed, it is more appropriate to choose LLA, but the accuracy of estimation and prediction is slightly worse than that of Lasso. When the sample size increases, the LLA and Oracle estimates perform equally well in the selection of important variables and the complexity of the model. As can be seen from Table 3, when the error term is set to a mixed normal distribution, the ability of the proposed method to correctly select zero variables is good. In the case of a small sample size, the ability of the Lasso method to select important variables is better.
From Tables 4-6 where we choose the Huber loss function, the LLA method we proposed behaves well both in variable selection and robustness. Compared with Table 1 and  Table 4, when the data has outliers, we should choose LAD as the loss function. Moreover, when the error term follows a mixed normally distribution, the LLA method behaves bet-   ter than the Lasso method. The reason for this is that the real data has a mixed normal distribution with high probability.

Conclusion
In this paper, we mainly study the M-estimation method for the high-dimensional linear regression model and discuss the properties of the M-estimator when the penalty term is the local linear approximation. We show that the proposed estimator possesses the good properties by applying certain assumptions. In the numerical simulation, we select the appropriate algorithm to show the good robustness of this method.