Variable Selection for Robust Mixture Regression Model with Skew Scale Mixtures of Normal Distributions

Tingzhu Chen, Wanzhou Ye

Department of Mathematics, College of Science, Shanghai University, Shanghai, China.

DOI: 10.4236/apm.2022.123010                

Abstract

In this paper, we propose a robust mixture regression model based on the skew scale mixtures of normal distributions (RMR-SSMN) which can accommodate asymmetric, heavy-tailed and contaminated data better. For the variable selection problem, the penalized likelihood approach with a new combined penalty function which balances the SCAD and l₂ penalty is proposed. The adjusted EM algorithm is presented to get parameter estimates of RMR-SSMN models at a faster convergence rate. As simulations show, our mixture models are more robust than general FMR models and the new combined penalty function outperforms SCAD for variable selection. Finally, the proposed methodology and algorithm are applied to a real data set and achieve reasonable results.

Keywords

Robust Mixture Regression Model, Skew Scale Mixtures of Normal Distributions, EM Algorithm, SCAD Penalty

1. Introduction

In applied statistics, the arc-sine laws for the Wiener process and the skew Brownian motion [1] are widely used in finance if the market is homogeneous. However, the problem of modeling heterogeneous data has been extensively studied in recent years and the Finite Mixture of Regression (FMR) model is an important tool for heterogeneous cases. A large number of applications associate a random response variable $Y$ with covariates $x$ through FMR models and the assumption is that for each observation data point $\left(x_1,Y_1\right),\cdots ,\left(x_n,Y_n\right)$, the regression coefficients are not the same. More details about the FMR model can be found in [2].

The Gaussian FMR model is the most common FMR model, which assumes that the random error of each subgroup follows the normal distribution. It is well known that using the normal distribution to model data with asymmetric and heavy-tailed behaviors is unsuitable, and the parameter estimates are sensitive to outliers. To overcome the potential shortcomings of Gaussian mixture models, McLachlan et al. [3] proposed to replace the mixtures of normal with mixtures of t-distribution which results in a more robust mixture model. Basso et al. [4] studied a finite mixture model based on scale mixtures of skew-normal distributions, and Franczak et al. [5] proposed a mixture model using shifted asymmetric Laplace distributions which parameterize the skewness as well as the location and the scale.

The problem of variable selection in FMR models has been widely discussed recently. There are generally two types of variable selection methods. One is the optimal subset selection method and the discontinuous penalty method based on the information criterion, including stepwise regression, best subset regression, BIC criterion, AIC criterion and so on. The other is the continuous penalty method. By imposing penalties on the parameters of the objective function, one can select significant variables and obtain the parameter estimates simultaneously. The Least Absolute Shrinkage and Selection Operator (LASSO), elastic net regularization [6], MCP penalty [7] and SCAD penalty [8] are penalty functions for variable selection. We utilize the SCAD and a new penalty function proposed in this paper which balance the SCAD and $l_2$ penalty to perform variable selection on a robust mixture regression model based on Skew Scale Mixtures of Normal (SSMN) distributions [9] and this robust model can accommodate asymmetric and heavy-tailed data better.

The paper is organized as follows. In Section 2, a robust mixture regression model using the skew scale mixtures of normal distributions (RMR-SSMN) is introduced. Then, variable selection methods with SCAD penalty function and a newly proposed penalty function are presented in Section 3. Section 4 outlines the adjusted EM algorithm for estimating and a BIC method for selecting turning parameters and components. In Section 5, we carry out simulation studies to compare the performances between FMR models and RMR-SSMN models, and show the effect of variable selection with penalty functions. An application to a real data set of the method is discussed in Section 6 and some conclusions are obtained in Section 7.

2. Robust Mixture Regression Model with SSMN Distributions

It is known that the FMR model can model heterogeneous data and the Skew Scale Mixtures of Normal (SSMN) distributions [9] cover both asymmetric and heavy-tailed distributions. Therefore, we propose a robust mixture regression model whose regression errors of components follow SSMN distributions. Unsurprisingly, this model is more robust than general FMR models for heterogeneous cases.

2.1. Skew Scale Mixtures of Normal Distributions

If a random variable $Y$ follows a skew-normal (SN) distribution with location parameter $\mu \in ℝ$, scale parameter ${\sigma }^2$, and skewness parameter $\lambda$, denoted $Y\sim SN\left(\mu ,{\sigma }^2,\lambda \right)$, then its density function is given as follows:

$f\left(y\right)=2\varphi \left(y;\mu ,{\sigma }^2\right)\Phi \left(\lambda \frac{y-\mu }{\sigma }\right)$. (1)

where $\varphi \left(y;\mu ,{\sigma }^2\right)$ and $\Phi \left(y;\mu ,{\sigma }^2\right)$ are the probability density function and the cumulative distribution function of $N\left(\mu ,{\sigma }^2\right)$ calculated at $y$, respectively.

Note that when $\lambda =0$, the $SN\left(\mu ,{\sigma }^2,\lambda \right)$ reduces to the $N\left(\mu ,{\sigma }^2\right)$, and as given in [9], the SN distribution’s marginal stochastic representation is presented by:

$Y=\mu +\sigma \left(\delta \left|T_0\right|+{\left(1-{\delta }^2\right)}^{\frac{1}{2}}{T}_1\right),$ (2)

where $\delta =\lambda /{\left(1+{\lambda }^2\right)}^{1/2}$, $T_0\sim N\left(0,1\right)$ and $T_1\sim N\left(0,1\right)$ are independent.

Furthermore, if a random variable $Y$ follows a SSMN distribution [9] with location parameter $\mu \in ℝ$, scale parameter ${\sigma }^2$, and skewness parameter $\lambda$, then its probability density function is given by:

$f\left(y\right)=2{\displaystyle {\int }_0^{\infty }\varphi \left(y;\mu ,\mathcal{l}\left(u\right){\sigma }^2\right)\text{d}H\left(u;\tau \right)\Phi }\left(\lambda \frac{y-\mu }{\sigma }\right),$ (3)

where $H\left(u;\tau \right)$ is the cumulative distribution function of $U$ who derived from the parameter vector $\tau$, and $U$ is a positive random variable, and $\mathcal{l}\left(u\right)$ is a strictly positive function. If the probability density function of the random variable $Y$ is shown as the Equation (3), it can be denoted as: $Y\sim SSMN\left(\mu ,{\sigma }^2,\lambda ,H;\mathcal{l}\right)$.

