# Variable Selection and Parameter Estimation with the Atan Regularization Method

Authors: Yanxin Wang; Li Zhu
Journal: Journal of Probability and Statistics
 (2016)

---

## Abstract

Variable selection is fundamental to high-dimensional statistical modeling. Many variable selection techniques may be implemented by penalized least squares using various penalty functions. In this paper, an arctangent type penalty which very closely resemblesl
0 penalty is proposed; we call it Atan penalty. The Atan-penalized least squares procedure is shown to consistently select the correct model and is asymptotically normal, provided the number of variables grows slower than the number of observations. The Atan procedure is efficiently implemented using an iteratively reweighted Lasso algorithm. Simulation results and data example show that the Atan procedure with BIC-type criterion performs very well in a variety of settings.


---

## Body

## 1. Introduction

High-dimensional data arise frequently in modern applications in biology, economics, chemometrics, neuroscience, and other scientific fields. To facilitate the analysis, it is often reasonable and useful to assume that only a small number of covariates are relevant for modeling the response variable. Under this sparsity assumption, a widely used approach for analyzing high-dimensional data is regularized or penalized regression. This approach estimates the unknown regression coefficients by solving the following penalized regression problem:(1)
m
i
n
β
∈
R
p
1
2
n
y
-
X
β
2
+
∑
j
=
1
p
p
λ
β
j
,where X
=
(
x
1
,
x
2
,
…
,
x
n
)
T is n
×
p design matrix, y
=
(
y
1
,
y
2
,
…
,
y
n
)
T is an n-dimensional response vector, β
=
(
β
1
,
…
,
β
p
)
T is the vector of unknown regression coefficients, · denotes L
2 norm (Euclidean norm), and p
λ
· is a penalty function which depends on a tuning parameter λ
>
0.In the above regularization framework, various penalty functions are used to perform variable selection by putting relatively large penalties on small coefficients, such as the best subset selection,l
1 penalized regression or Lasso [1], Bridge regression [2], SCAD [3], MCP [4], SICA [5], SELO [6], Dantzig selector [7], and Bayesian variable selection method [8]. The mainstream methods are the Lasso, Dantzig selector, and the folded concave penalization [9] such as the SCAD and MCP. The best subset selection, namely, l
0 penalty, along with the traditional model selection criteria such as AIC, BIC, and RIC [10–12] is attractive for variable selection since it directly penalizes the number of nonzero coefficients. However, one drawback of l
0 penalized least squares (PLS) procedure is instability of the resulting estimators [13]. This results from the fact that l
0 penalty is not continuous at 0. Another perhaps more significant drawback of l
0 penalty is that implementing l
0 PLS procedures is NP-hard and may involve an exhaustive search over all possible models. Thus, implementing these procedures is computationally infeasible when the number of potential predictors is even moderately large, let along the high-dimensional data.The Lasso penalized regression is computationally attractive and enjoys great performance in prediction. However, Lasso may not consistently select the correct model and is not necessarily asymptotically normal [14, 15]; a strong irrepresentable condition is necessary for the Lasso to be selection consistent [15, 16]. The folded concave penalization, unlike the Lasso, does not require the irrepresentable condition to achieve the variable selection consistency and can correct the intrinsic estimation bias of the Lasso. Fan and Li [3] first systematically studied nonconvex penalized likelihood for fixed finite dimension p. They recommended the SCAD penalty which enjoys the oracle property (a variable selection and estimation procedure is said to have the oracle property if it selects the true model A, with probability tending to one, and if the estimated coefficients are asymptotically normal, with the same asymptotic variance as the least squares estimator based on the true model) for variable selection. Fan and Peng [17] extended these results by allowing p to grow with n at the rate p
=
o
n
1
/
5 or p
=
o
n
1
/
3. Lv and Fan [5] introduced the weak oracle property, which means that the estimator enjoys the same sparsity as the oracle estimator with asymptotic probability one and has consistency, and established regularity conditions under which the PLS estimator given by folded concave penalties has a nonasymptotic weak oracle property when dimensionality p can grow nonpolynomially with sample size n. Theoretical properties enjoyed by SELO [6] estimators allow the number of predictors p to tend to infinity, along with the number of observations n, provided p
/
n
→
0. For high-dimensional nonconvex penalized regression with p
>
n, Kim et al. [18] proved that the oracle estimator itself is a local minimum of SCAD penalized least squares regression under very relaxed conditions; Zhang [4] proposed MCP and devised a novel PLUS algorithm which when used together can achieve the oracle property under certain regularity conditions. Recently, Fan et al. [9] have shown that the folded concave penalization methods enjoy the strong oracle property for high-dimensional sparse estimation. Important insight has also been gained through the recent work on theoretical analysis of the global solution [19–21].The practical performance of PLS procedures depends heavily on the choice of a tuning parameter. The theoretically optimal tuning parameter does not have an explicit representation and depends on unknown factors such as the variance of the unobserved random noise. Cross-validation is commonly adopted in practice to select the tuning parameter but is observed to often result in overfitting. In the case of fixedp, Wang et al. [22] proposed that one selects tuning parameter by minimizing the generalized BIC tuning parameter selector. Wang et al. [23] extended those results to the setting of a diverging number of parameters. Recently, Dicker et al. [6] proposed a BIC-like tuning parameter selector. Wang et al. [21] extended the work of [24, 25] for BIC on high-dimensional least squares regression; they proposed a high-dimensional BIC for a nonconvex penalized solution path.In this paper, we propose an arctangent type (Atan) penalty function which very closely approximatesl
0 penalty. Because the Atan penalty is continuous, the Atan estimator may be more stable than the estimators obtained through l
0 penalized methods. The Atan penalty is a smooth function on [
0
,
∞
) and we use an iteratively reweighted Lasso algorithm. We formally establish the model selection oracle property enjoyed by the Atan estimator. In particular, the asymptotic normality of the Atan is formally established. Our asymptotic framework allows the number of predictors, p
→
∞, along with the number of observations n, provided p
/
n
→
0. Furthermore, a BIC-like tuning parameter selection procedure is implemented for Atan.This paper is organized in the following way. In Section2, we introduce PLS estimators and give a brief overview of existing nonconvex penalty terms, and the Atan penalty is then presented. Then, we discuss some of its theoretical properties in Section 3. In Section 4, we describe a simple and efficient algorithm for obtaining Atan estimator. Simulation studies and an application of the proposed methodology are presented in Section 5. Conclusions are given in Section 6. The proofs are relegated to the Appendix.


## 2. The Atan-Penalized Least Squares Method

### 2.1. Linear Model and Penalized Least Squares

Suppose that{
(
y
i
,
x
i
)
}
i
=
1
n is a random sample from the linear regression model(2)
y
=
X
β
∗
+
ε
,where X
=
(
x
1
,
x
2
,
…
,
x
n
)
T is n
×
p design matrix, y
=
(
y
1
,
y
2
,
…
,
y
n
)
T is an n-dimensional response vector, and ε are the iid random errors with mean 0 and variance σ
2
n-dimensional noise vector.When discussing variable selection, it is convenient to have concise notation. Denote the columns ofX by x
1
,
…
,
x
p
∈
R
n and the rows of X by x
1
,
…
,
x
n
∈
R
p. Let A
=
{
j
;
β
j
∗
≠
0
} be the true model and suppose that p
0 is the size of the true model. That is, suppose that A
=
p
0, where A denotes the cardinality of A. In addition, for S
⊆
{
1,2
,
…
,
p
}, let β
S
=
(
β
j
)
j
∈
S be the S-dimensional subvector of β containing entries indexed by S and let X
S be n
×
S matrix obtained from X by extracting columns corresponding to S. Given p
×
p matrix C and subsets S
1
,
S
2
⊆
1,2
,
…
,
p, let C
S
1
,
S
2 be S
1
×
S
2 submatrix of C with rows determined by S
1 and columns determined by S
2.Various penalty functions have been used in the variable selection literature for linear regression model. Commonly used penalty functions includel
q, 0
≤
q
≤
2, the nonnegative garrotte [26], elastic-net [27, 28], SCAD [3], and MCP [4]. In particular, l
1 penalized least squares procedure is called the Lasso. However, Lasso estimates may be biased and inconsistent for model selection [3, 15]. This implies that the Lasso does not have the oracle property. The adaptive Lasso is a weighted version of Lasso which has the oracle property [15]. Slightly abusing notation is that the adaptive Lasso penalty is defined by p
λ
(
β
)
=
λ
ω
j
β, where ω
j is a data-dependent weight.Fan and Li [3] proposed a continuously differentiable penalty function called the SCAD penalty, which is defined by(3)
p
λ
′
β
=
λ
I
β
≤
λ
+
a
λ
-
β
+
a
-
1
λ
I
β
>
λ
,
for some
a
>
2
.Authors of [3] suggested using a
=
3.7 from a Bayesian perspective. The minimum concave penalty [4] translates the flat part of the derivative of the SCAD into the origin and is given by(4)
p
λ
′
β
=
a
λ
-
β
+
a
,which minimizes the maximum of the concavity. Zhang [4] proved that the MCP procedure may select the correct model with probability tending to 1 and that MCP estimator has good properties in terms of l
p-loss, provided λ and a satisfy certain conditions. Zhang’s results in fact allow for p
≫
n.


### 2.2. Atan Penalty

