# Penalized Maximum Likelihood Method to a Class of Skewness Data Analysis

Authors: Xuedong Chen; Qianying Zeng; Qiankun Song
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

An extension of some standard likelihood and variable selection criteria based on procedures of linear regression models under the skew-normal distribution or the skew-t distribution is developed.
This novel class of models provides a useful generalization of symmetrical linear regression models, since the random term distributions cover both symmetric as well as asymmetric and heavy-tailed distributions. A generalized expectation-maximization algorithm is developed for computing the l
1 penalized estimator. Efficacy of the proposed methodology and algorithm is demonstrated by simulated data.


---

## Body

## 1. Introduction

In recent years, nonnormal distributions have received substantial interest in the statistics literature. The growing need for more flexible tools to analyze data that exhibit nonnormal features, including asymmetry, multimodality, and heavy tails, has led to intense development in nonnormal model-based methods. Many such distributions have been successfully applied to numerous datasets from a wide range of fields, including medical sciences, bioinformatics, environmetrics, engineering, economics, and financial sciences. Some recent applications of skew-normal and skew-t models include [1–5].The rich literature and active discussion of skew distributions as well as regression analysis of skew distributions were initiated by the pioneering work of [6], in which the univariate skew distribution was introduced. Following its generalization to the multivariate skew-normal distribution in [7], the number of contributions has grown rapidly. The concept of introducing additional parameters to regulate skewness in a distribution was subsequently extended to other parametric families. For a comprehensive survey of skew distributions, see, for example, the articles [4, 8–10].Besides the skew-normal distribution, which plays a central role in these developments, the skew-t distribution has also received much attention. Being a natural extension of the t distribution, the skew-t distribution retains reasonable tractability and is more robust against outliers than the skew-normal distribution.On the other hand, the Lasso is a popular method for regression analysis that uses anl
1 penalty to achieve a spares solution. This idea has been broadly applied, for example, to generalized linear models by [11] and Cox’s proportional hazard models for survival data by [12]. Inspired by an important work of [13], we will introduce a new approach to regression analysis based on a class skew distribution by taking advantage of sparse penalization methods. We propose to fit an l
1 penalized linear model with an l
1 penalty imposed on the location parameters of the skew distribution. Due to the sparse shrinkage property of the l
1 penalty, some components of regression coefficients are estimated by exact zero when the tuning parameter is properly chosen.The rest of the paper is organized as follows. In Section2, we introduce the skew-normal and skew-t distribution and present some of their properties. Section 3 outlines the EM-type algorithm and a penalty estimate of the proposed model of the skew distribution. In Section 4, a method of choice of tuning parameter is developed. The methodology proposed is illustrated in Section 5 by analyzing a simulated study.


## 2. The Skew-Normal Distribution and the Skew-t Distribution

### 2.1. The Skew-Normal Distribution

As developed by [6], a random variable Y follows a univariate skew-normal distribution with location parameter μ, scale parameter ω, and shape parameter α
∈
R, denoted by Y
~
SN
(
μ
,
ω
,
α
), if it has the density
(1)
f
SN
(
y
∣
μ
,
ω
,
α
)
=
2
ω
ϕ
(
y
-
μ
ω
)
Φ
(
α
y
-
μ
ω
)
,
where ϕ
(
·
) and Φ
(
·
) denote the standard normal density function and cumulative function, respectively. Note that if α
=
0, the density of Y reduces to the N
(
μ
,
ω
) density.A stochastic representation of theSN
(
μ
,
ω
,
α
) distribution can be given by
(2)
Y
=
μ
+
ω
δ
z
+
ω
1
-
δ
2
ɛ
,
where z
~
T
N
[
0
,
∞
)
(
0,1
) and T
N
[
a
,
b
]
(
μ
,
σ
2
) denotes the univariate truncated normal distribution with mean μ, variance σ
2 in interval [
a
,
b
]. ɛ
~
N
(
0,1
), independent of z, and δ
=
α
/
1
+
α
2; hence E
(
Y
)
=
μ
+
ω
δ
2
/
π.There is an alternative definition as discussed by [2, 9], which can be viewed as a new parameterization of above parameters; that is, ψ
=
ω
δ, σ
2
=
ω
2
(
1
-
δ
2
); furthermore, suppose u
~
N
(
0,1
) and η
~
N
(
0
,
σ
2
), independently; then (2) can be represented as
(3)
Y
=
μ
+
ψ
|
u
|
+
η
,
which defines a univariate skew-normal distribution SN
(
μ
,
ω
,
α
), where ω
2
=
ψ
2
+
σ
2 and α
=
ψ
/
σ, so (3) can be viewed as a reparametrization of (2). So it follows that the skew-normal distribution admits a convenient hierarchical characterization given by(i)
u
~
N
(
0,1
);
(ii)
Y
∣
u
~
N
(
μ
+
ψ
|
u
|
,
σ
2
).