For $Y\sim SSMN\left(\mu ,{\sigma }^2,\lambda ,H;\mathcal{l}\right)$, its hierarchical representation has the form as follows:

$Y|U=u\sim SN\left(\mu ,\mathcal{l}\left(u\right){\sigma }^2,\lambda {\mathcal{l}}^{\frac{1}{2}}\left(u\right)\right),U\sim H\left(\tau \right).$ (4)

This paper will consider the following distributions in the SSMN distributions family:

· The skew Student-t-normal distribution (STN) [10] with $U\sim Gamma\left(v/2,v/2\right)$, $v>0$ and $\mathcal{l}\left(u\right)=1/u$, which follows probability density:

$f\left(y\right)=2\frac{1}{\sigma \sqrt{v\pi }}\frac{\Gamma \left(\left(v+1\right)/2\right)}{\Gamma \left(v/2\right)}{\left(1+\frac{d}{v}\right)}^{-\left(\frac{v+1}{2}\right)}\Phi \left(\lambda \frac{y-\mu }{\sigma }\right),$ (5)

where $d={\left(y-\mu \right)}^2/{\sigma }^2$ and $\Gamma (.)$ is the gamma function. We can obtain that $U|Y=y\sim Gamma\left(\left(v+1\right)/2,\left(v+d\right)/2\right)$.

· The skew contaminated normal distribution (SCN). $U$ is a discrete random variable taking one of two values and $\mathcal{l}\left(u\right)=1/u$. Given the parameter vector $\tau ={\left(v,\gamma \right)}^T$, $0<v<1$, $0<\gamma <1$, the density function of $U$ is $h\left(u;\tau \right)=v{\mathbb{I}}_{\left(u=\gamma \right)}+\left(1-v\right){\mathbb{I}}_{\left(u=1\right)}$. Naturally get as follows:

$f\left(y\right)=2\left\{v\varphi \left(y;\mu ,{\sigma }^2/\gamma \right)+\left(1-v\right)\varphi \left(y;\mu ,{\sigma }^2\right)\right\}\Phi \left(\lambda \frac{y-\mu }{\sigma }\right).$ (6)

Therefore, the conditional distribution $U|Y=y$ can be obtained as:

$f\left(u|Y=y\right)=\frac{1}{f_0\left(y\right)}\left\{v\varphi \left(y;\mu ,{\sigma }^2/\gamma \right){\mathbb{I}}_{\left(u=\gamma \right)}+\left(1-v\right)\varphi \left(y;\mu ,{\sigma }^2\right){\mathbb{I}}_{\left(u=1\right)}\right\},$ (7)

where $f_0\left(y\right)=v\varphi \left(y;\mu ,{\sigma }^2/\gamma \right){\mathbb{I}}_{\left(u=\gamma \right)}+\left(1-v\right)\varphi \left(y;\mu ,{\sigma }^2\right){\mathbb{I}}_{\left(u=1\right)}$.

· The skew power-exponential distribution (SPE) has following probability density:

$f\left(y\right)=2\frac{v}{2^{1/2v}\sigma \Gamma \left(1/2v\right)}{\text{e}}^{-d^v/2}\Phi \left(\lambda \frac{y-\mu }{\sigma }\right),0.5<v\le 1,$ (8)

Ferreira et al. [9] have proved that $E\left[{\mathcal{l}}^{-1}\left(U\right)|Y=y\right]=v{d}^{v-1}$.

2.2. Robust Mixture Regression Model with SSMN Distributions

Suppose we have $n$ independent random variables $y_1,y_2,\cdots ,y_n$, which are taken from a mixture of SSMN distributions. The conditional density function of the robust mixture regression model with SSMN distributions (RMR-SSMN) which has $K$ components is given by:

$f\left(y_i;x_i,\Psi \right)={\displaystyle \underset{k=1}{\overset{K}{\sum }}{\omega }_{k}SSMN}\left(y_i;{\alpha }_k+x_i^T{\beta }_k,{\sigma }_k^2,{\lambda }_k,{\tau }_k;\mathcal{l}\right),$ (9)

with covariate vector $x_i\in {ℝ}^q$ and q-dimensional unknown regression coefficients vector ${\beta }_k,k=1,\cdots ,K$. ${\omega }_k,k=1,\cdots ,K$ denote the mixing proportions satisfying ${\omega }_k\ge 0$, ${\displaystyle {\sum }_{k=1}^K{\omega }_k}=1$.

$\Psi ={\left({\omega }_1,\cdots ,{\omega }_{K-1},{\alpha }_1,\cdots ,{\alpha }_K,{\beta }_1^T,\cdots ,{\beta }_K^T,{\sigma }_1^2,\cdots ,{\sigma }_K^2,{\lambda }_1,\cdots ,{\lambda }_K,{\tau }_1^T,\cdots ,{\tau }_K^T\right)}^T$ is the parameter vector of the model. For convenience, let $\omega ={\left({\omega }_1,\cdots ,{\omega }_{K-1}\right)}^T$, $\alpha ={\left({\alpha }_1,\cdots ,{\alpha }_K\right)}^T$, $\beta ={\left({\beta }_1^T,\cdots ,{\beta }_K^T\right)}^T$, ${\sigma }^2={\left({\sigma }_1^2,\cdots ,{\sigma }_K^2\right)}^T$, $\lambda ={\left({\lambda }_1,\cdots ,{\lambda }_K\right)}^T$, and ${\tau }^{\ast }={\left({\tau }_1^T,\cdots ,{\tau }_K^T\right)}^T$. In this paper, RMR-SSMN models contain the robust mixture regression model with STN distribution (RMR-STN), SCN distribution (RMR-SCN), SPE distribution (RMR-SPE) and SN distribution (RMR-SN).

3. Variable Selection Method

If a component in the q-dimensional explanatory variable $x$ has no significant effect on the response variable $y$, the regression coefficient of this component estimated by the maximum likelihood method will close to 0 rather than 0. Thus, this covariate is not excluded from the model and makes the model unstable. To avoid this problem, we use a penalized likelihood approach [11] for selecting variables and estimating parameters simultaneously. Let $\left\{\left(y_i,x_i\right);i=1,\cdots ,n\right\}$ be sample observations from RMR-SSMN models. The log-likelihood function of $\Psi$ is given by:

$l\left(\Psi \right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\log {\displaystyle \underset{k=1}{\overset{K}{\sum }}{\omega }_k}}SSMN\left(y_i;{\alpha }_k+x_i^T{\beta }_k,{\sigma }_k^2,{\lambda }_k,{\tau }_k;\mathcal{l}\right).$ (10)