In this section we first propose a novel nonconvex penalty that we call Atan penalty. We then study its applications in sparse modeling.The Atan penalty is defined by(5)
p
λ
,
γ
β
=
λ
γ
+
2
π
arctan
⁡
β
γ
,for λ
≥
0 and γ
>
0. It is clear that this penalty is concave in β. Moreover, we can establish its relationship with l
0 and l
1 penalties. In particular, we have the following propositions.Proposition 1.
Letp
λ
,
γ
β be given in (5); then (6)
a
l
i
m
γ
→
∞
p
λ
,
γ
β
=
λ
β
,
b
l
i
m
γ
→
0
p
λ
,
γ
β
=
λ
,
i
f
β
≠
0
,
0
,
i
f
β
=
0
.The propositions show that the limits of Atan at 0 and ∞ are l
0 penalty and l
1 penalty, respectively. The first-order derivative of p
λ
,
γ
β with respect to β is (7)
p
λ
,
γ
′
β
=
λ
γ
γ
+
2
/
π
γ
2
+
β
2
.The Atan penalty function (γ
=
0.005) is plotted in Figure 1(a), along with the SCAD, Lasso, MCP, and l
0 penalty. Figure 1(b) depicts the Atan with different γ.Figure 1
(a) Plots of the penalties (includingl
0, Lasso, SCAD, MCP, and Atan). (b) Plots of the Atan penalties with different γ.
(a)
(b)


## 2.1. Linear Model and Penalized Least Squares

Suppose that{
(
y
i
,
x
i
)
}
i
=
1
n is a random sample from the linear regression model(2)
y
=
X
β
∗
+
ε
,where X
=
(
x
1
,
x
2
,
…
,
x
n
)
T is n
×
p design matrix, y
=
(
y
1
,
y
2
,
…
,
y
n
)
T is an n-dimensional response vector, and ε are the iid random errors with mean 0 and variance σ
2
n-dimensional noise vector.When discussing variable selection, it is convenient to have concise notation. Denote the columns ofX by x
1
,
…
,
x
p
∈
R
n and the rows of X by x
1
,
…
,
x
n
∈
R
p. Let A
=
{
j
;
β
j
∗
≠
0
} be the true model and suppose that p
0 is the size of the true model. That is, suppose that A
=
p
0, where A denotes the cardinality of A. In addition, for S
⊆
{
1,2
,
…
,
p
}, let β
S
=
(
β
j
)
j
∈
S be the S-dimensional subvector of β containing entries indexed by S and let X
S be n
×
S matrix obtained from X by extracting columns corresponding to S. Given p
×
p matrix C and subsets S
1
,
S
2
⊆
1,2
,
…
,
p, let C
S
1
,
S
2 be S
1
×
S
2 submatrix of C with rows determined by S
1 and columns determined by S
2.Various penalty functions have been used in the variable selection literature for linear regression model. Commonly used penalty functions includel
q, 0
≤
q
≤
2, the nonnegative garrotte [26], elastic-net [27, 28], SCAD [3], and MCP [4]. In particular, l
1 penalized least squares procedure is called the Lasso. However, Lasso estimates may be biased and inconsistent for model selection [3, 15]. This implies that the Lasso does not have the oracle property. The adaptive Lasso is a weighted version of Lasso which has the oracle property [15]. Slightly abusing notation is that the adaptive Lasso penalty is defined by p
λ
(
β
)
=
λ
ω
j
β, where ω
j is a data-dependent weight.Fan and Li [3] proposed a continuously differentiable penalty function called the SCAD penalty, which is defined by(3)
p
λ
′
β
=
λ
I
β
≤
λ
+
a
λ
-
β
+
a
-
1
λ
I
β
>
λ
,
for some
a
>
2
.Authors of [3] suggested using a
=
3.7 from a Bayesian perspective. The minimum concave penalty [4] translates the flat part of the derivative of the SCAD into the origin and is given by(4)
p
λ
′
β
=
a
λ
-
β
+
a
,which minimizes the maximum of the concavity. Zhang [4] proved that the MCP procedure may select the correct model with probability tending to 1 and that MCP estimator has good properties in terms of l
p-loss, provided λ and a satisfy certain conditions. Zhang’s results in fact allow for p
≫
n.


## 2.2. Atan Penalty

In this section we first propose a novel nonconvex penalty that we call Atan penalty. We then study its applications in sparse modeling.The Atan penalty is defined by(5)
p
λ
,
γ
β
=
λ
γ
+
2
π
arctan
⁡
β
γ
,for λ
≥
0 and γ
>
0. It is clear that this penalty is concave in β. Moreover, we can establish its relationship with l
0 and l
1 penalties. In particular, we have the following propositions.Proposition 1.
Letp
λ
,
γ
β be given in (5); then (6)
a
l
i
m
γ
→
∞
p
λ
,
γ
β
=
λ
β
,
b
l
i
m
γ
→
0
p
λ
,
γ
β
=
λ
,
i
f
β
≠
0
,
0
,
i
f
β
=
0
.The propositions show that the limits of Atan at 0 and ∞ are l
0 penalty and l
1 penalty, respectively. The first-order derivative of p
λ
,
γ
β with respect to β is (7)
p
λ
,
γ
′
β
=
λ
γ
γ
+
2
/
π
γ
2
+
β
2
.The Atan penalty function (γ
=
0.005) is plotted in Figure 1(a), along with the SCAD, Lasso, MCP, and l
0 penalty. Figure 1(b) depicts the Atan with different γ.Figure 1
(a) Plots of the penalties (includingl
0, Lasso, SCAD, MCP, and Atan). (b) Plots of the Atan penalties with different γ.
(a)
(b)


## 3. Theoretical Properties

In this section we study the theoretical properties of the Atan estimator proposed in Section2 in the situation where the number of parameters p tends to ∞ with increasing sample size n. We discuss some conditions of the penalty and loss functions in Section 3.1. Our main results are presented in Section 3.2.


### 3.1. Regularity Conditions

We need to place the following conditions on the penalty functions:(A)
n
→
∞ and p
σ
2
/
n
→
0.
(B)
ρ
n
/
(
p
σ
2
)
→
∞, where ρ
=
m
i
n
j
∈
A
β
j
∗.
(C)
λ
=
O
(
1
), λ
n
/
(
p
σ
2
)
→
∞, and γ
=
O
(
p
1
/
2
σ
3
n
-
3
/
2
).
(D)
There exist constantsC
1
,
C
2
∈
R such that C
1
<
λ
m
i
n
(
1
/
n
X
T
X
)
<
λ
m
a
x
(
1
/
n
X
T
X
)
<
C
2, where λ
m
i
n
(
1
/
n
X
T
X
) and λ
m
a
x
(
1
/
n
X
T
X
) are the smallest and largest eigenvalues of 1
/
n
X
T
X, respectively.
(E)
l
i
m
n
→
∞
n
-
1
m
a
x
1
≤
i
≤
n
∑
j
=
1
p
‍
x
i
j
2
=
0.
(F)
E
(
ε
i
/
σ
2
+
δ
)
<
M for some δ and M
<
∞.Sincep may vary with n, it is implicit that β
∗ may vary with n. Additionally, we allow model A and the distribution of ε (in particular, σ
2) to change with n. Condition (A) limits how p and σ
2 may grow with n. This condition is substantially weaker than that required in [17], which requires p
5
/
n
→
0, and slightly weaker than that required in [28], which requires log
⁡
p
/
log
⁡
n
→
ν
∈
[
0,1
) and the same as that required in [6]. As mentioned in Section 1, other authors have studied PLS methods in settings where p
>
n; that is, their growth condition on p is weaker than condition (A) [18]. Condition (B) gives a lower bound on the size of the smallest nonzero entry of β
∗. Notice that the smallest nonzero entry of β
∗ is allowed to vanish asymptotically, provided it does not do so faster than p
σ
2
/
n. Similar conditions are found in [17]. Condition (C) restricts the rates of tuning parameters λ and γ. Note that condition (C) does not constrain the minimum size of γ. Indeed, no such constraint is required for our asymptotic results about the Atan estimator. Since the Atan penalty approaches l
0 penalty as γ
→
0, this suggests that the Atan and l
0 penalized least squares estimator have similar asymptotic properties. On the other hand, in practice, we have found that one should not take γ too small, in order to preserve stability of the Atan estimator. Condition (D) is an identifiability condition. Conditions (E) and (F) are used to prove asymptotic normality of Atan estimators and are related to the Lindeberg condition, of the Lindeberg-Feller central limit theorem. Conditions (A)–(F) imply that Atan has the oracle property and may correctly identify model A as we will see in Theorem 4.


### 3.2. Oracle Properties

Let(8)
Q
n
β
=
1
2
n
y
-
X
β
2
+
∑
j
=
1
p
p
λ
,
γ
β
jbe the objective function, and p
λ
,
γ
β
j is the Atan penalty function.Theorem 2.
Suppose that conditions (A)–(D) hold; then, for everyr
∈
(
0,1
), there exists a constant C
0
>
0 such that(9)
l
i
m
inf
⁡
P
n
→
∞
arg
min
β
⁡
Q
n
β
⊆
β
∈
R
p
;
β
-
β
∗
≤
C
p
σ
2
n
>
1
-
r
,whenever C
≥
C
0. Consequently, there exists a sequence of local minimizers of Q
n
(
β
) and β
^, such that β
^
-
β
∗
=
O
P
(
p
σ
2
/
n
).Lemma 3.
Assume that (A)–(D) hold, and fixC
>
0; then (10)
l
i
m
n
→
∞
P
a
r
g
m
i
n
β
-
β
∗
≤
C
p
σ
2
/
n
Q
n
β
⊆
β
∈
R
p
;
β
A
c
=
0
=
1
,where A
c
=
1
,
…
,
p
∖
A is the complement of A in {
1
,
…
,
p
}.Theorem 4 (oracle properties).
Suppose that (A)–(F) hold; then there exists a sequence ofn
/
p
σ
2-consistent local minima of Atan, β
^, such that(i)
l
i
m
n
→
∞
P
j
;
β
^
j
≠
0
=
A
=
1,
(ii)
n
B
n
n
-
1
X
A
T
X
A
/
σ
2
1
/
2
β
^
A
-
β
A
∗
→
N
0
,
G,
in distribution, whereB
n is any arbitrary q
×
A matrix such that B
n
B
n
T
→
G, and G is q
×
q nonegative symmetric matrix.