### 2.2. The Skew-t Distribution

To achieve a higher degree of excess kurtosis and accommodate to data of a heavy tail, skew-t distributions have been introduced by [14]. A univariate random variable Y follows the skew-t distribution, Y
~
ST
(
μ
,
ω
,
α
,
ν
); the density function reads
(4)
f
ST
(
y
∣
μ
,
ω
,
α
,
ν
)
=
2
ω
t
ν
(
d
y
)
T
ν
+
1
(
α
d
y
ν
+
1
ν
+
d
y
2
)
,
where d
y
=
(
y
-
μ
)
/
ω and t
ν and T
ν denote, respectively, the density function and distribution function of a standard Student’s t distribution with ν degrees of freedom, and parameters μ, ω, and α are the same as above. Clearly, if Y
~
ST
(
μ
,
ω
,
α
,
ν
), it has the following stochastic representation:
(5)
Y
=
μ
+
ω
X
W
,
where X
~
SN
(
0,1
,
α
) and W
~
Γ
(
ν
/
2
,
ν
/
2
), independently and Γ
(
a
,
b
) denotes the Gamma distribution with parameters a and b.Similarly, the skew-t distribution has the convolution-type representation as
(6)
Y
=
μ
+
ψ
|
u
|
+
η
,
where u
~
t
1
(
0,1
,
ν
), η
~
t
1
(
0
,
σ
2
,
ν
), independently, and ψ and σ
2 are the same as above. Meanwhile, the skew-t distribution also has the hierarchical representation as(i)
w
~
Γ
(
ν
/
2
,
ν
/
2
);
(ii)
u
∣
w
~
N
(
0,1
/
w
);
(iii)
Y
∣
u
,
w
~
N
(
μ
+
ψ
|
u
|
,
σ
2
/
w
).


## 2.1. The Skew-Normal Distribution

As developed by [6], a random variable Y follows a univariate skew-normal distribution with location parameter μ, scale parameter ω, and shape parameter α
∈
R, denoted by Y
~
SN
(
μ
,
ω
,
α
), if it has the density
(1)
f
SN
(
y
∣
μ
,
ω
,
α
)
=
2
ω
ϕ
(
y
-
μ
ω
)
Φ
(
α
y
-
μ
ω
)
,
where ϕ
(
·
) and Φ
(
·
) denote the standard normal density function and cumulative function, respectively. Note that if α
=
0, the density of Y reduces to the N
(
μ
,
ω
) density.A stochastic representation of theSN
(
μ
,
ω
,
α
) distribution can be given by
(2)
Y
=
μ
+
ω
δ
z
+
ω
1
-
δ
2
ɛ
,
where z
~
T
N
[
0
,
∞
)
(
0,1
) and T
N
[
a
,
b
]
(
μ
,
σ
2
) denotes the univariate truncated normal distribution with mean μ, variance σ
2 in interval [
a
,
b
]. ɛ
~
N
(
0,1
), independent of z, and δ
=
α
/
1
+
α
2; hence E
(
Y
)
=
μ
+
ω
δ
2
/
π.There is an alternative definition as discussed by [2, 9], which can be viewed as a new parameterization of above parameters; that is, ψ
=
ω
δ, σ
2
=
ω
2
(
1
-
δ
2
); furthermore, suppose u
~
N
(
0,1
) and η
~
N
(
0
,
σ
2
), independently; then (2) can be represented as
(3)
Y
=
μ
+
ψ
|
u
|
+
η
,
which defines a univariate skew-normal distribution SN
(
μ
,
ω
,
α
), where ω
2
=
ψ
2
+
σ
2 and α
=
ψ
/
σ, so (3) can be viewed as a reparametrization of (2). So it follows that the skew-normal distribution admits a convenient hierarchical characterization given by(i)
u
~
N
(
0,1
);
(ii)
Y
∣
u
~
N
(
μ
+
ψ
|
u
|
,
σ
2
).