Following the idea in [11], we can get the estimates of $\Psi$ by maximizing the penalized log-likelihood function which is defined as:

$L\left(\Psi \right)=l\left(\Psi \right)-p\left(\Psi \right),$ (11)

with the penalty function is given by:

$p\left(\Psi \right)=n{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\omega }_k}{\displaystyle \underset{j=1}{\overset{q}{\sum }}{p}_{a_k}}\left(\left|{\beta }_{kj}\right|\right),$ (12)

where $p_{a_k}(.)$ is a nonnegative and non-decreasing function in $\left|{\beta }_{kj}\right|$ with the turning parameter $a_k,k=1,\cdots ,K$. The turning parameter controls the intensity of the penalty for the regression coefficients.

The SCAD penalty has a type of oracle property as discussed in [8]. In this work, we complete the variable selection procedure using the following SCAD penalty function:

$p_{a_k}\left(\left|{\beta }_{kj}\right|\right)=\{\begin{array}{ll}{a}_k\left|{\beta }_{kj}\right|,\hfill & \left|{\beta }_{kj}\right|\le a_k\hfill \\ -\frac{\left({\beta }_{kj}^2-2c{a}_k\left|{\beta }_{kj}\right|+a_k^2\right)}{2\left(c-1\right)},\hfill & a_k<\left|{\beta }_{kj}\right|\le c{a}_k\hfill \\ \frac{a_k^2\left(c+1\right)}{2},\hfill & \left|{\beta }_{kj}\right|>c{a}_k\hfill \end{array}$ (13)

Meanwhile, inspired by [12], we propose a combined penalty function which balance the SCAD penalty and $l_2$ penalty. This penalty function by introducing a connection parameter $b$ is more effective in variable selection than directly mixing SCAD and $l_2$, and the specific form is given by:

$p_{a_k}\left(\left|{\beta }_{kj}\right|\right)=\{\begin{array}{ll}{a}_k\left[b\left|{\beta }_{kj}\right|+\left(1-b\right){\beta }_{kj}^2\right],\hfill & \left|{\beta }_{kj}\right|\le a_k\hfill \\ -\frac{b\left({\beta }_{kj}^2-2c{a}_k\left|{\beta }_{kj}\right|+a_k^2\right)}{2\left(c-1\right)}+a_k\left(1-b\right){\beta }_{kj}^2,\hfill & a_k<\left|{\beta }_{kj}\right|\le c{a}_k\hfill \\ \frac{a_k^{2}b\left(c+1\right)}{2}+a_k\left(1-b\right){\beta }_{kj}^2,\hfill & \left|{\beta }_{kj}\right|>c{a}_k\hfill \end{array}$ (14)

We call this new penalty function as MIXL2-SCAD. Some asymptotic properties of the penalty function are showed in [12], and the constant $a_k>0$ and $c>2$. Following the idea of [8], let $c=3.7$. In particular, the constant $b$, $0\le b\le 1$ and $a_k$ in MIXL2-SCAD jointly control the speed of contraction of ${\beta }_{kj}$, and when $b=1$, MIXL2-SCAD penalty reduces to the SCAD penalty.

4. Numeric Solutions

The expectation-maximization (EM) algorithm can be applied to mixture regression models based on SSMN distributions for maximizing the penalized log-likelihood function. When the M-step of EM is analytically intractable for SSMN distributions, it can be replaced with a sequence of conditional maximization (CM) steps which is derived from ECM algorithm [13]. Furthermore, we also maximize the constrained actual marginal log-likelihood function which called CML steps [14] for simplicity.

4.1. Maximization of the Penalized Log-Likelihood Function

Let us introduce the latent vector $Z_i={\left(z_{i1},\cdots ,z_{iK}\right)}^T$ with the component indicator variable $z_{ik}$ which has the following form:

$z_{ik}=\{\begin{array}{ll}1,\hfill & \text{the}i\text{th sample comes from the latent}k\text{th component,}\hfill \\ 0,\hfill & \text{otherwise}\text{.}\hfill \end{array}$ (15)

Using the Equations (2) and (4), we can get the following hierarchical representation for the mixture of SSMN distributions.

$\begin{array}{l}{Y}_i|\left(T_i=t_i,U_i=u_i,z_{ik}=1\right)~N\left({\alpha }_k+x_i^T{\beta }_k+\frac{\mathcal{l}\left(u_i\right){\sigma }_k{\lambda }_k}{{\left(1+{\lambda }_k^2\mathcal{l}\left(u_i\right)\right)}^{1/2}}{t}_i,\frac{{\sigma }_k^2\mathcal{l}\left(u_i\right)}{1+{\lambda }_k^2\mathcal{l}\left(u_i\right)}\right),\\ U_i|z_{ik}=1~H\left({\tau }_k\right),\\ T_i|z_{ik}=1~TN\left(0,1;\left(0,\infty \right)\right),\\ Z_i~M\left(1;{\omega }_1,\cdots ,{\omega }_K\right).\end{array}$ (16)

$TN\left(0,1;\left(0,\infty \right)\right)$ denotes the truncated normal distribution. Let $t={\left(t_1,\cdots ,t_n\right)}^T$, $u={\left(u_1,\cdots ,u_n\right)}^T$, $Y={\left(y_1,\cdots ,y_n\right)}^T$ and $Z^{\ast }={\left(Z_1,\cdots ,Z_n\right)}^T$. Among them, $t$ and $u$ are also regarded as latent vectors. Then the complete log-likelihood function with complete-data $Y_c={\left(Y^T,u^T,t^T,Z^{{\ast }^{\text{T}}}\right)}^T$ is given by:

$l_c\left(\Psi \right)=l_c\left(\omega \right)+l_c\left(\alpha ,\beta ,{\sigma }^2,\lambda ,{\tau }^{\ast }\right),$ (17)

with:

$l_c\left(\omega \right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{z}_{ik}}\log {\omega }_k,}$ (18)

$\begin{array}{c}{l}_c\left(\alpha ,\beta ,{\sigma }^2,\lambda ,{\tau }^{\ast }\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{z}_{ik}\left[C-\log {\sigma }_k^2-\frac{t_i^2}{2{\sigma }_k^2}+\frac{t_i{\lambda }_k}{{\sigma }_k^2}\left(y_i-{\alpha }_k-x_i^T{\beta }_k\right)\right]}}\\ -{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}\frac{z_{ik}}{2{\sigma }_k^2}\left\{\left[{\mathcal{l}}^{-1}\left(u_i\right)+{\lambda }_k^2\right]{\left(y_i-{\alpha }_k-x_i^T{\beta }_k\right)}^2\right\}}}\\ +{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{z}_{ik}\log h}}\left(u_i;{\tau }_k\right).\end{array}$ (19)