## 3.1. Regularity Conditions

We need to place the following conditions on the penalty functions:(A)
n
→
∞ and p
σ
2
/
n
→
0.
(B)
ρ
n
/
(
p
σ
2
)
→
∞, where ρ
=
m
i
n
j
∈
A
β
j
∗.
(C)
λ
=
O
(
1
), λ
n
/
(
p
σ
2
)
→
∞, and γ
=
O
(
p
1
/
2
σ
3
n
-
3
/
2
).
(D)
There exist constantsC
1
,
C
2
∈
R such that C
1
<
λ
m
i
n
(
1
/
n
X
T
X
)
<
λ
m
a
x
(
1
/
n
X
T
X
)
<
C
2, where λ
m
i
n
(
1
/
n
X
T
X
) and λ
m
a
x
(
1
/
n
X
T
X
) are the smallest and largest eigenvalues of 1
/
n
X
T
X, respectively.
(E)
l
i
m
n
→
∞
n
-
1
m
a
x
1
≤
i
≤
n
∑
j
=
1
p
‍
x
i
j
2
=
0.
(F)
E
(
ε
i
/
σ
2
+
δ
)
<
M for some δ and M
<
∞.Sincep may vary with n, it is implicit that β
∗ may vary with n. Additionally, we allow model A and the distribution of ε (in particular, σ
2) to change with n. Condition (A) limits how p and σ
2 may grow with n. This condition is substantially weaker than that required in [17], which requires p
5
/
n
→
0, and slightly weaker than that required in [28], which requires log
⁡
p
/
log
⁡
n
→
ν
∈
[
0,1
) and the same as that required in [6]. As mentioned in Section 1, other authors have studied PLS methods in settings where p
>
n; that is, their growth condition on p is weaker than condition (A) [18]. Condition (B) gives a lower bound on the size of the smallest nonzero entry of β
∗. Notice that the smallest nonzero entry of β
∗ is allowed to vanish asymptotically, provided it does not do so faster than p
σ
2
/
n. Similar conditions are found in [17]. Condition (C) restricts the rates of tuning parameters λ and γ. Note that condition (C) does not constrain the minimum size of γ. Indeed, no such constraint is required for our asymptotic results about the Atan estimator. Since the Atan penalty approaches l
0 penalty as γ
→
0, this suggests that the Atan and l
0 penalized least squares estimator have similar asymptotic properties. On the other hand, in practice, we have found that one should not take γ too small, in order to preserve stability of the Atan estimator. Condition (D) is an identifiability condition. Conditions (E) and (F) are used to prove asymptotic normality of Atan estimators and are related to the Lindeberg condition, of the Lindeberg-Feller central limit theorem. Conditions (A)–(F) imply that Atan has the oracle property and may correctly identify model A as we will see in Theorem 4.


## 3.2. Oracle Properties

Let(8)
Q
n
β
=
1
2
n
y
-
X
β
2
+
∑
j
=
1
p
p
λ
,
γ
β
jbe the objective function, and p
λ
,
γ
β
j is the Atan penalty function.Theorem 2.
Suppose that conditions (A)–(D) hold; then, for everyr
∈
(
0,1
), there exists a constant C
0
>
0 such that(9)
l
i
m
inf
⁡
P
n
→
∞
arg
min
β
⁡
Q
n
β
⊆
β
∈
R
p
;
β
-
β
∗
≤
C
p
σ
2
n
>
1
-
r
,whenever C
≥
C
0. Consequently, there exists a sequence of local minimizers of Q
n
(
β
) and β
^, such that β
^
-
β
∗
=
O
P
(
p
σ
2
/
n
).Lemma 3.
Assume that (A)–(D) hold, and fixC
>
0; then (10)
l
i
m
n
→
∞
P
a
r
g
m
i
n
β
-
β
∗
≤
C
p
σ
2
/
n
Q
n
β
⊆
β
∈
R
p
;
β
A
c
=
0
=
1
,where A
c
=
1
,
…
,
p
∖
A is the complement of A in {
1
,
…
,
p
}.Theorem 4 (oracle properties).
Suppose that (A)–(F) hold; then there exists a sequence ofn
/
p
σ
2-consistent local minima of Atan, β
^, such that(i)
l
i
m
n
→
∞
P
j
;
β
^
j
≠
0
=
A
=
1,
(ii)
n
B
n
n
-
1
X
A
T
X
A
/
σ
2
1
/
2
β
^
A
-
β
A
∗
→
N
0
,
G,
in distribution, whereB
n is any arbitrary q
×
A matrix such that B
n
B
n
T
→
G, and G is q
×
q nonegative symmetric matrix.


## 4. Implementation

### 4.1. Iteratively Reweighted Lasso Algorithm

The computation for the Atan-penalized method is much more involved, because the resulting optimization problem is usually nonconvex and has multiple local minimizers. Several algorithms have been developed for computing the folded concave penalized estimators, such as the local quadratic approximation (LQA) algorithm and the local linear approximation (LLA) algorithm [3, 29]. Both LQA and LLA are related to the majorization-minorization (MM) principle [30]. Recently, coordinate descent algorithm was applied to solve the folded concave penalized least squares [31, 32]. Reference [4] devised a PLUS algorithm for solving the PLS using the MCP. Zhang [33] analyzed the capped-l
1 penalty for solving the PLS.In this section, we present an iteratively reweighted Lasso (IRL) algorithm to solve the following minimization problem:(11)
m
i
n
β
∈
R
p
1
2
n
y
-
X
β
2
+
λ
∑
j
=
1
p
γ
+
2
π
arctan
⁡
β
j
γ
.We show that the solution of the penalty least squares problem can be transformed into that of a series of convex weighted Lasso estimators, to which the existing Lasso algorithms can be efficiently applied.Now, taking a first-order Taylor-series approximation of the Atan penalty about the current valueβ
j
=
β
j
∗, we obtain the overapproximation of (11):(12)
m
i
n
β
∈
R
p
1
2
n
y
-
X
β
2
+
λ
∑
j
=
1
p
γ
γ
+
2
/
π
γ
2
+
β
j
∗
2
β
j
.Note that the linear approximation in (12) is analogous to one proposed in Zou and Li [29] and they argued strongly in favor of a one-step approximation. Instead, we offer the following algorithm for Atan-penalized least squares. Assume all the covariates have been standardized to have mean zero and unit variance. Without loss of generality, the span of regularization parameters is λ
m
i
n
=
λ
0
<
⋯
<
λ
L
=
λ
m
a
x, for l
=
0
,
…
,
L [34]. Here, we summarize the details of IRL algorithm as in Algorithm 1, where the active set is defined as A
=
j
∣
β
j
≠
0
,
j
=
1,2
,
…
,
p.Algorithm 1:Iteratively reweighted Lasso (IRL) algorithm.
(1)Start with λ
m
a
x
=
max
1
≤
j
≤
p
⁡
|
X
T
y
|
/
n. Set l
=
L and β
(
l
)
=
0.
Outer loop:
(2)Set k
=
0 and β
~
(
0
)
=
β
(
l
);
(3)Increment k
=
k
+
1 and β
~
(
k
+
1
)
=
β
~
(
k
);
(4)Update the weights: ω
~
j
=
γ
(
γ
+
2
/
π
)
/
(
γ
2
+
(
β
~
j
k
)
2
), j
=
1
,
…
,
p;
Inner loop:
Solve the Karush-Kuhn-Tucker (KKT) conditions for fixed ω
~
j:
x
j
T
(
y
-
X
β
)
-
λ
l
ω
~
j
sgn
⁡
(
β
j
)
=
0, if j
∈
A,
|
x
j
T
(
y
-
X
β
)
|
<
λ
l
ω
~
j, if j
∉
A,
(5)Goto Step (3);
(6)Repeat Steps (3)–(5) until ‖
β
~
(
k
+
1
)
-
β
~
(
k
)
‖
2
<
τ.
The estimate β
(
l
) is the limit point of the outer loop, β
~
∞;
(7)Decrement l
=
l
-
1 and λ
l
=
λ
l
-
1. Return to (2) using β
(
l
) as a warm start.In the situation withn
>
p, we set λ
m
i
n
=
0. However, when n
<
p, one attains the saturated model long before the regularization parameter reaches zero. In this case, we suggest λ
m
i
n
=
1
0
-
4
λ
m
a
x [34], and we used τ
=
1
0
-
4. In practice, we have found that if the columns of X are standardized so that x
j
2
=
n, for j
=
1
,
…
,
p, then taking γ
=
0.005 works well.Remark 5.
Algorithm1 is very much like MM algorithms. Hunter and Li [30] were the first to advocate MM algorithms for variable selection when n
>
p. However, we should notice that Algorithm 1 differs from the algorithm described in [30]. The use of MM algorithms in [30] is to justify a quadratic overapproximation to penalty functions with singularities at the origin. In the inner loop of Algorithm 1, we avoid such approximation by solving the KKT conditions precisely and efficiently [35]. Our use of MM algorithms in Algorithm 1 is to justify the local linear approximation to Atan penalty in the outer loop.Remark 6.
LLA to SCAD penalty were also proposed by Zou and Li [29]. Algorithm 1 differs from that proposed in [29] in that it constructs the entire coefficient path, even when p
>
n, whereas the procedure in [29] computes the coefficient estimates for fixed regularization parameter λ starting with a root-n consistent estimator of true coefficient β
∗ at the initial step. Moreover, even with the same initial estimate, the limit point from Algorithm 1 will differ from [29] after one iteration.Remark 7.
In the more general case wherep
>
n, we used coordinate-wise optimization to compute the entire regularized coefficient path via Algorithm 1. Our experience is that coordinate-wise optimization works well in practice and converges very quickly.