## 2.2. The Skew-t Distribution

To achieve a higher degree of excess kurtosis and accommodate to data of a heavy tail, skew-t distributions have been introduced by [14]. A univariate random variable Y follows the skew-t distribution, Y
~
ST
(
μ
,
ω
,
α
,
ν
); the density function reads
(4)
f
ST
(
y
∣
μ
,
ω
,
α
,
ν
)
=
2
ω
t
ν
(
d
y
)
T
ν
+
1
(
α
d
y
ν
+
1
ν
+
d
y
2
)
,
where d
y
=
(
y
-
μ
)
/
ω and t
ν and T
ν denote, respectively, the density function and distribution function of a standard Student’s t distribution with ν degrees of freedom, and parameters μ, ω, and α are the same as above. Clearly, if Y
~
ST
(
μ
,
ω
,
α
,
ν
), it has the following stochastic representation:
(5)
Y
=
μ
+
ω
X
W
,
where X
~
SN
(
0,1
,
α
) and W
~
Γ
(
ν
/
2
,
ν
/
2
), independently and Γ
(
a
,
b
) denotes the Gamma distribution with parameters a and b.Similarly, the skew-t distribution has the convolution-type representation as
(6)
Y
=
μ
+
ψ
|
u
|
+
η
,
where u
~
t
1
(
0,1
,
ν
), η
~
t
1
(
0
,
σ
2
,
ν
), independently, and ψ and σ
2 are the same as above. Meanwhile, the skew-t distribution also has the hierarchical representation as(i)
w
~
Γ
(
ν
/
2
,
ν
/
2
);
(ii)
u
∣
w
~
N
(
0,1
/
w
);
(iii)
Y
∣
u
,
w
~
N
(
μ
+
ψ
|
u
|
,
σ
2
/
w
).


## 3. Regression Model of Skew Distribution and EM Algorithm