$C$ is a constant that does not depend on any unknown parameter, and $h\left(u_i;{\tau }_k\right)$ is the density function of the latent variable $u_i$.

Replacing $l\left(\Psi \right)$ with $l_c\left(\Psi \right)$ in the penalized log-likelihood function, the complete penalized log-likelihood function is given by:

$L_c\left(\Psi \right)=l_c\left(\Psi \right)-p\left(\Psi \right).$ (20)

Refer to the method of Fan and Li [8], given the initial parameter value ${\Psi }^{\left(0\right)}$, $p\left(\Psi \right)$ can be replaced by the following local quadratic function:

$p\left(\Psi \right)\approx n{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\omega }_k}{\displaystyle \underset{j=1}{\overset{q}{\sum }}\left[p_{a_k}\left(\left|{\beta }_{kj}^{\left(0\right)}\right|\right)+\frac{{p^{\prime }}_{a_k}\left(\left|{\beta }_{kj}^{\left(0\right)}\right|\right)}{2\left|{\beta }_{kj}^{\left(0\right)}\right|}\left({\beta }_{kj}^2-{\beta }_{kj}^{{\left(0\right)}^2}\right)\right]}$ (21)

This approximation will be applied in the CM-step of the algorithm at each iteration. The adjusted EM algorithm proceeds with the following three steps. The E-step calculates the conditional expectation of the complete penalized log-likelihood function, the CM-step and CML-step obtain the closed-form of parameter estimates.

· The E-step. Given the current estimates ${\hat{\Psi }}^{\left(m\right)}$, calculate the Q-function, $Q\left(\Psi |{\hat{\Psi }}^{\left(m\right)}\right)=E\left[L_c\left(\Psi \right)|Y,{\hat{\Psi }}^{\left(m\right)}\right]$, obtained as:

$\begin{array}{c}Q\left(\Psi |{\hat{\Psi }}^{\left(m\right)}\right)=Q_1\left(\omega |{\hat{\Psi }}^{\left(m\right)}\right)+Q_2\left(\alpha ,\beta ,{\sigma }^2,\lambda |{\hat{\Psi }}^{\left(m\right)}\right)\\ +Q_3\left({\tau }^{\ast }|{\hat{\Psi }}^{\left(m\right)}\right)-p\left(\Psi |{\hat{\Psi }}^{\left(m\right)}\right),\end{array}$ (22)

with:

$Q_1\left(\omega |{\hat{\Psi }}^{\left(m\right)}\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\hat{z}}_{ik}{}^{\left(m\right)}}\log {\omega }_k,}$ (23)

$\begin{array}{l}{Q}_2\left(\alpha ,\beta ,{\sigma }^2,\lambda |{\hat{\Psi }}^{\left(m\right)}\right)\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\hat{z}}_{ik}{}^{\left(m\right)}}\left[-\log {\sigma }_k^2-\frac{{\hat{t}}_{ik}^2{}^{\left(m\right)}}{2{\sigma }_k^2}+\frac{{\hat{t}}_{ik}{}^{\left(m\right)}{\lambda }_k}{{\sigma }_k^2}\left(y_i-{\alpha }_k-x_i^T{\beta }_k\right)\right]}\\ -{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}\frac{{\hat{z}}_{ik}{}^{(m)}}{2{\sigma }_k^2}\left[\left({\hat{\mathcal{l}}}_{ik}^{-1}{}^{\left(m\right)}+{\lambda }_k^2\right){\left(y_i-{\alpha }_k-x_i^T{\beta }_k\right)}^2\right],}}\end{array}$ (24)

$Q_3\left({\tau }^{\ast }|{\hat{\Psi }}^{\left(m\right)}\right)=E\left[{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{z}_{ik}\log h\left(u_i;{\tau }_k\right)}}|Y,{\hat{\Psi }}^{\left(m\right)}\right].$ (25)

The required expressions are ${\hat{z}}_{ik}{}^{\left(m\right)}$, ${\hat{t}}_{ik}^2{}^{\left(m\right)}$, ${\hat{t}}_{ik}{}^{\left(m\right)}$ and ${\hat{\mathcal{l}}}_{ik}^{-1}{}^{\left(m\right)}$.

First, the conditional expectation ${\hat{z}}_{ik}{}^{\left(m\right)}=E\left[z_{ik}|y_i,{\hat{\Psi }}^{\left(m\right)}\right]$ is given by:

${\hat{z}}_{ik}{}^{\left(m\right)}=\frac{{\hat{\omega }}_k{}^{\left(m\right)}SSMN\left(y_i;{\hat{\alpha }}_k{}^{\left(m\right)}+x_i^T{\hat{\beta }}_k{}^{\left(m\right)},{\hat{\sigma }}_k^2{}^{\left(m\right)},{\hat{\lambda }}_k{}^{\left(m\right)},{\hat{\tau }}_k{}^{\left(m\right)};\mathcal{l}\right)}{{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\hat{\omega }}_k{}^{\left(m\right)}SSMN\left(y_i;{\hat{\alpha }}_k{}^{\left(m\right)}+x_i^T{\hat{\beta }}_k{}^{\left(m\right)},{\hat{\sigma }}_k^2{}^{\left(m\right)},{\hat{\lambda }}_k{}^{\left(m\right)},{\hat{\tau }}_k{}^{\left(m\right)};\mathcal{l}\right)}}.$ (26)

Then, refer to [15], ${\hat{t}}_{ik}{}^{\left(m\right)}=E\left[t_i|y_i,{\hat{\Psi }}^{\left(m\right)},z_{ik}=1\right]$ and ${\hat{t}}_{ik}^2{}^{\left(m\right)}=E\left[t_i^2|y_i,{\hat{\Psi }}^{\left(m\right)},z_{ik}=1\right]$ can be evaluated by:

${\hat{t}}_{ik}{}^{\left(m\right)}={\hat{\lambda }}_k{}^{\left(m\right)}{\hat{e}}_{ik}{}^{\left(m\right)}+{\hat{\sigma }}_k{}^{\left(m\right)}{W}_{\Phi }\left(\frac{{\hat{\lambda }}_k{}^{\left(m\right)}{\hat{e}}_{ik}{}^{\left(m\right)}}{{\hat{\sigma }}_k{}^{\left(m\right)}}\right),$ (27)