### 4.2. Regularity Parameter Selection

Tuning parameter selection is an important issue in most PLS procedures. There are relatively few studies on the choice of penalty parameters. Traditional model selection criteria, such as AIC [10] and BIC [11], suffer from a number of limitations. Their major drawback arises because parameter estimation and model selection are two different processes, which can result in instability [13] and complicated stochastic properties. To overcome the deficiency of traditional methods, Fan and Li [3] proposed the SCAD method, which estimates parameters while simultaneously selecting important variables. They selected tuning parameter by minimizing GCV [1, 3, 26].However, it is well known that GCV and AIC-based methods are not consistent for model selection in the sense that asn
→
∞, they may select irrelevant predictors with nonvanishing probability [22, 36]. On the other hand, BIC-based tuning parameter selection roughly corresponds to maximizing the posterior probability of selecting the true model in an appropriate Bayesian formulation and has been shown to be consistent for model selection in several settings [22, 37, 38]. The BIC tuning parameter selector is defined by (13)
B
I
C
=
log
⁡
y
-
X
β
^
2
n
+
D
F
^
log
⁡
n
n
,where D
F
^ is the estimate of the degrees of freedom given by (14)
D
F
^
=
t
r
X
X
T
+
n
Σ
λ
T
X
T
,and Σ
λ
=
diag
⁡
p
λ
′
β
^
1
/
β
^
1
,
…
,
p
λ
′
β
^
p
/
β
^
p. The diagonal elements of Σ
λ are coefficients of quadratic terms in the local quadratic approximation to SCAD penalty function p
λ
· [3].Dicker et al. [6] proposed BIC-like procedures implemented by minimizing(15)
B
I
C
0
=
log
⁡
y
-
X
β
^
2
n
-
p
^
0
+
log
⁡
n
n
p
^
0
.To estimate the residual variance, they use σ
^
2
=
(
n
-
p
^
0
)
-
1
y
-
X
β
^
2. This differs from other estimates of the residual variance used in PLS methods, where the denominator n
-
p
^
0 is replaced by n [22]; here, n
-
p
^
0 is used to account for degrees of freedom lost to estimation.More works on the high-dimensional BIC for the least squares regression to tuning parameter selection for nonconvex penalized regression can be seen in [24, 25, 39]. Here we used BIC statistic as (15).


### 4.3. A Standard Error Formula

The standard errors for the estimated parameters can be obtained directly because we are estimating parameters and selecting variables at the same time. Letβ
^
=
β
^
λ
,
γ be a local minimizer of Atan. Following [3, 17], standard errors of β
^ may be estimated by using quadratic approximations to Atan. Indeed, the approximation(16)
p
λ
,
a
β
j
≈
p
λ
,
γ
β
j
0
+
1
2
β
j
0
p
λ
,
γ
′
β
j
0
β
j
2
-
β
j
0
2
,
for
β
j
≈
β
j
0
.suggests that Atan may be replaced by the quadratic minimization problem(17)
m
i
n
1
n
y
-
X
β
2
+
∑
j
=
1
p
p
λ
,
a
′
β
j
0
β
j
0
β
j
2
,at least for the purposes of obtaining standard errors. Using this expression, we obtain a sandwich formula for the estimated standard error of β
^
A
^, where A
^
=
{
j
;
β
^
j
≠
0
}. Consider(18)
c
o
v
^
β
^
A
^
=
σ
^
2
X
A
^
T
X
A
^
+
n
Δ
A
^
,
A
^
β
^
-
1
X
A
^
T
X
A
^
X
A
^
T
X
A
^
+
n
Δ
A
^
,
A
^
β
^
-
1
,where Δ
(
β
)
=
d
i
a
g
p
λ
,
a
′
β
1
/
β
1
,
…
,
p
λ
,
a
′
β
p
/
β
p, σ
^
2
=
n
-
1
y
-
X
β
^
2, and p
^
0
=
A
^ is the number of elements in A
^. Under the conditions of Theorem 4, (19)
B
n
X
A
^
T
X
A
^
c
o
v
^
β
^
A
^
B
n
T
σ
2
⟶
G
.This is a consistency result for c
o
v
^
(
β
^
A
^
).


## 4.1. Iteratively Reweighted Lasso Algorithm

The computation for the Atan-penalized method is much more involved, because the resulting optimization problem is usually nonconvex and has multiple local minimizers. Several algorithms have been developed for computing the folded concave penalized estimators, such as the local quadratic approximation (LQA) algorithm and the local linear approximation (LLA) algorithm [3, 29]. Both LQA and LLA are related to the majorization-minorization (MM) principle [30]. Recently, coordinate descent algorithm was applied to solve the folded concave penalized least squares [31, 32]. Reference [4] devised a PLUS algorithm for solving the PLS using the MCP. Zhang [33] analyzed the capped-l
1 penalty for solving the PLS.In this section, we present an iteratively reweighted Lasso (IRL) algorithm to solve the following minimization problem:(11)
m
i
n
β
∈
R
p
1
2
n
y
-
X
β
2
+
λ
∑
j
=
1
p
γ
+
2
π
arctan
⁡
β
j
γ
.We show that the solution of the penalty least squares problem can be transformed into that of a series of convex weighted Lasso estimators, to which the existing Lasso algorithms can be efficiently applied.Now, taking a first-order Taylor-series approximation of the Atan penalty about the current valueβ
j
=
β
j
∗, we obtain the overapproximation of (11):(12)
m
i
n
β
∈
R
p
1
2
n
y
-
X
β
2
+
λ
∑
j
=
1
p
γ
γ
+
2
/
π
γ
2
+
β
j
∗
2
β
j
.Note that the linear approximation in (12) is analogous to one proposed in Zou and Li [29] and they argued strongly in favor of a one-step approximation. Instead, we offer the following algorithm for Atan-penalized least squares. Assume all the covariates have been standardized to have mean zero and unit variance. Without loss of generality, the span of regularization parameters is λ
m
i
n
=
λ
0
<
⋯
<
λ
L
=
λ
m
a
x, for l
=
0
,
…
,
L [34]. Here, we summarize the details of IRL algorithm as in Algorithm 1, where the active set is defined as A
=
j
∣
β
j
≠
0
,
j
=
1,2
,
…
,
p.Algorithm 1:Iteratively reweighted Lasso (IRL) algorithm.
(1)Start with λ
m
a
x
=
max
1
≤
j
≤
p
⁡
|
X
T
y
|
/
n. Set l
=
L and β
(
l
)
=
0.
Outer loop:
(2)Set k
=
0 and β
~
(
0
)
=
β
(
l
);
(3)Increment k
=
k
+
1 and β
~
(
k
+
1
)
=
β
~
(
k
);
(4)Update the weights: ω
~
j
=
γ
(
γ
+
2
/
π
)
/
(
γ
2
+
(
β
~
j
k
)
2
), j
=
1
,
…
,
p;
Inner loop:
Solve the Karush-Kuhn-Tucker (KKT) conditions for fixed ω
~
j:
x
j
T
(
y
-
X
β
)
-
λ
l
ω
~
j
sgn
⁡
(
β
j
)
=
0, if j
∈
A,
|
x
j
T
(
y
-
X
β
)
|
<
λ
l
ω
~
j, if j
∉
A,
(5)Goto Step (3);
(6)Repeat Steps (3)–(5) until ‖
β
~
(
k
+
1
)
-
β
~
(
k
)
‖
2
<
τ.
The estimate β
(
l
) is the limit point of the outer loop, β
~
∞;
(7)Decrement l
=
l
-
1 and λ
l
=
λ
l
-
1. Return to (2) using β
(
l
) as a warm start.In the situation withn
>
p, we set λ
m
i
n
=
0. However, when n
<
p, one attains the saturated model long before the regularization parameter reaches zero. In this case, we suggest λ
m
i
n
=
1
0
-
4
λ
m
a
x [34], and we used τ
=
1
0
-
4. In practice, we have found that if the columns of X are standardized so that x
j
2
=
n, for j
=
1
,
…
,
p, then taking γ
=
0.005 works well.Remark 5.
Algorithm1 is very much like MM algorithms. Hunter and Li [30] were the first to advocate MM algorithms for variable selection when n
>
p. However, we should notice that Algorithm 1 differs from the algorithm described in [30]. The use of MM algorithms in [30] is to justify a quadratic overapproximation to penalty functions with singularities at the origin. In the inner loop of Algorithm 1, we avoid such approximation by solving the KKT conditions precisely and efficiently [35]. Our use of MM algorithms in Algorithm 1 is to justify the local linear approximation to Atan penalty in the outer loop.Remark 6.
LLA to SCAD penalty were also proposed by Zou and Li [29]. Algorithm 1 differs from that proposed in [29] in that it constructs the entire coefficient path, even when p
>
n, whereas the procedure in [29] computes the coefficient estimates for fixed regularization parameter λ starting with a root-n consistent estimator of true coefficient β
∗ at the initial step. Moreover, even with the same initial estimate, the limit point from Algorithm 1 will differ from [29] after one iteration.Remark 7.
In the more general case wherep
>
n, we used coordinate-wise optimization to compute the entire regularized coefficient path via Algorithm 1. Our experience is that coordinate-wise optimization works well in practice and converges very quickly.


## 4.2. Regularity Parameter Selection