In this section, we consider the linear regression model where the observation follows skew-normal distributions or skew-t distributions. In general, a linear regression model is defined as
(7)
y
i
=
x
i
β
+
ω
δ
2
π
+
ɛ
i
,
i
=
1
,
…
,
n
,
or
(8)
y
i
=
x
i
β
+
ω
δ
2
π
Γ
(
(
ν
-
1
)
/
2
)
Γ
(
ν
/
2
)
+
ɛ
i
,
i
=
1
,
…
,
n
,
where y
i are responses, parameter vector β
=
(
β
1
,
…
,
β
p
)
T, x
i is a vector of explanatory variable values, and the random error ɛ
i follows SN
(
0
,
ω
,
α
) or ST
(
0
,
ω
,
α
,
ν
), respectively, which corresponds to the regression model with an error distribution of which the mean is zero.It follows that the log-likelihood function for original parameterθ
=
(
β
,
ω
,
α
) or θ
=
(
β
,
ω
,
α
,
ν
), given the observed sample y
=
(
y
1
,
…
,
y
n
)
t from SN or ST, respectively, is given by
(9)
LL
(
θ
)
=
n
log
⁡
(
2
)
-
n
2
log
⁡
(
2
π
)
-
n
2
log
⁡
(
ω
2
)
+
∑
i
=
1
n
‍
log
⁡
(
K
i
)
,
where K
i is the function of ϕ
(
·
)
Φ
(
·
) or t
(
·
)
T
(
·
) in (1) or (4), according to the SN or ST distribution, respectively. Note that the Newton-Raphson method can be directly applied to get the ML estimate of the above log-likelihood, but in order to obtain a more efficient approach of estimation and variable selection of skew models, we discuss a more elaborate technique to find the ML estimate based on an EM algorithm.Firstly, for skew-normal distributions, letu
=
(
u
1
,
…
,
u
n
)
t; then under the hierarchical representation of Section 2.1, it follows that the complete log-likelihood function of modified parameter φ
=
(
β
,
σ
,
ψ
) associated with y
c
=
(
y
,
u
)
t is
(10)
l
c
(
φ
)
=
l
(
φ
∣
y
c
)
=
C
-
n
2
log
⁡
(
σ
2
)
-
1
2
σ
2
∑
i
=
1
n
‍
[
(
y
i
-
x
i
β
-
ψ
u
i
)
2
-
σ
2
u
i
2
]
,
where C is a constant that is independent of φ. Then consider the penalized log-likelihood defined by
(11)
PLL
(
φ
)
=
l
c
(
φ
)
-
1
2
∑
j
=
1
p
‍
P
λ
(
|
β
j
|
)
=
C
-
n
2
log
⁡
⁡
(
σ
2
)
-
1
2
σ
2
∑
i
=
1
n
‍
[
(
y
i
-
x
i
β
-
ψ
u
i
)
2
-
σ
2
u
i
2
]
-
1
2
∑
j
=
1
p
‍
P
λ
(
|
β
j
|
)
,
where P
λ
(
·
) is a nonnegative penalty function and we use P
λ
(
|
β
j
|
)
=
λ
|
β
j
|, which is the Lasso penalty, where λ is the tuning parameter not included in φ.The Lasso estimator is then defined asargmaxPLL
(
φ
); that is,
(12)
φ
^
=
argmin
{
λ
2
∑
j
=
1
p
‍
|
β
j
|
∑
i
=
1
n
‍
n
2
log
⁡
⁡
(
σ
2
)
-
1
2
σ
2
∑
i
=
1
n
‍
[
(
y
i
-
x
i
β
-
ψ
u
i
)
2
-
σ
2
u
i
2
]
+
λ
2
∑
j
=
1
p
‍
|
β
j
|
}
.Clearly,u
=
(
u
1
,
…
,
u
n
) can be seen as “missing data”; in the following, we will derive a generalized expectation-maximization algorithm for computing the l
1 penalized estimator. It turns out that the penalized estimator can be computed by iterative Lasso-penalized least square.Set(13)
Q
(
φ
∣
φ
^
(
k
)
)
=
E
[
l
c
(
φ
)
∣
y
,
φ
^
(
k
)
]
=
E
u
∣
y
[
l
(
φ
∣
y
c
)
∣
y
,
φ
^
(
k
)
]
as the Q-function over φ at iteration k step and where E
u
∣
y denotes the conditional expectation of u given y.At theE-step, by using the properties of conditional expectation of the truncated normal distribution, we obtain
(14)
e
^
1
,
i
(
k
)
=
E
[
u
i
∣
y
i
,
φ
^
(
k
)
]
=
μ
^
1
,
i
(
k
)
+
σ
^
1
(
k
)
ϕ
(
μ
^
1
,
i
(
k
)
/
σ
^
1
(
k
)
)
Φ
(
μ
^
1
,
i
(
k
)
/
σ
^
1
(
k
)
)
,
e
^
2
,
i
(
k
)
=
E
[
u
i
2
∣
y
i
,
φ
^
(
k
)
]
=
(
μ
^
1
,
i
(
k
)
)
2
+
(
σ
^
1
(
k
)
)
2
+
μ
^
1
,
i
(
k
)
σ
^
1
(
k
)
ϕ
(
μ
^
1
,
i
(
k
)
/
σ
^
1
(
k
)
)
Φ
(
μ
^
1
,
i
(
k
)
/
σ
^
1
(
k
)
)
,
where μ
^
1
,
i
(
k
)
=
ψ
^
(
k
)
(
y
i
-
x
i
β
^
(
k
)
)
/
[
(
ψ
^
(
k
)
)
2
+
(
σ
^
(
k
)
)
2
] and σ
^
1
(
k
)
=
σ
^
(
k
)
/
(
ψ
^
(
k
)
)
2
+
(
σ
^
(
k
)
)
2.Then, combining theQ-function and the Lasso penalty, we set
(15)
Q
~
(
φ
∣
φ
^
(
k
)
)
=
Q
(
φ
∣
φ
^
(
k
)
)
-
λ
2
∑
j
=
1
p
‍
|
β
j
|
,
as the modified Q-function.At theM-step update φ by maximizing Q
~
(
φ
∣
φ
^
(
k
)
) over φ in a sequence of conditional maximization steps. Firstly, we get the following expression:
(16)
β
^
(
k
+
1
)
=
argmin
{
λ
2
∑
j
=
1
p
‍
|
β
j
|
1
2
(
σ
^
(
k
)
)
2
∑
i
=
1
n
‍
[
y
i
-
x
i
β
-
ψ
^
(
k
)
e
^
1
,
i
(
k
)
]
+
λ
2
∑
j
=
1
p
‍
|
β
j
|
1
2
(
σ
^
(
k
)
)
2
∑
i
=
1
n
‍
}
;
then, we can obtain the closed expression of other parameters,
(17)
ψ
^
(
k
+
1
)
=
∑
i
=
1
n
(
y
i
-
x
i
β
^
(
k
+
1
)
)
e
^
1
,
i
(
k
)
∑
i
=
1
n
e
^
2
,
i
(
k
)
,
σ
2
^
(
k
+
1
)
=
1
n
∑
i
=
1
n
‍
[
(
y
i
-
x
i
β
^
(
k
+
1
)
)
2
+
(
ψ
^
(
k
+
1
)
)
2
e
^
2
,
i
(
k
)
-
2
(
y
i
-
x
i
β
^
(
k
+
1
)
)
ψ
^
(
k
+
1
)
e
^
2
,
i
(
k
)
(
y
i
-
x
i
β
^
(
k
+
1
)
)
2
]
.
So, the above procedures can be summarized as follows.Step 1:
set initial values forβ, ψ, σ, and tuning parameter λ.
Step 2:
compute the value ofe
^
1
,
i and e
^
2
,
i at current iteration step as the “working observed value.”
Step 3:
solve the penalized least squares problem about parameterβ and obtain its new Lasso estimator.
Step 4:
compute the new estimator ofψ and σ.
Step 5:
repeat steps 2–4 till convergence.Secondly, for skew-t distributions, let Y
=
(
y
1
,
…
,
y
n
), u
=
(
u
1
,
…
,
u
n
), w
=
(
w
1
,
…
,
w
n
); then, based on the hierarchical representation of Section 2.2, we can get the complete log-likelihood function associated with y
c
=
(
y
,
u
,
w
) as
(18)
l
S
T
c
(
φ
)
=
l
ST
(
φ
∣
y
c
)
=
C
-
n
2
log
⁡
(
σ
2
)
-
1
2
σ
2
∑
i
=
1
n
‍
w
i
(
y
i
-
x
i
β
-
ψ
u
i
)
2
-
1
2
∑
i
=
1
n
‍
w
i
u
i
2
+
n
ν
2
log
⁡
⁡
(
ν
2
)
-
n
log
⁡
Γ
(
ν
2
)
+
ν
2
∑
i
=
1
n
‍
log
⁡
(
w
i
)
-
ν
2
∑
i
=
1
n
‍
w
i
.Similarly, setting(19)
PL
L
ST
(
φ
)
=
l
S
T
c
(
φ
)
-
1
2
∑
j
=
1
p
‍
P
λ
(
|
β
j
|
)
,
Q
(
φ
∣
φ
^
(
k
)
)
=
E
[
l
S
T
c
(
φ
)
∣
y
,
φ
^
(
k
)
]
,
then we can combine the E-step, M-step, and the Lasso estimator together as
(20)
β
^
(
k
+
1
)
=
argmin
{
1
2
(
σ
^
(
k
)
)
2
∑
i
=
1
n
‍
s
^
1
,
i
(
k
)
[
y
i
-
x
i
β
-
ψ
^
(
k
)
s
^
e
1
(
k
)
]
+
λ
2
∑
j
=
1
p
‍
|
β
j
|
}
,
ψ
^
(
k
+
1
)
=
∑
i
=
1
n
(
y
i
-
x
i
β
^
(
k
+
1
)
)
s
^
e
1
,
i
(
k
)
∑
i
=
1
n
s
^
e
2
,
i
(
k
)
,
σ
2
^
(
k
+
1
)
=
1
n
∑
i
=
1
n
‍
[
s
^
1
,
i
(
k
)
(
y
i
-
x
i
β
^
(
k
+
1
)
)
2
+
(
ψ
^
(
k
+
1
)
)
2
s
^
e
2
,
i
(
k
)
-
2
(
y
i
-
x
i
β
^
(
k
+
1
)
)
ψ
^
(
k
+
1
)
e
^
2
,
i
(
k
)
s
^
1
,
i
(
k
)
(
y
i
-
x
i
β
^
(
k
+
1
)
)
2
]
,
and parameter ν
(
k
+
1
) can be obtained via estimating equation as follows:
(21)
n
[
log
⁡
(
ν
2
)
-
digamma
(
ν
2
)
+
1
]
+
∑
i
=
1
n
‍
(
s
^
2
,
i
(
k
)
-
s
^
1
,
i
(
k
)
)
=
0
,
where
(22)
s
^
1
,
i
(
k
)
=
E
φ
[
w
i
∣
y
i
,
φ
^
(
k
)
]
=
(
ν
+
1
ν
+
d
y
i
)
T
ν
+
3
(
y
2
,
i
/
α
)
T
ν
+
1
(
y
2
,
i
/
α
)
,
s
^
2
,
i
(
k
)
=
E
φ
[
log
⁡
(
w
i
)
∣
y
i
,
φ
^
(
k
)
]
=
s
^
1
,
i
(
k
)
-
log
⁡
(
ν
+
d
y
i
2
)
-
digamma
(
ν
+
1
2
)
-
1
,
s
^
e
1
,
i
(
k
)
=
E
φ
[
w
i
u
i
∣
y
i
,
φ
^
(
k
)
]
=
α
ν
+
1
ν
+
d
y
i
t
ν
+
1
(
y
2
,
i
/
α
)
T
ν
+
1
(
y
2
,
i
/
α
)
+
s
^
2
,
i
q
i
,
s
^
e
2
,
i
(
k
)
=
E
φ
[
w
i
u
i
2
∣
y
i
,
φ
^
(
k
)
]
=
(
α
)
2
+
α
q
i
t
ν
+
1
(
y
2
,
i
/
α
)
T
ν
+
1
(
y
2
,
i
/
α
)
+
s
^
2
,
i
(
q
i
)
2
,
and d
y
i
=
(
y
i
-
x
i
β
)
2
/
(
σ
2
+
ψ
2
), q
i
=
ψ
(
y
i
-
x
i
β
)
/
(
σ
2
+
ψ
2
), α
=
ψ
/
σ
2
+
ψ
2, y
1
,
i
=
q
i
(
ν
+
1
)
/
(
ν
+
d
y
i
), and y
2
,
i
=
q
i
(
ν
+
3
)
/
(
ν
+
d
y
i
).Similarly, we can also summarize the above procedures as follows.Step 1:
set initial values forβ, ψ, σ, ν, and tuning parameter λ.
Step 2:
compute the value ofs
^
1
,
i, s
^
2
,
i and s
^
e
1
,
i, s
^
e
2
,
i at current iteration step as the “working observed value.”
Step 3:
solve the penalized least squares problem about parameterβ and obtain its new Lasso estimator.
Step 4:
compute the new estimates ofψ and σ directive and then get the new estimates of ν.
Step 5:
get the tuning parametersλ by GCV.
Step 6:
repeat steps 2–5 till convergence.