${\hat{t}}_{ik}^2{}^{\left(m\right)}={\left({\hat{\lambda }}_k{}^{\left(m\right)}{\hat{e}}_{ik}{}^{\left(m\right)}\right)}^2+{\hat{\sigma }}_k^2{}^{\left(m\right)}+{\hat{\lambda }}_k{}^{\left(m\right)}{\hat{\sigma }}_k{}^{\left(m\right)}{\hat{e}}_{ik}{}^{\left(m\right)}{W}_{\Phi }\left(\frac{{\hat{\lambda }}_k{}^{\left(m\right)}{\hat{e}}_{ik}{}^{\left(m\right)}}{{\hat{\sigma }}_k{}^{\left(m\right)}}\right).$ (28)

with $W_{\Phi }\left(u\right)=\varphi \left(u\right)/\Phi \left(u\right)$ and ${\hat{e}}_{ik}{}^{\left(m\right)}=y_i-{\hat{\alpha }}_k{}^{\left(m\right)}-x_i^T{\hat{\beta }}_k{}^{\left(m\right)}$.

Further, ${\hat{\mathcal{l}}}_{ik}^{-1}{}^{\left(m\right)}$ has different expressions for RMR-SSMN models with different distributions in the SSMN family, obtained as:

${\hat{\mathcal{l}}}_{ik}^{-1}{}^{\left(m\right)}=\{\begin{array}{ll}1,\hfill & \text{forRMR-SNmodel}\hfill \\ \frac{{\hat{v}}_k{}^{\left(m\right)}+1}{{\hat{v}}_k{}^{\left(m\right)}+d_{ik}},\hfill & \text{forRMR-STNmodel}\hfill \\ \frac{1-{\hat{v}}_k{}^{\left(m\right)}+{\hat{v}}_k{}^{\left(m\right)}{\hat{\gamma }}_k{}^{{\left(m\right)}^{\frac{3}{2}}}\exp \left[\left(1-{\hat{\gamma }}_k{}^{\left(m\right)}\right)d_{ik}/2\right]}{1-{\hat{v}}_k{}^{\left(m\right)}+{\hat{v}}_k{}^{\left(m\right)}{\hat{\gamma }}_k{}^{{\left(m\right)}^{\frac{1}{2}}}\exp \left[\left(1-{\hat{\gamma }}_k{}^{\left(m\right)}\right)d_{ik}/2\right]},\hfill & \text{forRMR-SCNmodel}\hfill \\ {\hat{v}}_k{}^{\left(m\right)}{d}_{ik}{}^{{\hat{v}}_k{}^{\left(m\right)}-1},\hfill & \text{forRMR-SPEmodel}\hfill \end{array}$ (29)

with $d_{ik}={\left(y_i-{\hat{\alpha }}_k{}^{\left(m\right)}-x_i^T{\hat{\beta }}_k{}^{\left(m\right)}\right)}^2/{\hat{\sigma }}_k^2{}^{\left(m\right)}$.

· The CM-step. Maximize $Q\left(\Psi |{\hat{\Psi }}^{\left(m\right)}\right)$ with respect to $\Psi$ on the $\left(m+1\right)\text{th}$ iteration. As in [11], the mixing proportions are updated by:

${\hat{\omega }}_k{}^{\left(m+1\right)}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\hat{z}}_{ik}{}^{\left(m\right)}},$ (30)

which are the approximate iterated values. Maximizing $Q_1\left(\omega |{\hat{\Psi }}^{\left(m\right)}\right)$ with respect to the $\omega$ instead of maximizing $Q\left(\Psi |{\hat{\Psi }}^{\left(m\right)}\right)$ will simplify the computation of ${\hat{\omega }}_k{}^{\left(m+1\right)}$ and this updating scheme works well in our simulations.

We now consider that $\omega$ is constant, and maximize $Q\left(\Psi |{\hat{\Psi }}^{\left(m\right)}\right)$ with respect to the rest parameters in $\Psi$. The updates of ${\left({\alpha }_k,{\sigma }_k^2,{\lambda }_k,{\beta }_k\right)}^T$ are given by:

${\hat{\alpha }}_k{}^{\left(m+1\right)}=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\hat{z}}_{ik}{}^{\left(m\right)}\left[-{\hat{t}}_{ik}{}^{\left(m\right)}{\hat{\lambda }}_k{}^{\left(m\right)}+\left({\hat{\mathcal{l}}}_{ik}^{-1}{}^{\left(m\right)}+{\hat{\lambda }}_k{{}^{\left(m\right)}}^2\right)\left(y_i-x_i^T{\hat{\beta }}_k{}^{\left(m\right)}\right)\right]}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\hat{z}}_{ik}{}^{\left(m\right)}\left({\hat{\mathcal{l}}}_{ik}^{-1}{}^{\left(m\right)}+{\hat{\lambda }}_k{{}^{\left(m\right)}}^2\right)}},$ (31)

${\hat{\sigma }}_k^2{}^{\left(m+1\right)}=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\hat{z}}_{ik}{}^{\left(m\right)}\left[{\hat{t}}_{ik}^2{}^{\left(m\right)}-2{\hat{t}}_{ik}{}^{\left(m\right)}{\hat{\lambda }}_k{}^{\left(m\right)}{\hat{e}}_{ik}{}^{\left(m\right)}+\left({\hat{\mathcal{l}}}_{ik}^{-1}{}^{\left(m\right)}+{\hat{\lambda }}_k{{}^{\left(m\right)}}^2\right){\hat{e}}_{ik}{{}^{\left(m\right)}}^2\right]}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}2{\hat{z}}_{ik}{}^{\left(m\right)}}},$ (32)

${\hat{\lambda }}_k{}^{\left(m+1\right)}={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\hat{z}}_{ik}{}^{\left(m\right)}{\hat{t}}_{ik}{}^{\left(m\right)}{\hat{e}}_{ik}{}^{\left(m\right)}}/{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\hat{z}}_{ik}{}^{\left(m\right)}{\hat{e}}_{ik}{{}^{\left(m\right)}}^2},$ (33)

${\hat{\beta }}_k{}^{\left(m+1\right)}={\left[X^T{A}_{k}X+n{\hat{\sigma }}_k^2{}^{\left(m\right)}{\hat{\omega }}_k{}^{\left(m\right)}{\Delta }_a\left({\hat{\beta }}_k{}^{\left(m\right)}\right)\right]}^{-1}{X}^T{A}_k{B}_k,$ (34)