Tuning parameter selection is an important issue in most PLS procedures. There are relatively few studies on the choice of penalty parameters. Traditional model selection criteria, such as AIC [10] and BIC [11], suffer from a number of limitations. Their major drawback arises because parameter estimation and model selection are two different processes, which can result in instability [13] and complicated stochastic properties. To overcome the deficiency of traditional methods, Fan and Li [3] proposed the SCAD method, which estimates parameters while simultaneously selecting important variables. They selected tuning parameter by minimizing GCV [1, 3, 26].However, it is well known that GCV and AIC-based methods are not consistent for model selection in the sense that asn
→
∞, they may select irrelevant predictors with nonvanishing probability [22, 36]. On the other hand, BIC-based tuning parameter selection roughly corresponds to maximizing the posterior probability of selecting the true model in an appropriate Bayesian formulation and has been shown to be consistent for model selection in several settings [22, 37, 38]. The BIC tuning parameter selector is defined by (13)
B
I
C
=
log
⁡
y
-
X
β
^
2
n
+
D
F
^
log
⁡
n
n
,where D
F
^ is the estimate of the degrees of freedom given by (14)
D
F
^
=
t
r
X
X
T
+
n
Σ
λ
T
X
T
,and Σ
λ
=
diag
⁡
p
λ
′
β
^
1
/
β
^
1
,
…
,
p
λ
′
β
^
p
/
β
^
p. The diagonal elements of Σ
λ are coefficients of quadratic terms in the local quadratic approximation to SCAD penalty function p
λ
· [3].Dicker et al. [6] proposed BIC-like procedures implemented by minimizing(15)
B
I
C
0
=
log
⁡
y
-
X
β
^
2
n
-
p
^
0
+
log
⁡
n
n
p
^
0
.To estimate the residual variance, they use σ
^
2
=
(
n
-
p
^
0
)
-
1
y
-
X
β
^
2. This differs from other estimates of the residual variance used in PLS methods, where the denominator n
-
p
^
0 is replaced by n [22]; here, n
-
p
^
0 is used to account for degrees of freedom lost to estimation.More works on the high-dimensional BIC for the least squares regression to tuning parameter selection for nonconvex penalized regression can be seen in [24, 25, 39]. Here we used BIC statistic as (15).


## 4.3. A Standard Error Formula

The standard errors for the estimated parameters can be obtained directly because we are estimating parameters and selecting variables at the same time. Letβ
^
=
β
^
λ
,
γ be a local minimizer of Atan. Following [3, 17], standard errors of β
^ may be estimated by using quadratic approximations to Atan. Indeed, the approximation(16)
p
λ
,
a
β
j
≈
p
λ
,
γ
β
j
0
+
1
2
β
j
0
p
λ
,
γ
′
β
j
0
β
j
2
-
β
j
0
2
,
for
β
j
≈
β
j
0
.suggests that Atan may be replaced by the quadratic minimization problem(17)
m
i
n
1
n
y
-
X
β
2
+
∑
j
=
1
p
p
λ
,
a
′
β
j
0
β
j
0
β
j
2
,at least for the purposes of obtaining standard errors. Using this expression, we obtain a sandwich formula for the estimated standard error of β
^
A
^, where A
^
=
{
j
;
β
^
j
≠
0
}. Consider(18)
c
o
v
^
β
^
A
^
=
σ
^
2
X
A
^
T
X
A
^
+
n
Δ
A
^
,
A
^
β
^
-
1
X
A
^
T
X
A
^
X
A
^
T
X
A
^
+
n
Δ
A
^
,
A
^
β
^
-
1
,where Δ
(
β
)
=
d
i
a
g
p
λ
,
a
′
β
1
/
β
1
,
…
,
p
λ
,
a
′
β
p
/
β
p, σ
^
2
=
n
-
1
y
-
X
β
^
2, and p
^
0
=
A
^ is the number of elements in A
^. Under the conditions of Theorem 4, (19)
B
n
X
A
^
T
X
A
^
c
o
v
^
β
^
A
^
B
n
T
σ
2
⟶
G
.This is a consistency result for c
o
v
^
(
β
^
A
^
).


## 5. Numerical Studies

### 5.1. Simulation Studies