## 4. Choice of Tuning Parameters

We now consider the choice of tuning parameters. In application, the cross-validation (CV), or generalized cross-validation (GCV), is often used for choosing tuning parameters. Following the examples of [11, 15], we develop a componentwise deviance-based GCV criterion for the proposed models.Letφ
~ be the EM estimates of parameters φ without penalty functions. For a given value of tuning parameter of λ, let β
^ be the maximum penalized likelihood estimates of the parameters of the model by fixing the rest of components of φ at φ
~. Denote the deviance function as
(23)
D
(
β
^
)
=
E
[
l
c
(
φ
)
∣
β
^
]
-
E
[
l
c
(
φ
)
∣
φ
~
]
.
Further, let l
'
'
(
β
^
) be the second derivative of the log-likelihood function with respect to β evaluated at β
^. We define a GCV criterion of the model as
(24)
GCV
(
λ
)
=
D
(
β
^
)
n
(
1
-
e
(
λ
)
/
n
)
2
,
where e
(
λ
) is the effective number of regression coefficients. It is given by
(25)
e
(
λ
)
=
trace
{
[
l
′′
(
β
^
)
-
Σ
(
β
^
)
]
-
1
l
′′
(
β
^
)
}
,
where Σ
(
β
^
)
=
diag
⁡
{
1
/
|
β
^
1
|
,
…
,
1
/
|
β
^
p
|
}. The tuning parameter λ is chosen by minimizing GCV
(
λ
).