with:

$A_k=\left[\text{diag}({\hat{\mathcal{l}}}_{1k}^{-1}{}^{\left(m\right)},\cdots ,{\hat{\mathcal{l}}}_{nk}^{-1}{}^{\left(m\right)})+{\hat{\lambda }}_k{{}^{\left(m\right)}}^2{\mathbb{I}}_n\right]\text{diag}\left({\hat{z}}_{1k}{}^{\left(m\right)},\cdots ,{\hat{z}}_{nk}{}^{\left(m\right)}\right),$

$B_k={\left({\hat{b}}_{1k}{}^{\left(m\right)},\cdots ,{\hat{b}}_{nk}{}^{\left(m\right)}\right)}^T,{\hat{b}}_{ik}{}^{\left(m\right)}=y_i-{\hat{\alpha }}_k{}^{\left(m\right)}-\frac{{\hat{t}}_{ik}{}^{\left(m\right)}{\hat{\lambda }}_k{}^{\left(m\right)}}{{\hat{\mathcal{l}}}_{ik}^{-1}{}^{\left(m\right)}+{\hat{\lambda }}_k{{}^{\left(m\right)}}^2},$

${\Delta }_a\left({\hat{\beta }}_k{}^{\left(m\right)}\right)=\text{diag}\left(\frac{{p^{\prime }}_{a_k}\left(\left|{\hat{\beta }}_{k1}{}^{\left(m\right)}\right|\right)}{\left|{\hat{\beta }}_{k1}{}^{\left(m\right)}\right|},\frac{{p^{\prime }}_{a_k}\left(\left|{\hat{\beta }}_{k2}{}^{\left(m\right)}\right|\right)}{\left|{\hat{\beta }}_{k2}{}^{\left(m\right)}\right|},\cdots ,\frac{{p^{\prime }}_{a_k}\left(\left|{\hat{\beta }}_{kq}{}^{\left(m\right)}\right|\right)}{\left|{\hat{\beta }}_{kq}{}^{\left(m\right)}\right|}\right),$

and ${\mathbb{I}}_n$ is an identity matrix of order $n$, and $X={\left(x_1,\cdots ,x_n\right)}^T$ is a matrix of order $n\times q$.

· The CML-step. Fix ${\hat{\Psi }}_p{}^{\left(m+1\right)}={\left({\hat{\alpha }}_k{}^{\left(m+1\right)},{\hat{\beta }}_k{}^{\left(m+1\right)},{\hat{\sigma }}_k^2{}^{\left(m+1\right)},{\hat{\lambda }}_k{}^{\left(m+1\right)}\right)}^T$ and ${\hat{\omega }}_k{}^{\left(m+1\right)}$, update ${\tau }^{\ast }$ to get ${\hat{\tau }}^{\ast }{}^{\left(m+1\right)}={\left({\hat{\tau }}_1{}^{\left(m+1\right)},\cdots ,{\hat{\tau }}_K{}^{\left(m+1\right)}\right)}^T$ by optimizing the constrained log-likelihood function:

${\hat{\tau }}^{\ast }{}^{\left(m+1\right)}={\text{argmax}}_{{\tau }_1,\cdots ,{\tau }_K}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\log {\displaystyle \underset{k=1}{\overset{K}{\sum }}{\hat{\omega }}_k{}^{\left(m+1\right)}SSMN}}\left(y_i;{\hat{\Psi }}_p{}^{\left(m+1\right)},{\tau }_k;\mathcal{l}\right).$ (35)

The above iterations are repeated alternately until the maximum number of iterations is reached or a suitable stopping rule is met. In this work, the iterations will be completed when $\Vert {\hat{\Psi }}^{\left(m+1\right)}-{\hat{\Psi }}^{\left(m\right)}\Vert$ is sufficiently small, such as 10⁻⁵.

4.2. Selection of Turning Parameters and Components

When using the methods proposed in this paper, we also need to consider how to determine the components $K$ and the size of the turning parameters in the penalty function. Cross-Validation (CV), Generalized Cross-Validation (GCV), AIC and BIC are commonly used criteria for the selection of turning parameters.

As showed in [12], the final selected model will be overfitting if the turning parameter selected by GCV and they use the BIC to choose. In this paper, we also propose a proper BIC criterion for RMR-SSMN models to select turning parameters $a={\left(a_1,\cdots ,a_K\right)}^{\text{T}}$, the constant $b$ and the components $K$.

Let $\theta ={\left(a,b,K\right)}^{\text{T}}$, we should take a set of $\theta$ at a time over a suitable range and use the proposed adjusted EM algorithm to obtain the corresponding parameter estimates $\hat{\Psi }$. The optimal set of $\theta$ is selected by minimizing the following BIC criterion:

$BIC\left(\theta \right)=-2l\left(\hat{\Psi }\right)+\left(\tilde{p}K-1+{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\eta }_k}\right)\times \log \left(n\right).$ (36)

where ${\eta }_k$ represents the number of non-zero regression coefficients of ${\beta }_k$ and $\tilde{p}$ is either equal to 4 (RMR-SN model), 5 (RMR-STN and RMR-SPE models) or 6 (RMR-SCN model).

5. Simulation Studies

We perform Monte Carlo simulations to evaluate the performance of the proposed robust mixture model and adjusted EM algorithm. To evaluate the effect of variable selection and the accuracy of parameter estimates, we use the correctly estimated zero coefficients (S1), correctly estimated non-zero coefficients (S2), the mean estimate over all falsely identified non-zero predictors ( $M_{NZ}$ ) [16] of $\beta$ and the mean squared error (MSE) of regression coefficients ( $\text{MSE}\left(\hat{\beta }\right)$ ),

$\text{MSE}\left(\hat{\beta }\right)=E{\left({\hat{\beta }}_k-{\beta }_k\right)}^T\left({\hat{\beta }}_k-{\beta }_k\right)$.

5.1. Simulation 1

The first simulation uses the SCAD penalty function to select significant variables for RMR-STN, RMR-SPE and RMR-SCN models, and compare the simulation results with the Gaussian FMR model and RMR-SN model.

We set $K=2$ for the simulation so that the sample data set $\left\{\left(y_i,x_i\right);i=1,\cdots ,n\right\}$ for the mixture regression model is derived from the following model:

$y=\{\begin{array}{ll}{\alpha }_1+x^T{\beta }_1+{\epsilon }_1,\hfill & Z=1\hfill \\ {\alpha }_2+x^T{\beta }_2+{\epsilon }_2,\hfill & Z=2\hfill \end{array}$ (37)

where $Z$ is used to identify the subgroup that the sample belongs to. ${\alpha }_1=2$, ${\beta }_1={\left(4,1,-2,0,0,0,0,0,0,0\right)}^T$, ${\alpha }_2=-2$, ${\beta }_2={\left(1,-3,0,2,3,0,0,0,0,0\right)}^T$, ${\omega }_1=0.6$ and ${\omega }_2=0.4$.

The covariate $x$ is generated from a multivariate normal with mean 0, variance 1, and two correlation structures: ${\rho }_{ij}={0.5}^{\left|i-j\right|}$, $1\le i,j\le n$. The simulation considers the following three error distributions cases: 1) the random errors ${\epsilon }_1$ and ${\epsilon }_2$ follow the t-distribution with 3 degrees of freedom ( $t\left(3\right)$ ); 2) the random errors ${\epsilon }_1$ and ${\epsilon }_2$ follow the chi-square distribution with 3 degrees of freedom ( ${\chi }^2\left(3\right)$ ); 3) The random errors ${\epsilon }_1$ and ${\epsilon }_2$ follow the mixture distribution of normal $0.9N\left(0,1\right)+0.1N\left(0,5^2\right)$. So that there are 15 sets of combination, and for each combination, we respectively performed 100 repetitions for the simulation with $n=300$.

From Table 1, the value of S1 in Com1 and Com2 from RMR-STN, RMR-SPE and RMR-SCN are all bigger than the value in FMR model for all three cases, respectively. In case (1), the S1 in Com2 from RMR-SPE is biggest (S1 = 0.9533), however, the S1 in Com2 from FMR is smallest (S1 = 0.8533). In case (2), the RMR-SCN model has the biggest S1 (S1 = 0.9033) in Com2 while the S1 in Com2 from FMR is 0.8167. In case (3), when the S2 in Com1 and Com2 from FMR are 0.9933 and 0.9950, respectively, the values of S2 in both components from RMR-STN are 1.00.

Furthermore, the value of $\text{MSE}\left(\hat{\beta }\right)$ in Com1 and Com2 from RMR-STN, RMR-SPE and RMR-SCN are much smaller than the value in FMR model. When errors follow ${\chi }^2\left(3\right)$ distribution, RMR-SN performs well with the smallest $\text{MSE}\left(\hat{\beta }\right)$ and S2 = 1.00 in both Com1 and Com2 that indicates the non-zero coefficients are all identified correctly. Overall, RMR-SSMN models are more robust than FMR for variable selection when the data set is asymmetric ( ${\chi }^2\left(3\right)$ ), heavy-tailed ( $t\left(3\right)$ ) and contaminated ( $0.9N\left(0,1\right)+0.1N\left(0,5^2\right)$ ).

5.2. Simulation 2

Simulation 2 uses the MIXL2-SCAD penalty function to select significant variables for RMR-STN, RMR-SPE and RMR-SCN models. By comparing the results of Simulation 1 and Simulation 2, the effects of SCAD and MIXL2-SCAD penalty function on variable selection are analyzed. In addition, the generation of the sample data set and the distributions of random errors in this simulation are the same as in Simulation 1, and both $n=300$ and $n=500$ cases are considered.

Table 1. Results of FMR, RMR-SN, RMR-STN, RMR-SPE and RMR-SCN using SCAD penalty function on 100 replicates.

From Table 2, we can know that as the sample size $n$ increases, the values of S1 and S2 in Com1 and Com2 are getting closer and closer to 1, and the value of $\text{MSE}\left(\hat{\beta }\right)$ is getting smaller and smaller, indicating the asymptotic property of parameter estimates. When $n=500$ and errors follow $t\left(3\right)$ distribution, the values of S1 and S2 in Com1 from RMR-SPE model are equal to 1.00, which indicates that the MIXL2-SCAD penalty ensures the non-zero and zero coefficients

Table 2. Results of RMR-STN, RMR-SPE and RMR-SCN using MIXL2-SCAD penalty function on 100 replicates.

can be identified completely. When $n=500$ and errors follow $0.9N\left(0,1\right)+0.1N\left(0,5^2\right)$, the same result appears in Com2 from RMR-SPE and RMR-SCN model. The absolute values of mean estimate over all falsely identified non-zero predictors ( $M_{NZ}$ ) are smaller than 0.01 from MIXL2-SCAD when $n=500$.

By comparing Table 1 and Table 2, we can see that the values of S1 and S2 in Com1 and Com2 from MIXL2-SCAD are all greater than or equal to the values from SCAD penalty for all cases when $n=300$. It is worth noting that in case (3), when n = 300, the values of S2 in Com1 and Com2 from RMR-STN, RMR-SPE and RMR-SCN using MIXL2-SCAD penalty are all 1.00, however, the values of S2 in Com1 from RMR-SPE and RMR-SCN using SCAD penalty are 0.9933. From these comparisons of experimental data, we can know that MIXL2-SCAD performs better than SCAD penalty in variable selection.

6. Real Data Analysis

In this section, we obtain the Seoul bike sharing demand data set from the website http://archive.ics.uci.edu/ml/datasets.php. From this dataset, we screen out the total number of bikes rented from 10:00 am to 11:00 am every functional day of bike rental system in Seoul from December 1, 2017 to November 30, 2018 with 12 features that may affect the demand of rental bikes. There are 353 observations in total. The 12 features are: temperature ( $x_1$ ), humidity ( $x_2$ ), wind-speed ( $x_3$ ), visibility ( $x_4$ ), dew point temperature ( $x_5$ ), solar radiation ( $x_6$ ), rainfall ( $x_7$ ), snowfall ( $x_8$ ), holiday (holiday = 1, else = 0; $x_9$ ), spring (spring = 1, else = 0; $x_{10}$ ), summer (summer = 1, else = 0; $x_{11}$ ) and autumn (autumn = 1, else = 0; $x_{12}$ ). $x_9\text{-}{x}_{12}$ are dummy variables and $x_{10}\text{-}{x}_{12}$ indicate different seasons. Considering that there may be further differential effects between seasons and holiday, we continue to introduce 3 interaction terms between dummy variables, namely $x_9\ast x_{10}$, $x_9\ast x_{11}$, $x_9\ast x_{12}$. This leads to a set of 15 potential covariates affecting rented bike count (RBC) from 10:00 am to 11:00 am.