We now investigate the sparsity recovery and estimation properties of the proposed estimator via numerical simulations. We compare the following estimators: the Lasso estimator (implemented using R package glmnet [40]); the adaptive Lasso estimator (denoted by Alasso [15]), the SCAD estimator from the CD algorithm without calibration [31]; the MCP estimator from the CD algorithm with a
=
1.5 [31]; the Dantzig selector [7]; and Bayesian variable select method [8]. For the proposed Atan estimator, we take γ
=
0.005 and BIC statistic (15) is used to select the tuning parameter λ. In the following, we report simulation results from four examples.Example 8.
In this example, simulation data are generated from the linear regression model:(20)
y
=
x
T
β
∗
+
σ
ε
,where β
∗
=
(
3,1.5,0
,
0,2
,
0,0
,
0,0
,
0,0
,
0
)
T, ε
~
N
(
0,1
), and x is multivariate normal distribution with zero mean and covariance between the ith and jth elements being ρ
i
-
j with ρ
=
0.5. In our simulation, sample size n is set to be 100 and 200, σ
=
2. For each case, we repeated the simulation 200 times.For linear model, model error forμ
^
=
x
T
β
^ is M
E
(
μ
^
)
=
(
β
^
-
β
)
T
E
(
x
x
T
)
(
β
^
-
β
). Simulation results are summarized in Table 1, in which MRME stands for median of ratios of ME of a selected model to that of the unpenalized minimum square estimate under the full model. Both the columns of “C” and “IC” are measures of model complexity. Column “C” shows the average number of nonzero coefficients correctly estimated to be nonzero, and column “IC” presents the average number of zero coefficients incorrectly estimated to be nonzero. In the column labeled “underfit,” we present the proportion of excluding any nonzero coefficients in 200 replications. Likewise, we report the probability of selecting the exact subset model and the probability of including all three significant variables and some noise variables in the columns “correct-fit” and “overfit,” respectively.Table 1
Simulation results for linear regression models of Example8.
Method
MRME
Number of zeros
Proportion of
C
IC
Underfit
Correct-fit
Overfit
n
=
100, σ
=
2
Lasso
0.6213
3.0000
1.2050
0.0000
0.3700
0.6700
Alasso
0.3074
3.0000
0.3500
0.0000
0.7300
0.2700
SCAD
0.2715
3.0000
1.0650
0.0000
0.4400
0.5600
MCP
0.3041
3.0000
0.5650
0.0000
0.5800
0.4200
Dantzig
0.4623
3.0000
0.6546
0.0000
0.5700
0.4300
Bayesian
0.3548
3.0000
0.5732
0.0000
0.6300
0.3700
Atan
0.2550
3.0000
0.1750
0.0000
0.8450
0.1550
n
=
200, σ
=
2
Lasso
0.6027
3.0000
1.0700
0.0000
0.3550
0.6450
Alasso
0.2781
3.0000
0.1600
0.0000
0.8650
0.1350
SCAD
0.2900
3.0000
0.8550
0.0000
0.5250
0.4750
MCP
0.2752
3.0000
0.3650
0.0000
0.6850
0.3150
Dantzig
0.3863
3.0000
0.8576
0.0000
0.4920
0.5080
Bayesian
0.2563
3.0000
0.4754
0.0000
0.7150
0.2850
Atan
0.2508
3.0000
0.1000
0.0000
0.9050
0.0950As can be seen from Table1, all variable selection procedures dramatically reduce model error. Atan has the smallest model error among all competitors, followed by Alasso, MCP, SCAD, Bayesian, and Dantzig. In terms of sparsity, Atan also has the highest probability of correct-fit. The Atan penalty performs better than the other penalties. Also, Atan has some advantages when dimensional p is high which can be seen in Example 9.We now test the accuracy of our standard error formula (18). The median absolute deviation divided by 0.6745, denoted by SD in Table 2, of 3 estimated coefficients in the 200 simulations can be regarded as the true standard error. The median of the 200 estimated SDs, denoted by S
D
m, and the median absolute deviation error of the 200 estimated standard errors divided by 0.6745, denoted by S
D
m
a
d, gauge the overall performance of standard error formula (18). Table 2 presents the results for nonzero coefficients when sample size n
=
200. The results for the other case with n
=
100 are similar. Table 2 suggests that the sandwich formula performs surprisingly well.Table 2
Standard deviations of estimators for the linear regression model (n
=
200).
Method
⁢β
^
1
⁢β
^
2
⁢β
^
5
S
D
S
D
m (S
D
m
a
d)
S
D
S
D
m (S
D
m
a
d)
S
D
S
D
m (S
D
m
a
d)
Lasso
0.1753
0.1453 (0.0100)
0.1730
0.1688 (0.0094)
0.1591
0.1301 (0.0079)
Alasso
0.1483
0.1638 (0.0095)
0.1642
0.1636 (0.0098)
0.1475
0.1439 (0.0073)
SCAD
0.1797
0.1634 (0.0105)
0.1819
0.1608 (0.0104)
0.1398
0.1438 (0.0076)
MCP
0.1602
0.1643 (0.0096)
0.1861
0.1656 (0.0097)
0.1464
0.1435 (0.0069)
Dantzig
0.1734
0.1645 (0.0115)
0.1723
0.1665 (0.0094)
0.1581
0.1538 (0.0074)
Bayesian
0.1568
0.1635 (0.0089)
0.1678
0.1649 (0.0092)
0.1367
0.1375 (0.0078)
Atan
0.1510
0.1643 (0.0108)
0.1591
0.1658 (0.0096)
0.1609
0.1434 (0.0083)Example 9.
The example is from Wang et al. [23]. More specifically, we take β
=
(
11
/
4
,
-
23
/
6,37
/
12
,
-
13
/
9,1
/
3,0
,
0
,
…
,
0
)
T
∈
R
p and p
=
[
4
n
1
/
4
]
-
5 and [
t
] stands for the largest integer no larger than t. For this example, predictor dimension p is diverging but the dimension of the true model is fixed to be 5. Results from the simulation study are found in Table 3. A similar conclusion as in Example 8 can be found.Table 3
Simulation results for linear regression models of Example9.
Method
MRME
Number of zeros
Proportion of
C
IC
Underfit
Correct-fit
Overfit
n
=
200, p
=
10
Lasso
0.9955
4.9900
2.0100
0.0100
0.1150
0.8750
Alasso
0.5740
4.8650
0.1100
0.1350
0.8150
0.0500
SCAD
0.5659
4.9800
0.5600
0.0200
0.5600
0.4200
MCP
0.6177
4.9100
0.1200
0.0900
0.8250
0.0850
Dantzig
0.6987
4.8900
0.6700
0.1100
0.4360
0.4540
Bayesian
0.5656
4.8650
0.2500
0.1350
0.6340
0.2310
Atan
0.5447
4.8900
0.1150
0.1100
0.8250
0.0650
n
=
400, p
=
12
Lasso
1.2197
5.0000
2.0650
0.0000
0.1400
0.8600
Alasso
0.4458
4.9950
0.0900
0.0050
0.9250
0.0700
SCAD
0.4481
5.0000
0.4850
0.0000
0.6350
0.3650
MCP
0.4828
5.0000
0.1150
0.0000
0.8950
0.1050
Dantzig
0.7879
5.0000
0.5670
0.0000
0.3200
0.6800
Bayesian
0.4237
5.0000
0.1800
0.0000
0.7550
0.2450
Atan
0.4125
4.9950
0.0250
0.0050
0.9700
0.0250
n
=
800, p
=
16
Lasso
1.2004
5.0000
2.5700
0.0000
0.0900
0.9100
Alasso
0.3156
5.0000
0.0700
0.0000
0.9300
0.0700
SCAD
0.3219
5.0000
0.6550
0.0000
0.5950
0.4050
MCP
0.3220
5.0000
0.0750
0.0000
0.9300
0.0700
Dantzig
0.5791
5.0000
0.5470
0.0000
0.3400
0.6600
Bayesian
0.3275
5.0000
0.2800
0.0000
0.6750
0.3250
Atan
0.3239
5.0000
0.0750
0.0000
0.9350
0.0650Example 10.
In this simulation study presented here, we examined the performance of the various PLS methods forp substantially larger than in the previous studies. In particular, we took p
=
339, n
=
500, σ
5
=
5, and β
∗
=
2
I
37
T
,
-
3
I
37
T
,
I
37
T
,
0
228
T, where I
k
∈
R
k is the vector with all entries equal to 1. Thus, p
0
=
111. We simulated 200 independent datasets {
(
y
1
,
x
1
T
)
,
…
,
(
y
n
,
x
n
T
)
} in this study and, for each dataset, we computed estimates of β
∗. Results from this simulation study are found in Table 4.Table 4
Simulation results for linear regression models of Example10.
Method
MRME
Number of zeros
Proportion of
C
IC
Underfit
Correct-fit
Overfit
n
=
500, σ
=
5
Lasso
0.3492
110.2000
24.7450
0.5650
0.0000
0.4350
Alasso
0.2783
103.5250
6.1250
1.0000
0.0000
0.0000
SCAD
0.3012
106.0400
35.7150
0.9900
0.0000
0.0150
MCP
0.2791
103.3600
8.1100
1.0000
0.0000
0.0900
Dantzig
0.3120
108.3400
18.4650
0.7570
0.0000
0.2430
Bayesian
0.2876
104.4350
7.4700
0.9800
0.0000
0.0200
Atan
0.2794
101.4000
4.0600
1.0000
0.0000
0.0000Perhaps the most striking aspect of the results presented in Table4 is that hardly no method ever selected the correct model in this simulation study. However, given that p, p
0, and β
∗ are substantially larger in this study than in the previous simulation studies, this may not be too surprising. Notice that, on average, Atan selects the most parsimonious models of all methods and has the smallest model error. Atan’s nearest competitor in terms of model error is Alasso. This implementation of Alasso has mean model error 0.2783, but its average selected model size is 103.5250 larger than Atan’s. Since p
0
=
111, it is clear that Atan underfits in some instances. In fact, all of the methods in this study underfit to some extent. This may be due to the fact that many of the nonzero entries in β
∗ are small relative to the noise level σ
2
=
5.Example 11.
As an extension of this method, we consider the problem of simultaneous variable selection and estimation in the partially linear model:(21)
Y
=
X
′
β
+
g
T
+
ε
,where Y is a scalar response variate, X is a p-vector covariate, T is a scalar covariate and takes values in a compact interval (for simplicity, we assume this interval to be [
0,1
]), β is p
×
1 column vector of unknown regression parameter, function g
(
·
) is unknown, and model error ε is independent of (
X
,
T
) with mean 0. Traditionally, it has generally been assumed that β is finite dimension; several standard approaches, such as the kernel method, the spline method [41], and the local linear estimation [17], have been proposed.In this study, we simulaten
=
100,200 points T
i, i
=
1
,
…
,
100
(
200
), from the uniform distribution on 0,1. For each i, e
i
j
′
s are simulated to be normally distributed with autocorrelated variance structure A
R
(
ρ
), such that (22)
c
o
v
e
i
j
,
e
i
l
=
ρ
i
-
j
,
1
≤
i
,
j
≤
10
.
X
i
j’s are then formed as follows: (23)
X
i
1
=
sin
⁡
2
T
i
+
e
i
1
,
X
i
2
=
0.5
+
T
i
-
2
+
e
i
2
,
X
i
3
=
exp
⁡
T
i
+
e
i
3
,
X
i
5
=
T
i
-
0.7
4
+
e
i
5
,
X
i
6
=
T
i
1
+
T
i
2
-
1
+
e
i
6
,
X
i
5
=
1
+
T
i
+
e
i
7
,
X
i
8
=
log
⁡
3
T
i
+
8
+
e
i
8
.We investigate the scenario: p
=
10, X
i
j
=
e
i
j, j
=
4,9
,
10. In the scenario, we have β
j
=
j, 1
≤
j
≤
4, and β
j
=
0 with others, ε
~
N
(
0,1
). For each case, we repeated the simulation 200 times. We investigate g
(
·
) functions: g
(
t
)
=
cos
⁡
2
π
t. For comparison, we apply the different penalties in the parametric component with the B-spline in the nonparametric component.The results are summarized in Table5. Columns 2–5 in Table 5 are the averages of the estimates of β
j, j
=
1
,
…
,
4, respectively. Column 6 is the number of estimates of β
j, 5
≤
j
≤
p, which are 0, averaged over 200 simulations, and their medians are given in column 7. Model errors are computed as M
E
(
μ
^
)
=
(
β
^
-
β
)
T
E
(
x
x
T
)
(
β
^
-
β
). Their medians are listed in the 8th column, followed by the model errors standard deviations in parentheses.Table 5
Example11: comparison of estimators.
Estimator
β
1
β
2
β
3
β
4
C
-
C
~
MME (SD)
RASE (g
(
T
))
n
=
100, p
=
10
SCAD
0.9617
2.0128
2.8839
4.0611
4.7450
5.0000
0.0695 (0.0625)
0.6269
Lasso
0.8278
1.9517
2.7984
3.9274
5.5200
6.0000
0.1727 (0.0790)
0.8079
Dantzig
0.8767
1.9544
2.7894
3.9302
5.4590
6.0000
0.1696 (0.0680)
0.7278
Atan
0.9710
2.0121
2.9865
4.0472
4.9820
5.0000
0.0658 (0.0630)
0.6214
n
=
200, p
=
10
SCAD
0.9836
1.9989
2.9283
4.0278
5.2000
5.0000
0.0230 (0.0210)
0.3604
Lasso
0.9219
1.9529
2.9124
3.9580
5.1450
5.0000
0.0500 (0.0315)
0.4220
Dantzig
0.9534
1.9675
2.9096
3.9345
5.1460
5.0000
0.0510 (0.0290)
0.4154
Atan
0.9904
1.9946
2.9548
4.0189
5.1250
5.0000
0.0190 (0.0245)
0.3460The performance ofg
^
(
T
) is assessed by the square root of average squared errors (RASE):(24)
R
A
S
E
2
=
1
n
∑
i
=
1
n
g
^
T
i
-
g
T
i
2
,where {
T
i
,
i
=
1
,
…
,
n
} are the observed data points at which function g is estimated. Column 9 summarizes the RASE values for the different situations.We can see the following from Table5: (1) Comparing with the Lasso and Dantzig estimator, the SCAD and Atan estimators have a smaller model error, which is due to the fact that the SCAD and Atan estimators are unbiased estimators while the Lasso and Dantzig are biased especially for the larger coefficient. Moreover, the Atan and SCAD estimators are more stable. (2) Each method is able to select important variables, but it is obvious that the Atan estimator has slightly stronger sparsity. (3) For the nonparametric component, Atan and SCAD estimator have smaller RASE values.


### 5.2. Real Data Analysis