## 5. Simulation Study

In order to investigate the experimental performance of our methodology, we undertake a simulation study to investigate the effect of estimation and variable selection. The dataset is generated as follows: letx
i
=
(
x
1
i
,
…
,
x
p
i
)
T be covariate, where p
=
30 and x
1
i
,
…
,
p
1 are independent and distributed as uniform distribution U
(
-
1,1
). Moreover, assume response y
i is formulated as
(26)
y
i
=
x
i
T
β
+
ɛ
i
,
where ɛ
i is distributed as the following five cases.Case(
1
): ɛ
i
~
N
(
0
,
σ
2
); Case (
2
): ɛ
i
~
SN
(
0
,
σ
1
2
,
α
1
); Case (
3
): ɛ
i
~
SN
(
0
,
σ
1
2
,
α
2
); Case (
4
): ɛ
i
~
ST
(
0
,
σ
2
2
,
α
1
,
ν
); and Case (
5
): ɛ
i
~
ST
(
0
,
σ
2
2
,
α
2
,
ν
).We take the sample sizen to be 100. The true values for β are β
1
=
β
2
=
β
3
=
1 and the other components of β are equal to zero. Further, let σ
1
=
2.5, σ
2
=
1.0, α
1
=
1.5, α
2
=
-
0.8, and ν
=
3.Under above settings, we use three linear regression models based on normal, skew-normal, and skew-t distribution to fit each dataset, respectively.Table1 shows the results of the simulation study, where No. denotes the number of nonzero components of β, which can determine the effect of variable selection, and where AE denotes the average of absolute error; that is, (
1
/
p
)
∑
j
=
1
p
‍
|
β
j
0
-
β
^
j
|, which can show the accuracy of the estimates, and the λ is the tuning parameter.Table 1
Results of the simulation study.
Case
N
SN
ST
No.
AE
λ
No.
AE
λ
No.
AE
λ
1
3
0.6238
9.45
3
0.7323
7.46
3
0.7431
6.91
2
7
2.3536
12.20
3
1.0913
6.38
3
1.0396
6.22
3
11
3.3893
14.80
3
1.0488
6.37
3
1.1158
6.12
4
8
2.8939
18.92
3
1.1261
6.12
3
1.0135
6.25
5
10
3.8273
16.82
3
1.1542
5.12
3
1.0363
5.05Examination of Table1 reveals that (
i
) when the true distribution is normal, the difference of the three models is very small and the coefficient β can be found correctly; (
ii
) when the true distribution is nonnormal, the models based on normal fit the data badly; not only does the proportion of nonzero increase but also the average error becomes significant; (
iii
) for the cases of the true distribution being nonnormal, the models based on the skew distribution can improve the fitting effect better than the ones of the normal distribution.Figures1(a)–1(c) present the trace plots of coefficients versus tuning parameters: (a) Case (
5
) fitted by normal distribution (N) model; (b) Case (
5
) fitted by skew-normal distributions (SN) model; (c) Case (
5
) fitted by skew-t distributions (ST) model. Moreover, (d) shows profiles of the tuning parameter of (c).(a)–(c) Trace plot of regression coefficients and profile of tuning parameter.
(a)
(b)
(c)
(d)From Figure1(c), it is clear that the ST model is most appropriate for fitting this dataset.


## 6. Conclusion

In this paper, we have proposed a class of linear regression models based on the asymmetric distribution. An EM-type algorithm and the Lasso penalty estimates are developed by exploring the statistical properties of the model that can be implemented efficiently in types of software such as SAS, R, and Matlab. Meanwhile, we have shown that thel
1 penalized maximum likelihood estimation can be done via a generalized EM algorithm that is equivalent to iterative l
1 penalized least square for some parameters. We have observed that if the data are generated from a true model of the skew distribution, the l
1 penalized models based on EM algorithm not only identify zero coefficient but also are significantly more accurate than the general MLE models.On the other hand, due to recent advances in computational technology, it is worthwhile to carry out Bayesian treatments via Markov chain Monte Carlo (MCMC) sampling methods in the context of this paper.



---

*Source: 824816-2014-09-29.xml*