Let $Y=\text{RBC}/sd\left(\text{RBC}\right)$ be the response variable, where $sd\left(\text{RBC}\right)$ is the standard deviation of RBC. Figure 1 shows the histogram and density estimate of $Y$, we can see that the data set has obvious heterogeneity, so that the RMR-STN model is applicable. We also apply RMR-SPE and RMR-SCN models to this real data set, the outcomes are worse than RMR-STN’s result, thus we do not report the results here.

The parameter estimates under FMR, RMR-STN ( $K=2$ ) and RMR-STN ( $K=3$ ) with BIC method and MIXL2-SCAD penalty function are given in Table 3. The $K=3$ RMR-STN model has the lowest BIC (542.5) and the $K=2$ RMR-STN model ranks second (BIC = 544.7) when FMR model has the biggest BIC (562.8). Furthermore, the predicted rented bike count from the $K=3$ RMR-STN model has the smallest MSE of 0.09 and the biggest regression ${\tilde{R}}^2$ of 0.90.

From Table 3, the bike rented demand can be divided into three categories: “low”, “medium” and “high” during the time period from 10:00 am to 11:00 am

Figure 1. Histogram and density estimate for $Y=\text{RBC}/sd\left(\text{RBC}\right)$.

Table 3. Summary of FMR, RMR-STN ( $K=2$ ) and RMR-STN ( $K=3$ ) model with BIC method and MIXL2-SCAD penalty for Seoul bike sharing demand data set.

with $K=3$ RMR-STN model. Humidity is a negative factor for all three types of demand. When the bike rented demand is “medium”, warmer temperature and increased solar radiation help increase bike demand, while rainfall, snowfall, and holidays reduce the demand. In contrast, when bike rented demand is “high”, the positive effect of dew point temperature on bike demand is greatest, while the negative effects of holidays and snowfall disappear. In addition, we can also find that the rented bike count has a strong seasonality and the rented count will be more in other seasons than in winter.

7. Conclusion

In this paper, we mainly propose a robust mixture regression model based on the skew scale mixtures of normal distributions (RMR-SSMN) which can avoid the potential limitation of normal mixtures. A new penalty function (MIXL2-SCAD) which combines SCAD and $l_2$ penalties is presented for variable selection. Through simulations, we find that the RMR-SSMN models are more robust than general FMR models for heterogeneous data with asymmetry and heavy-tailed properties, and outliers. Furthermore, the capability of MIXL2-SCAD to select the most parsimonious FMR model is obviously better than SCAD. The proposed methodology is applied to a real data set and achieves reasonable results. However, this paper only focuses on the mixture of the simple linear model, and further research can focus on the mixture of the semiparametric model or nonparametric model.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Krykun, I. (2018) The Arc-Sine Laws for the Skew Brownian Motion and Their Interpretation. Journal of Applied Mathematics and Physics, 6, 347-357. https://doi.org/10.4236/jamp.2018.62033
[2]	McLachlan, G.J. and Peel, D. (2000) Finite Mixture Models. Wiley, New York. https://doi.org/10.1002/0471721182
[3]	McLachlan, G.J. and Peel, D. (1998) Robust Cluster Analysis via Mixtures of Multivariate T-Distributions. Joint IAPR International Workshops on Statistical Techniques in Pattern Recognition (SPR) and Structural and Syntactic Pattern Recognition (SSPR), Sydney, 11-13 August 1998, 658-666. https://doi.org/10.1007/BFb0033290
[4]	Basso, R.M., Lachos, V.H., Cabral, C.R.B. and Ghosh, P. (2011) Robust Mixture Modeling based on Scale Mixtures of Skew-Normal Distributions. Computational Statistics & Data Analysis, 54, 2926-2941. https://doi.org/10.1016/j.csda.2009.09.031
[5]	Franczak, B.C., Browne, R.P. and McNicholas, P.D. (2014) Mixtures of Shifted Asymmetric Laplace Distributions. IEEE Transactions on Pattern Analysis and Machine Intelligence, 36, 1149-1157. https://doi.org/10.1109/TPAMI.2013.216
[6]	Zou, H. and Hastie, T. (2005) Regularization and Variable Selection via Elastic Net. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67, 301-320. https://doi.org/10.1111/j.1467-9868.2005.00503.x
[7]	Zhang, C.H. (2010) Nearly Unbiased Variable Selection under Minimax Concave Penalty. The Annals of Statistics, 38, 894-942. https://doi.org/10.1214/09-AOS729
[8]	Fan, J. and Li, R. (2001) Variable Selection via Nonconcave Penalized Likelihood and Its Oracle Properties. Journal of the American Statistical Association, 96, 1348-1360. https://doi.org/10.1198/016214501753382273
[9]	Ferreira, C.S., Bolfarine, H. and Lachos, V.H. (2011) Skew Scale Mixtures of Normal Distributions: Properties and Estimation. Statistical Methodology, 8, 154-171. https://doi.org/10.1016/j.stamet.2010.09.001
[10]	Gomez, H.W., Venegas, O. and Bolfarine, H. (2007) Skew-Symmetric Distributions Generated by the Normal Distribution Function. Environmetrics, 18, 395-407. https://doi.org/10.1002/env.817
[11]	Khalili, A. and Chen, J. (2007) Variable Selection in Finite Mixture of Regression Models. Journal of the American Statistical Association, 102, 1025-1038. https://doi.org/10.1198/016214507000000590
[12]	Khalili, A. (2010) New Estimation and Feature Selection Methods in Mixture-of-Experts Models. Canadian Journal of Statistics, 38, 519-539. https://doi.org/10.1002/cjs.10083
[13]	Meng, X.L. and Rubin, D.B. (1993) Maximum Likelihood Estimation via the ECM Algorithm: A General Framework. Biometrika, 80, 267-278. https://doi.org/10.1093/biomet/80.2.267
[14]	Liu, C. and Rubin, D.B. (1994) A Simple Extension of EM and ECM with Faster Monotone Convergence. Biometrika, 81, 633-648. https://doi.org/10.1093/biomet/81.4.633
[15]	Ferreira, C.S. and Lachos, V.H. (2016) Nonlinear Regression Models under Skew Scale Mixtures of Normal Distributions. Statistical Methodology, 33, 131-146. https://doi.org/10.1016/j.stamet.2016.08.004
[16]	Lloyd-Jones, L.R., Nguyen, H.D. and McLachlan, G.J. (2018) A Globally Convergent Algorithm for Lasso-Penalized Mixture of Linear Regression Models. Computational Statistics and Data Analysis, 119, 19-38. https://doi.org/10.1016/j.csda.2017.09.003