In this section, we apply the Atan regularization scheme to a prostate cancer example. The dataset in this example is derived from a study of prostate cancer in [42]. The dataset consists of the medical records of 97 patients who were about to receive a radical prostatectomy. The predictors are eight clinical measures: log
⁡(cancer volume) (lcavol), log (prostate weight) (lweight), age, the logarithm of the amount of benign prostatic hyperplasia (lbph), seminal vesicle invasion (svi), log
⁡(capsular penetration) (lcp), gleason score (gleason), and percentage gleason score 4 or 5 (pgg45). The response is the logarithm of prostate-specific antigen (lpsa). One of the main aims here is to identify which predictors are more important in predicting the response.The Lasso, Alasso, SCAD, MCP, and Atan are all applied to the data. We also compute the OLS estimate of the prostate cancer data. Results are summarized in Table6. The OLS estimator does not perform variable selection. Lasso selects five variables in the final model; SCAD selects lcavol, lweight, lbph, and svi in the final model, while Alasso, MCP, and Atan select lcavol, lweight, and svi. Thus, Atan selects a substantially simpler model than Lasso, SCAD. Furthermore, as indicated by the columns labeled R
2 (R
2 is equal to one minus the residual sum of squares divided by the total sum of squares) and R
2
/
R
O
L
S
2 in Table 6, the Atan estimator describes more variability in the data than Alasso and MCP and nearly as much as OLS estimator.Table 6
Prostate cancer data: comparing different methods.
Method
R
2
R
2
/
R
O
L
S
2
Variables selected
OLS
0.6615
1.0000
All
Lasso
0.5867
0.8870
(1, 2, 4, 5, 8)
Alasso
0.5991
0.9058
(1, 2, 5)
SCAD
0.6140
0.9283
(1, 2, 4, 5)
MCP
0.5999
0.9069
(1, 2, 5)
Atan
0.6057
0.9158
(1, 2, 5)


## 5.1. Simulation Studies

We now investigate the sparsity recovery and estimation properties of the proposed estimator via numerical simulations. We compare the following estimators: the Lasso estimator (implemented using R package glmnet [40]); the adaptive Lasso estimator (denoted by Alasso [15]), the SCAD estimator from the CD algorithm without calibration [31]; the MCP estimator from the CD algorithm with a
=
1.5 [31]; the Dantzig selector [7]; and Bayesian variable select method [8]. For the proposed Atan estimator, we take γ
=
0.005 and BIC statistic (15) is used to select the tuning parameter λ. In the following, we report simulation results from four examples.Example 8.
In this example, simulation data are generated from the linear regression model:(20)
y
=
x
T
β
∗
+
σ
ε
,where β
∗
=
(
3,1.5,0
,
0,2
,
0,0
,
0,0
,
0,0
,
0
)
T, ε
~
N
(
0,1
), and x is multivariate normal distribution with zero mean and covariance between the ith and jth elements being ρ
i
-
j with ρ
=
0.5. In our simulation, sample size n is set to be 100 and 200, σ
=
2. For each case, we repeated the simulation 200 times.For linear model, model error forμ
^
=
x
T
β
^ is M
E
(
μ
^
)
=
(
β
^
-
β
)
T
E
(
x
x
T
)
(
β
^
-
β
). Simulation results are summarized in Table 1, in which MRME stands for median of ratios of ME of a selected model to that of the unpenalized minimum square estimate under the full model. Both the columns of “C” and “IC” are measures of model complexity. Column “C” shows the average number of nonzero coefficients correctly estimated to be nonzero, and column “IC” presents the average number of zero coefficients incorrectly estimated to be nonzero. In the column labeled “underfit,” we present the proportion of excluding any nonzero coefficients in 200 replications. Likewise, we report the probability of selecting the exact subset model and the probability of including all three significant variables and some noise variables in the columns “correct-fit” and “overfit,” respectively.Table 1
Simulation results for linear regression models of Example8.
Method
MRME
Number of zeros
Proportion of
C
IC
Underfit
Correct-fit
Overfit
n
=
100, σ
=
2
Lasso
0.6213
3.0000
1.2050
0.0000
0.3700
0.6700
Alasso
0.3074
3.0000
0.3500
0.0000
0.7300
0.2700
SCAD
0.2715
3.0000
1.0650
0.0000
0.4400
0.5600
MCP
0.3041
3.0000
0.5650
0.0000
0.5800
0.4200
Dantzig
0.4623
3.0000
0.6546
0.0000
0.5700
0.4300
Bayesian
0.3548
3.0000
0.5732
0.0000
0.6300
0.3700
Atan
0.2550
3.0000
0.1750
0.0000
0.8450
0.1550
n
=
200, σ
=
2
Lasso
0.6027
3.0000
1.0700
0.0000
0.3550
0.6450
Alasso
0.2781
3.0000
0.1600
0.0000
0.8650
0.1350
SCAD
0.2900
3.0000
0.8550
0.0000
0.5250
0.4750
MCP
0.2752
3.0000
0.3650
0.0000
0.6850
0.3150
Dantzig
0.3863
3.0000
0.8576
0.0000
0.4920
0.5080
Bayesian
0.2563
3.0000
0.4754
0.0000
0.7150
0.2850
Atan
0.2508
3.0000
0.1000
0.0000
0.9050
0.0950As can be seen from Table1, all variable selection procedures dramatically reduce model error. Atan has the smallest model error among all competitors, followed by Alasso, MCP, SCAD, Bayesian, and Dantzig. In terms of sparsity, Atan also has the highest probability of correct-fit. The Atan penalty performs better than the other penalties. Also, Atan has some advantages when dimensional p is high which can be seen in Example 9.We now test the accuracy of our standard error formula (18). The median absolute deviation divided by 0.6745, denoted by SD in Table 2, of 3 estimated coefficients in the 200 simulations can be regarded as the true standard error. The median of the 200 estimated SDs, denoted by S
D
m, and the median absolute deviation error of the 200 estimated standard errors divided by 0.6745, denoted by S
D
m
a
d, gauge the overall performance of standard error formula (18). Table 2 presents the results for nonzero coefficients when sample size n
=
200. The results for the other case with n
=
100 are similar. Table 2 suggests that the sandwich formula performs surprisingly well.Table 2
Standard deviations of estimators for the linear regression model (n
=
200).
Method
⁢β
^
1
⁢β
^
2
⁢β
^
5
S
D
S
D
m (S
D
m
a
d)
S
D
S
D
m (S
D
m
a
d)
S
D
S
D
m (S
D
m
a
d)
Lasso
0.1753
0.1453 (0.0100)
0.1730
0.1688 (0.0094)
0.1591
0.1301 (0.0079)
Alasso
0.1483
0.1638 (0.0095)
0.1642
0.1636 (0.0098)
0.1475
0.1439 (0.0073)
SCAD
0.1797
0.1634 (0.0105)
0.1819
0.1608 (0.0104)
0.1398
0.1438 (0.0076)
MCP
0.1602
0.1643 (0.0096)
0.1861
0.1656 (0.0097)
0.1464
0.1435 (0.0069)
Dantzig
0.1734
0.1645 (0.0115)
0.1723
0.1665 (0.0094)
0.1581
0.1538 (0.0074)
Bayesian
0.1568
0.1635 (0.0089)
0.1678
0.1649 (0.0092)
0.1367
0.1375 (0.0078)
Atan
0.1510
0.1643 (0.0108)
0.1591
0.1658 (0.0096)
0.1609
0.1434 (0.0083)Example 9.
The example is from Wang et al. [23]. More specifically, we take β
=
(
11
/
4
,
-
23
/
6,37
/
12
,
-
13
/
9,1
/
3,0
,
0
,
…
,
0
)
T
∈
R
p and p
=
[
4
n
1
/
4
]
-
5 and [
t
] stands for the largest integer no larger than t. For this example, predictor dimension p is diverging but the dimension of the true model is fixed to be 5. Results from the simulation study are found in Table 3. A similar conclusion as in Example 8 can be found.Table 3
Simulation results for linear regression models of Example9.
Method
MRME
Number of zeros
Proportion of
C
IC
Underfit
Correct-fit
Overfit
n
=
200, p
=
10
Lasso
0.9955
4.9900
2.0100
0.0100
0.1150
0.8750
Alasso
0.5740
4.8650
0.1100
0.1350
0.8150
0.0500
SCAD
0.5659
4.9800
0.5600
0.0200
0.5600
0.4200
MCP
0.6177
4.9100
0.1200
0.0900
0.8250
0.0850
Dantzig
0.6987
4.8900
0.6700
0.1100
0.4360
0.4540
Bayesian
0.5656
4.8650
0.2500
0.1350
0.6340
0.2310
Atan
0.5447
4.8900
0.1150
0.1100
0.8250
0.0650
n
=
400, p
=
12
Lasso
1.2197
5.0000
2.0650
0.0000
0.1400
0.8600
Alasso
0.4458
4.9950
0.0900
0.0050
0.9250
0.0700
SCAD
0.4481
5.0000
0.4850
0.0000
0.6350
0.3650
MCP
0.4828
5.0000
0.1150
0.0000
0.8950
0.1050
Dantzig
0.7879
5.0000
0.5670
0.0000
0.3200
0.6800
Bayesian
0.4237
5.0000
0.1800
0.0000
0.7550
0.2450
Atan
0.4125
4.9950
0.0250
0.0050
0.9700
0.0250
n
=
800, p
=
16
Lasso
1.2004
5.0000
2.5700
0.0000
0.0900
0.9100
Alasso
0.3156
5.0000
0.0700
0.0000
0.9300
0.0700
SCAD
0.3219
5.0000
0.6550
0.0000
0.5950
0.4050
MCP
0.3220
5.0000
0.0750
0.0000
0.9300
0.0700
Dantzig
0.5791
5.0000
0.5470
0.0000
0.3400
0.6600
Bayesian
0.3275
5.0000
0.2800
0.0000
0.6750
0.3250
Atan
0.3239
5.0000
0.0750
0.0000
0.9350
0.0650Example 10.
In this simulation study presented here, we examined the performance of the various PLS methods forp substantially larger than in the previous studies. In particular, we took p
=
339, n
=
500, σ
5
=
5, and β
∗
=
2
I
37
T
,
-
3
I
37
T
,
I
37
T
,
0
228
T, where I
k
∈
R
k is the vector with all entries equal to 1. Thus, p
0
=
111. We simulated 200 independent datasets {
(
y
1
,
x
1
T
)
,
…
,
(
y
n
,
x
n
T
)
} in this study and, for each dataset, we computed estimates of β
∗. Results from this simulation study are found in Table 4.Table 4
Simulation results for linear regression models of Example10.
Method
MRME
Number of zeros
Proportion of
C
IC
Underfit
Correct-fit
Overfit
n
=
500, σ
=
5
Lasso
0.3492
110.2000
24.7450
0.5650
0.0000
0.4350
Alasso
0.2783
103.5250
6.1250
1.0000
0.0000
0.0000
SCAD
0.3012
106.0400
35.7150
0.9900
0.0000
0.0150
MCP
0.2791
103.3600
8.1100
1.0000
0.0000
0.0900
Dantzig
0.3120
108.3400
18.4650
0.7570
0.0000
0.2430
Bayesian
0.2876
104.4350
7.4700
0.9800
0.0000
0.0200
Atan
0.2794
101.4000
4.0600
1.0000
0.0000
0.0000Perhaps the most striking aspect of the results presented in Table4 is that hardly no method ever selected the correct model in this simulation study. However, given that p, p
0, and β
∗ are substantially larger in this study than in the previous simulation studies, this may not be too surprising. Notice that, on average, Atan selects the most parsimonious models of all methods and has the smallest model error. Atan’s nearest competitor in terms of model error is Alasso. This implementation of Alasso has mean model error 0.2783, but its average selected model size is 103.5250 larger than Atan’s. Since p
0
=
111, it is clear that Atan underfits in some instances. In fact, all of the methods in this study underfit to some extent. This may be due to the fact that many of the nonzero entries in β
∗ are small relative to the noise level σ
2
=
5.Example 11.
As an extension of this method, we consider the problem of simultaneous variable selection and estimation in the partially linear model:(21)
Y
=
X
′
β
+
g
T
+
ε
,where Y is a scalar response variate, X is a p-vector covariate, T is a scalar covariate and takes values in a compact interval (for simplicity, we assume this interval to be [
0,1
]), β is p
×
1 column vector of unknown regression parameter, function g
(
·
) is unknown, and model error ε is independent of (
X
,
T
) with mean 0. Traditionally, it has generally been assumed that β is finite dimension; several standard approaches, such as the kernel method, the spline method [41], and the local linear estimation [17], have been proposed.In this study, we simulaten
=
100,200 points T
i, i
=
1
,
…
,
100
(
200
), from the uniform distribution on 0,1. For each i, e
i
j
′
s are simulated to be normally distributed with autocorrelated variance structure A
R
(
ρ
), such that (22)
c
o
v
e
i
j
,
e
i
l
=
ρ
i
-
j
,
1
≤
i
,
j
≤
10
.
X
i
j’s are then formed as follows: (23)
X
i
1
=
sin
⁡
2
T
i
+
e
i
1
,
X
i
2
=
0.5
+
T
i
-
2
+
e
i
2
,
X
i
3
=
exp
⁡
T
i
+
e
i
3
,
X
i
5
=
T
i
-
0.7
4
+
e
i
5
,
X
i
6
=
T
i
1
+
T
i
2
-
1
+
e
i
6
,
X
i
5
=
1
+
T
i
+
e
i
7
,
X
i
8
=
log
⁡
3
T
i
+
8
+
e
i
8
.We investigate the scenario: p
=
10, X
i
j
=
e
i
j, j
=
4,9
,
10. In the scenario, we have β
j
=
j, 1
≤
j
≤
4, and β
j
=
0 with others, ε
~
N
(
0,1
). For each case, we repeated the simulation 200 times. We investigate g
(
·
) functions: g
(
t
)
=
cos
⁡
2
π
t. For comparison, we apply the different penalties in the parametric component with the B-spline in the nonparametric component.The results are summarized in Table5. Columns 2–5 in Table 5 are the averages of the estimates of β
j, j
=
1
,
…
,
4, respectively. Column 6 is the number of estimates of β
j, 5
≤
j
≤
p, which are 0, averaged over 200 simulations, and their medians are given in column 7. Model errors are computed as M
E
(
μ
^
)
=
(
β
^
-
β
)
T
E
(
x
x
T
)
(
β
^
-
β
). Their medians are listed in the 8th column, followed by the model errors standard deviations in parentheses.Table 5
Example11: comparison of estimators.
Estimator
β
1
β
2
β
3
β
4
C
-
C
~
MME (SD)
RASE (g
(
T
))
n
=
100, p
=
10
SCAD
0.9617
2.0128
2.8839
4.0611
4.7450
5.0000
0.0695 (0.0625)
0.6269
Lasso
0.8278
1.9517
2.7984
3.9274
5.5200
6.0000
0.1727 (0.0790)
0.8079
Dantzig
0.8767
1.9544
2.7894
3.9302
5.4590
6.0000
0.1696 (0.0680)
0.7278
Atan
0.9710
2.0121
2.9865
4.0472
4.9820
5.0000
0.0658 (0.0630)
0.6214
n
=
200, p
=
10
SCAD
0.9836
1.9989
2.9283
4.0278
5.2000
5.0000
0.0230 (0.0210)
0.3604
Lasso
0.9219
1.9529
2.9124
3.9580
5.1450
5.0000
0.0500 (0.0315)
0.4220
Dantzig
0.9534
1.9675
2.9096
3.9345
5.1460
5.0000
0.0510 (0.0290)
0.4154
Atan
0.9904
1.9946
2.9548
4.0189
5.1250
5.0000
0.0190 (0.0245)
0.3460The performance ofg
^
(
T
) is assessed by the square root of average squared errors (RASE):(24)
R
A
S
E
2
=
1
n
∑
i
=
1
n
g
^
T
i
-
g
T
i
2
,where {
T
i
,
i
=
1
,
…
,
n
} are the observed data points at which function g is estimated. Column 9 summarizes the RASE values for the different situations.We can see the following from Table5: (1) Comparing with the Lasso and Dantzig estimator, the SCAD and Atan estimators have a smaller model error, which is due to the fact that the SCAD and Atan estimators are unbiased estimators while the Lasso and Dantzig are biased especially for the larger coefficient. Moreover, the Atan and SCAD estimators are more stable. (2) Each method is able to select important variables, but it is obvious that the Atan estimator has slightly stronger sparsity. (3) For the nonparametric component, Atan and SCAD estimator have smaller RASE values.


## 5.2. Real Data Analysis

In this section, we apply the Atan regularization scheme to a prostate cancer example. The dataset in this example is derived from a study of prostate cancer in [42]. The dataset consists of the medical records of 97 patients who were about to receive a radical prostatectomy. The predictors are eight clinical measures: log
⁡(cancer volume) (lcavol), log (prostate weight) (lweight), age, the logarithm of the amount of benign prostatic hyperplasia (lbph), seminal vesicle invasion (svi), log
⁡(capsular penetration) (lcp), gleason score (gleason), and percentage gleason score 4 or 5 (pgg45). The response is the logarithm of prostate-specific antigen (lpsa). One of the main aims here is to identify which predictors are more important in predicting the response.The Lasso, Alasso, SCAD, MCP, and Atan are all applied to the data. We also compute the OLS estimate of the prostate cancer data. Results are summarized in Table6. The OLS estimator does not perform variable selection. Lasso selects five variables in the final model; SCAD selects lcavol, lweight, lbph, and svi in the final model, while Alasso, MCP, and Atan select lcavol, lweight, and svi. Thus, Atan selects a substantially simpler model than Lasso, SCAD. Furthermore, as indicated by the columns labeled R
2 (R
2 is equal to one minus the residual sum of squares divided by the total sum of squares) and R
2
/
R
O
L
S
2 in Table 6, the Atan estimator describes more variability in the data than Alasso and MCP and nearly as much as OLS estimator.Table 6
Prostate cancer data: comparing different methods.
Method
R
2
R
2
/
R
O
L
S
2
Variables selected
OLS
0.6615
1.0000
All
Lasso
0.5867
0.8870
(1, 2, 4, 5, 8)
Alasso
0.5991
0.9058
(1, 2, 5)
SCAD
0.6140
0.9283
(1, 2, 4, 5)
MCP
0.5999
0.9069
(1, 2, 5)
Atan
0.6057
0.9158
(1, 2, 5)


## 6. Conclusion and Discussion

In this paper, a new Atan penalty which very closely resemblesL
0 penalty is proposed. First, we establish the theory of the Atan estimator under mild conditions. The theory indicates that the Atan-penalized least squares procedure enjoys oracle properties even when the number of variables grows slower than the number of observations. Second, the iteratively reweighted Lasso algorithm makes our proposed estimator implementation fast. Third, we suggest a BIC-like tuning parameter selector to identify the true model consistently. Numerical studies further endorse our theoretical results and the advantage of the Atan estimator for model selection.We do not address the situation wherep
≫
n in this paper. In fact, the proposed Atan method can be easily extended for variable selection in the situation of p
≫
n. Also, as it is shown in Example 11, the Atan method can be applied to semiparametric model and nonparametric model [43, 44]. Furthermore, there is a recent field of applications of variable selection which is to look for impact points in functional data analysis [45, 46]. The possible combination of Atan estimator with the questions in [45, 46] would be interesting. These problems are beyond the scope of this paper and will be interesting topics for future research.



---

*Source: 6495417-2016-03-16.xml*