Adaptive Sparse Group Variable Selection for a Robust Mixture Regression Model Based on Laplace Distribution

Jiangtao Wang, Wanzhou Ye

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

DOI: 10.4236/apm.2020.101004                

Abstract

The traditional estimation of Gaussian mixture model is sensitive to heavy-tailed errors; thus we propose a robust mixture regression model by assuming that the error terms follow a Laplace distribution in this article. And for the variable selection problem in our new robust mixture regression model, we introduce the adaptive sparse group Lasso penalty to achieve sparsity at both the group-level and within-group-level. As numerical experiments show, compared with other alternative methods, our method has better performances in variable selection and parameter estimation. Finally, we apply our proposed method to analyze NBA salary data during the period from 2018 to 2019.

Keywords

Robust Mixture Regression, Laplace Distribution, Adaptive Sparse Group Lasso

1. Introduction

The mixture regression model is a powerful tool to explain the relationships between the response variable and the covariates when the population is heterogeneous and consists of several homogeneous components, and the early research can trace back to [1]. In 1977, EM algorithm was first proposed by [2] ; it greatly simplified the solution procedure of the mixture regression model. Then the mixture regression model attracted a lot of interest from statisticians; it was widely applied in many fields, such as business, marketing and social sciences. See [3] [4] [5] for example.

Recently, the research about the mixture regression model is becoming more and more detailed. On the one hand, statisticians paid attention to improving the robustness of mixture regression model, [6] used the t-distribution for overcoming the influence of outliers and [7] introduced a robust mixture regression model by assuming the error terms follow a Laplace distribution. Further, Wu et al. [8] dropped any parametric assumption about the error densities and proposed the mixture of quantile regressions model. On the other hand, variable selection became a research hotspot in mixture regression modeling. Khalili and Chen [9] considered a class of penalization methods, including L₁-norm penalty, SCAD penalty and Hard penalty. Furthermore, the adaptive Lasso was introduced as a penalty function in [10], and [11] suggested the use of the Lasso-penalized mixture regression model as a screening mechanism in a two-step procedure.

However, the above regularization methods are more focused on individual variable selection. They all ignore the grouping structures which describe the inherent interconnections among predictors. It may lead to inefficient models. In order to achieve both the robustness of mixture regression model and correct identification of group structures, we assume random errors follow a Laplace distribution and consider a situation that covariates have natural grouping structures, where those in the same group are correlated. In this case, variable selection should be conducted at both the group-level and within-group-level; thus we use the adaptive sparse group Lasso [12] as a penalty function of our proposed mixture regression model and adopt EM algorithm to estimate the mixture regression parameters. Moreover, under some mild conditions, we can prove that the maximum penalized log-likelihood estimators are both sparse and $\sqrt{n}$ -consistent simultaneously.

The rest of this article is organized as follows. In Section 2, we introduce the robust mixture regression model based on Laplace distribution and adopt the adaptive sparse group Lasso for variable selection. In Section 3, we prove some asymptotic properties for our proposed method. In Section 4, we solve the problem of tuning parameters and components selection. Section 5 conducts a numerical simulation to evaluate the performance of our method. In Section 6, we apply our proposed method to NBA salary data. Finally, the conclusion of this paper is given in Section 7.

2. Model Overview

2.1. Robust Mixture Regression with Laplace Distribution

Let $\left\{\left(x_1\mathrm{,}{y}_1\right)\mathrm{,}\cdots \mathrm{,}\left(x_n\mathrm{,}{y}_n\right)\right\}$ be a random sample of observations from the population $\left(x\mathrm{,}y\right)$, where $x={\left(x_1\mathrm{,}\cdots \mathrm{,}{x}_p\right)}^{\text{T}}\in {ℝ}^p$ is a p-dimensional covariate vector, and $y\in ℝ$ is the response variable which is dependent on corresponding $x$. Furthermore, for g mixture components, we can say that $\left(x\mathrm{,}y\right)$ follows a mixture linear regression model based on a normal distribution if the conditional density of y given $x$ is

$f\left(y\mathrm{<mo>\; |\; </mo>}x\mathrm{,}\Psi \right)={\displaystyle \underset{j=1}{\overset{g}{\sum }}}\frac{{\pi }_j}{\sqrt{2\pi {\sigma }_j^2}}\exp \left(-\frac{{\left(y-{\alpha }_j-x^{\text{T}}{\beta }_j\right)}^2}{2{\sigma }_j^2}\right)\mathrm{,}$ (1)

where the mixing probabilities satisfy ${\displaystyle {\sum }_{j=1}^g}{\pi }_j=1$, ${\pi }_j>0$ and the parameter $\Psi =\left({\pi }_1\mathrm{,}{\alpha }_1\mathrm{,}{\beta }_1\mathrm{,}{\sigma }_1^2\mathrm{,}\cdots \mathrm{,}{\pi }_g\mathrm{,}{\alpha }_g\mathrm{,}{\beta }_g\mathrm{,}{\sigma }_g^2\right)$. For the jth mixture component, there are intercept ${\alpha }_j\in ℝ$, regression coefficients ${\beta }_j={\left({\beta }_{j1}\mathrm{,}\cdots \mathrm{,}{\beta }_{jp}\right)}^{\text{T}}$ and variance ${\sigma }_j^2\left({\sigma }_j>0\right)$.

It is known that the mixture linear regression model is sensitive to outliers or heavy-tailed error distributions, and outliers impact more heavily on the mixture linear regression model than on the usual linear regression model, since outliers not only affect the estimation of the regression parameters, but also possibly totally blur the mixture structure. In order to improve the robustness of the estimation procedure, we introduce a robust mixture regression model with a Laplace distribution

$f\left(y\mathrm{<mo>\; |\; </mo>}x\mathrm{,}\Psi \right)={\displaystyle \underset{j=1}{\overset{g}{\sum }}}\frac{{\pi }_j}{\sqrt{2{\sigma }_j^2}}\exp \left(-\frac{\sqrt{2}\left|y-{\alpha }_j-x^{\text{T}}{\beta }_j\right|}{{\sigma }_j}\right)\mathrm{.}$ (2)

Then we can estimate the unknown parameter $\Psi$ by maximizing the log-likelihood function

$L_n\left(\Psi \right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log {\displaystyle \underset{j=1}{\overset{g}{\sum }}}\frac{{\pi }_j}{\sqrt{2{\sigma }_j^2}}\exp \left(-\frac{\sqrt{2}\left|y_i-{\alpha }_j-x_i^{\text{T}}{\beta }_j\right|}{{\sigma }_j}\right)\mathrm{.}$ (3)

2.2. Adaptive Sparse Group Lasso for Variable Selection

Now, we consider a situation that covariates have natural grouping structures and can be divided into K groups as $x={\left(x^{\left[1\right]\text{T}}\mathrm{,}\cdots \mathrm{,}{x}^{\left[k\right]\text{T}}\mathrm{,}\cdots \mathrm{,}{x}^{\left[K\right]\text{T}}\right)}^{\text{T}}\in {ℝ}^p$ by some known rules, where $x^{\left[k\right]}={\left(x_1^{\left[k\right]}\mathrm{,}\cdots \mathrm{,}{x}_{p_k}^{\left[k\right]}\right)}^{\text{T}}$ is a group which contains $p_k$ variables and ${\displaystyle {\sum }_{k=1}^K}{p}_k=p$. Then, the log-likelihood function can be written as

$L_n\left(\Psi \right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log {\displaystyle \underset{j=1}{\overset{g}{\sum }}}\frac{{\pi }_j}{\sqrt{2{\sigma }_j^2}}\exp \left(-\frac{\sqrt{2}\left|y_i-{\alpha }_j-{\displaystyle {\sum }_{k=1}^K{x}_i^{\left[k\right]\text{T}}{\beta }_j^{\left[k\right]}}\right|}{{\sigma }_j}\right)\mathrm{,}$ (4)

where ${\beta }_j^{\left[k\right]}={\left({\beta }_{j1}^{\left[k\right]}\mathrm{,}\cdots \mathrm{,}{\beta }_{j{p}_k}^{\left[k\right]}\right)}^{\text{T}}$.

In order to exploit the grouping structures of covariates, we apply the adaptive sparse group Lasso (adaSGL) to the robust mixture regression model, the penalized log-likelihood function

$F_n\left(\Psi \right)=L_n\left(\Psi \right)-P_n\left(\Psi \right)\mathrm{.}$ (5)

Here

$P_n\left(\Psi \right)=n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{t=1}{\overset{p}{\sum }}}{\lambda }_{jt1}\left|{\beta }_{jt}\right|+n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{k=1}{\overset{K}{\sum }}}{\lambda }_{jk2}\Vert {\beta }_j^{\left[k\right]}\Vert \mathrm{,}$ (6)

where $\Vert \mathrm{.}\Vert$ represents the Euclidean norm, ${\lambda }_{jt1}={\lambda }_{j1}{\omega }_{jt}$ and ${\lambda }_{jk2}={\lambda }_{j2}{\xi }_{jk}$. Moreover, the weights are defined based on the maximum penalized log-likelihood estimator $\tilde{\Psi }$ when $P_n\left(\Psi \right)$ is a Lasso penalty,

${\omega }_{jt}={\left|{\tilde{\beta }}_{jt}\right|}^{-1},\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}{\xi }_{jk}={\Vert {\tilde{\beta }}_j^{\left[k\right]}\Vert }^{-1}.$ (7)

Next, we follow the approach of Hunter and Li [13] and consider to maximize the $\epsilon$ -approximate penalized log-likelihood function

$F_{n\mathrm{,}\epsilon }\left(\Psi \right)=L_n\left(\Psi \right)-P_{n\mathrm{,}\epsilon }\left(\Psi \right)\mathrm{.}$ (8)

Here

$P_{n\mathrm{,}\epsilon }\left(\Psi \right)=n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{t=1}{\overset{p}{\sum }}}{\lambda }_{jt1}\sqrt{{\beta }_{jt}^2+{\epsilon }^2}+n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{k=1}{\overset{K}{\sum }}}{\lambda }_{jk2}\sqrt{{\displaystyle \underset{t=1}{\overset{p_k}{\sum }}}{\beta }_{jt}^{2\left[k\right]}+{\epsilon }^2}$ (9)

for some small $\epsilon >0$, and the weights are

${\omega }_{jt}={\left({\tilde{\beta }}_{jt}^2+{\epsilon }^2\right)}^{-1/2},\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}{\xi }_{jk}={\left({\displaystyle \underset{t=1}{\overset{p_k}{\sum }}}{\tilde{\beta }}_{jt}^{2\left[k\right]}+{\epsilon }^2\right)}^{-1/2}.$ (10)

Following Hunter and Li [13], we can similarly show that $\left|F_{n\mathrm{,}\epsilon }\left(\Psi \right)-F_n\left(\Psi \right)\right|\to 0$ uniformly as $\epsilon \to 0$, over any compact subset of the parameter space.

2.3. EM Algorithm for Robust Mixture Regression

However, the above penalized log-likelihood does not have an explicit maximizer. We introduce an EM algorithm to simplify the computation and denote $Z_{ij}$ as a latent Bernoulli variable such that

$Z_{ij}=(\begin{array}{cc}\mathrm{1,}& \text{if}i\text{th}\text{observation}\left(x_i,y_i\right)\text{is}\text{from}j\text{th}\text{component}\mathrm{<mo>\; ;\; </mo>}\\ \mathrm{0,}& \text{otherwise}\mathrm{.}\end{array}$ (11)

If the complete data set $T=\left(X\mathrm{,}Y\mathrm{,}Z\right)$ is observable, the complete log-likelihood function is

$L_n\left(\Psi \mathrm{<mo>\; ;\; </mo>}T\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{Z}_{ij}\log \left[\frac{{\pi }_j}{\sqrt{2{\sigma }_j^2}}\exp \left(-\frac{\sqrt{2}\left|y_i-{\alpha }_j-{\displaystyle {\sum }_{k=1}^K{x}_i^{\left[k\right]\text{T}}{\beta }_j^{\left[k\right]}}\right|}{{\sigma }_j}\right)\right]\mathrm{,}$ (12)

where $X={\left(x_1\mathrm{,}\cdots \mathrm{,}{x}_n\right)}^{\text{T}}$, $Y=\left(y_1\mathrm{,}\cdots \mathrm{,}{y}_n\right)$ and $Z=\left(Z_{11}\mathrm{,}\cdots \mathrm{,}{Z}_{ng}\right)$.

According to Andrews and Mallows [14], we know that a Laplace distribution can be expressed as a mixture of a normal distribution and another distribution related to the exponential distribution. To be specific, there are latent scale variables $V=\left(V_1\mathrm{,}\cdots \mathrm{,}{V}_n\right)$ such that we have the complete log-likelihood function

$\begin{array}{c}{L}_n\left(\Psi \mathrm{<mo>\; ;\; </mo>}D\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{Z}_{ij}\log \left\{\frac{{\pi }_j{V}_i}{\sqrt{\pi {\sigma }_j^2}}\exp \left[-\frac{V_i^2{\left(y_i-{\alpha }_j-{\displaystyle {\sum }_{k=1}^K{x}_i^{\left[k\right]\text{T}}{\beta }_j^{\left[k\right]}}\right)}^2}{{\sigma }_j^2}\right]\right\}\\ +{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{Z}_{ij}\log \left[\frac{1}{V_i^3}\exp \left(-\frac{1}{2V_i^2}\right)\right],\end{array}$ (13)

where $D=\left(T\mathrm{,}V\right)$. Naturally, we can obtain the penalized complete log-likelihood $F_{n\mathrm{,}\epsilon }\left(\Psi \mathrm{<mo>\; ;\; </mo>}D\right)=L_n\left(\Psi \mathrm{<mo>\; ;\; </mo>}D\right)-P_{n\mathrm{,}\epsilon }\left(\Psi \right)$.

Suppose that ${\Psi }^{\left(r\right)}$ is a parameter estimate for the rth iteration. In E step of EM algorithm, we can get $E\left[F_{n\mathrm{,}\epsilon }\left(\Psi \mathrm{<mo>\; ;\; </mo>}D\right)\mathrm{<mo>\; |\; </mo>}S\mathrm{,}{\Psi }^{\left(r\right)}\right]$ by calculating

${\tau }_{ij}^{\left(r\right)}=E\left[Z_{ij}\mathrm{<mo>\; |\; </mo>}S\mathrm{,}{\Psi }^{\left(r\right)}\right]\mathrm{,\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}{\delta }_{ij}^{\left(r\right)}=E\left[V_i^2\mathrm{<mo>\; |\; </mo>}S\mathrm{,}{\Psi }^{\left(r\right)}\mathrm{,}{Z}_{ij}=1\right]\mathrm{,}$ (14)

where $S=\left(X\mathrm{,}Y\right)$. And we can show that

${\tau }_{ij}^{\left(r\right)}=\frac{{\pi }_j^{\left(r\right)}{\sigma }_j^{-1\left(r\right)}\exp \left(-\sqrt{2}\left|y_i-{\alpha }_j^{\left(r\right)}-{\displaystyle {\sum }_{k=1}^K{x}_i^{\left[k\right]\text{T}}{\beta }_j^{\left[k\right]\left(r\right)}}\right|{\sigma }_j^{-1\left(r\right)}\right)}{{\displaystyle {\sum }_{j=1}^g{\pi }_j^{\left(r\right)}{\sigma }_j^{-1\left(r\right)}\exp \left(-\sqrt{2}\left|y_i-{\alpha }_j^{\left(r\right)}-{\displaystyle {\sum }_{k=1}^K{x}_i^{\left[k\right]\text{T}}{\beta }_j^{\left[k\right]\left(r\right)}}\right|{\sigma }_j^{-1\left(r\right)}\right)}}$ (15)

and

${\delta }_{ij}^{\left(r\right)}=\frac{{\sigma }_j^{\left(r\right)}}{\sqrt{2}\left|y_i-{\alpha }_j^{\left(r\right)}-{\displaystyle {\sum }_{k=1}^K{x}_i^{\left[k\right]\text{T}}{\beta }_j^{\left[k\right]\left(r\right)}}\right|}.$ (16)

The calculation for ${\delta }_{ij}^{\left(r\right)}$ follows the same argument as in Phillips [15].

In M step, we will maximize $E\left[F_{n\mathrm{,}\epsilon }\left(\Psi \mathrm{<mo>\; ;\; </mo>}D\right)\mathrm{<mo>\; |\; </mo>}S\mathrm{,}{\Psi }^{\left(r\right)}\right]$ for updating $\Psi$. Now, we follow the tactic of [16] and find a local quadratic approximation of $\sqrt{{\psi }^2+{\epsilon }^2}$,

$\sqrt{{\psi }^2+{\epsilon }^2}\simeq \sqrt{{\psi }_0^2+{\epsilon }^2}+\frac{{\psi }^2-{\psi }_0^2}{2\sqrt{{\psi }_0^2+{\epsilon }^2}}$ (17)

in a neighborhood of ${\psi }_0$. Then, we can replace the penalty function $P_{n\mathrm{,}\epsilon }\left(\Psi \right)$ in $\left(r+1\right)\text{th}$ iteration by

$\begin{array}{c}{\tilde{P}}_{n\mathrm{,}\epsilon }\left(\Psi \mathrm{<mo>\; ;\; </mo>}{\Psi }^{\left(r\right)}\right)=n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{t=1}{\overset{p}{\sum }}}{\lambda }_{jt1}\left[{\eta }_{jt}^{\left(r\right)}+\frac{{\beta }_{jt}^2-{\beta }_{jt}^{2\left(r\right)}}{2{\eta }_{jt}^{\left(r\right)}}\right]\\ +n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{k=1}{\overset{K}{\sum }}}{\lambda }_{jk2}\left[{\gamma }_{jk}^{\left(r\right)}+\frac{{\displaystyle {\sum }_{t=1}^{p_k}\left({\beta }_{jt}^{2\left[k\right]}-{\beta }_{jt}^{2\left[k\right]\left(r\right)}\right)}}{2{\gamma }_{jk}^{\left(r\right)}}\right],\end{array}$ (18)

where ${\eta }_{jt}^{\left(r\right)}=\sqrt{{\beta }_{jt}^{2\left(r\right)}+{\epsilon }^2}$ and ${\gamma }_{jk}^{\left(r\right)}=\sqrt{{\displaystyle {\sum }_{t=1}^{p_k}}{\beta }_{jt}^{2\left[k\right]\left(r\right)}+{\epsilon }^2}$. Similarly, from Lange [17], we have

${\left(y-\alpha -x^{\text{T}}\psi \right)}^2\simeq \frac{1}{p}{\displaystyle \underset{t=1}{\overset{p}{\sum }}}{\left[y-\alpha -x^{\text{T}}{\psi }_0-p{x}_t\left({\psi }_t-{\psi }_{0t}\right)\right]}^2$ (19)

in a neighborhood of ${\psi }_0={\left({\psi }_{01}\mathrm{,}\cdots \mathrm{,}{\psi }_{0p}\right)}^{\text{T}}$, where p is the dimensionality of $\psi$. And we apply (19) to $E\left[L_{n\mathrm{,}\epsilon }\left(\Psi \mathrm{<mo>\; ;\; </mo>}D\right)\mathrm{<mo>\; |\; </mo>}S\mathrm{,}{\Psi }^{\left(r\right)}\right]$, there is a local approximation $Q\left(\Psi \mathrm{<mo>\; ;\; </mo>}{\Psi }^{\left(r\right)}\right)$ of $E\left[F_{n\mathrm{,}\epsilon }\left(\Psi \mathrm{<mo>\; ;\; </mo>}D\right)\mathrm{<mo>\; |\; </mo>}S\mathrm{,}{\Psi }^{\left(r\right)}\right]$

$\begin{array}{l}Q\left(\Psi \mathrm{<mo>\; ;\; </mo>}{\Psi }^{\left(r\right)}\right)\\ =-{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\tau }_{ij}^{\left(r\right)}\frac{{\delta }_{ij}^{\left(r\right)}}{p{\sigma }_j^2}{\displaystyle \underset{t=1}{\overset{p}{\sum }}}{\left[y_i-{\alpha }_j-x_i^{\text{T}}{\beta }_j^{\left(r\right)}-p{x}_{it}\left({\beta }_{jt}-{\beta }_{jt}^{\left(r\right)}\right)\right]}^2\\ -\frac{n}{2}{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j\left({\displaystyle \underset{t=1}{\overset{p}{\sum }}}\frac{{\lambda }_{jt1}}{{\eta }_{jt}^{\left(r\right)}}{\beta }_{jt}^2+{\displaystyle \underset{k=1}{\overset{K}{\sum }}}\frac{{\lambda }_{jk2}}{{\gamma }_{jk}^{\left(r\right)}}{\displaystyle \underset{t=1}{\overset{p_k}{\sum }}}{\beta }_{jt}^{2\left[k\right]}\right)+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\tau }_{ij}^{\left(r\right)}\left(\log {\pi }_j-\frac{1}{2}\log {\sigma }_j^2\right)\\ -n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j\left[{\displaystyle \underset{t=1}{\overset{p}{\sum }}}{\lambda }_{jt1}\left({\eta }_{jt}^{\left(r\right)}-\frac{{\beta }_{jt}^{2\left(r\right)}}{2{\eta }_{jt}^{\left(r\right)}}\right)+{\displaystyle \underset{k=1}{\overset{K}{\sum }}}{\lambda }_{jk2}\left({\gamma }_{jk}^{\left(r\right)}-\frac{{\displaystyle {\sum }_{t=1}^{p_k}{\beta }_{jt}^{2\left[k\right]\left(r\right)}}}{2{\gamma }_{jk}^{\left(r\right)}}\right)\right]\\ -{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\tau }_{ij}^{\left(r\right)}\left(\frac{1}{2}\log \pi +\log {\delta }_{ij}^{\left(r\right)}+\frac{1}{2{\delta }_{ij}^{\left(r\right)}}\right)\end{array}$ (20)

in a neighborhood of ${\Psi }^{\left(r\right)}$. Note that (20) can be block-wise maximized in the coordinates of the parameter components $\pi$, $\alpha$, $\beta$ and ${\sigma }^2$. Here, $\pi =\left({\pi }_1,\cdots ,{\pi }_g\right)$, $\alpha =\left({\alpha }_1\mathrm{,}\cdots \mathrm{,}{\alpha }_g\right)$, $\beta =\left({\beta }_{11}\mathrm{,}\cdots \mathrm{,}{\beta }_{gp}\right)$ and ${\sigma }^2=\left({\sigma }_1^2\mathrm{,}\cdots \mathrm{,}{\sigma }_g^2\right)$.

Under the constraints that ${\displaystyle {\sum }_{j=1}^g}{\pi }_j=1$ and ${\pi }_j>0$, we adopt Lagrangian multiplier to update $\pi$ by solving

$\nabla Q\left(\pi \mathrm{,}{\alpha }^{\left(r\right)}\mathrm{,}{\beta }^{\left(r\right)}\mathrm{,}{\sigma }^{2\left(r\right)}\mathrm{<mo>\; ;\; </mo>}{\Psi }^{\left(r\right)}\right)=0\mathrm{,}$ (21)

where $\nabla$ is the gradient operator, $\zeta$ is a positive scalar and $0$ is a zero vector. Then we have the set of simultaneous equations

$a_j/{\pi }_j-b_j-\zeta =0,$ (22)

where $a_j={\displaystyle {\sum }_{i=1}^n}{\tau }_{ij}^{\left(r\right)}$ and $b_j=n{\displaystyle {\sum }_{t=1}^p}{\lambda }_{jt1}{\eta }_{jt}^{\left(r\right)}+n{\displaystyle {\sum }_{k=1}^K}{\lambda }_{jk2}{\gamma }_{jk}^{\left(r\right)}$, for each j.

According to ${\pi }_j=a_j/\left(b_j+\zeta \right)$ and ${\displaystyle {\sum }_{j=1}^g}{\pi }_j=1$, we can obtain the unique root ${\zeta }^{\mathrm{<mo>\; *\; </mo>}}$ by solving the equation

${\displaystyle \underset{j=1}{\overset{g}{\sum }}}\frac{a_j}{b_j+\zeta }-1=0.$ (23)

Therefore, the $\left(r+1\right)\text{th}$ iterate

${\pi }_j^{\left(r+1\right)}=\frac{a_j}{b_j+{\zeta }^{*}}.$ (24)

Furthermore, by solving $\nabla Q\left({\pi }^{\left(r+1\right)}\mathrm{,}\alpha \mathrm{,}{\beta }^{\left(r\right)}\mathrm{,}{\sigma }^2\mathrm{<mo>\; ;\; </mo>}{\Psi }^{\left(r\right)}\right)=0$, we have the updates

${\alpha }_j^{\left(r+1\right)}=\frac{{\displaystyle {\sum }_{i=1}^n{\tau }_{ij}^{\left(r\right)}{\delta }_{ij}^{\left(r\right)}\left(y_i-x_i^{\text{T}}{\beta }_j^{\left(r\right)}\right)}}{{\displaystyle {\sum }_{i=1}^n{\tau }_{ij}^{\left(r\right)}{\delta }_{ij}^{\left(r\right)}}}$ (25)

and

${\sigma }_j^{2\left(r+1\right)}=\frac{2{\displaystyle {\sum }_{i=1}^n{\tau }_{ij}^{\left(r\right)}{\delta }_{ij}^{\left(r\right)}{\left(y_i-{\alpha }_j^{\left(r+1\right)}-x_i^{\text{T}}{\beta }_j^{\left(r\right)}\right)}^2}}{{\displaystyle {\sum }_{i=1}^n{\tau }_{ij}^{\left(r\right)}}}.$ (26)

Similarly, for the parameter ${\beta }_{jt}$ in kth group, we obtain the updated formula

${\beta }_{jt}^{\left(r+1\right)}=\frac{{\displaystyle {\sum }_{i=1}^{n}2{\tau }_{ij}^{\left(r\right)}{\delta }_{ij}^{\left(r\right)}{x}_{it}\left(y_i-{\alpha }_j^{\left(r+1\right)}-x_i^{\text{T}}{\beta }_j^{\left(r\right)}+p{x}_{it}{\beta }_{jt}^{\left(r\right)}\right)}}{n{\pi }_j^{\left(r+1\right)}{\sigma }_j^{2\left(r+1\right)}\left({\lambda }_{jt1}/{\eta }_{jt}^{\left(r\right)}+{\lambda }_{jk2}/{\gamma }_{jk}^{\left(r\right)}\right)+2p{\displaystyle {\sum }_{i=1}^n{\tau }_{ij}^{\left(r\right)}{\delta }_{ij}^{\left(r\right)}{x}_{it}^2}}$ (27)

by solving $\nabla Q\left({\pi }^{\left(r+1\right)}\mathrm{,}{\alpha }^{\left(r+1\right)}\mathrm{,}\beta \mathrm{,}{\sigma }^{2\left(r+1\right)}\mathrm{<mo>\; ;\; </mo>}{\Psi }^{\left(r\right)}\right)=0$.

Based on the above, we propose the following EM algorithm.

1) Choose an initial value ${\Psi }^{\left(0\right)}=\left({\pi }_1^{\left(0\right)}\mathrm{,}{\alpha }_1^{\left(0\right)}\mathrm{,}{\beta }_1^{\left(0\right)}\mathrm{,}{\sigma }_1^{2\left(0\right)}\mathrm{,}\cdots \mathrm{,}{\pi }_g^{\left(0\right)}\mathrm{,}{\alpha }_g^{\left(0\right)}\mathrm{,}{\beta }_g^{\left(0\right)}\mathrm{,}{\sigma }_g^{2\left(0\right)}\right)$.

2) E-Step: at the $\left(r+1\right)\text{th}$ iteration, calculate ${\tau }_{ij}^{\left(r\right)}$ and ${\delta }_{ij}^{\left(r\right)}$ by (15) and (16).

3) M-Step: at the $\left(r+1\right)\text{th}$ iteration, update ${\pi }_j$, ${\alpha }_j$, ${\sigma }_j^2$ and ${\beta }_{jt}$ by (24), (25), (26) and (27).

4) Repeat E-Step and M-Step until convergence is obtained.

Note that if a perfect least absolute deviation (LAD) fit occurs in EM algorithm, i.e. $y_i-{\alpha }_j^{\left(r\right)}-x_i^{\text{T}}{\beta }_j^{\left(r\right)}\approx 0$ for some i, j and r. As a result, ${\delta }_{ij}^{\left(r+1\right)}$ will become very large and numerical instability. In this article, we simply introduce a hard threshold to control the extremely small LAD residuals, ${\delta }_{ij}^{\left(r+1\right)}$ will be assigned a value of 10⁶ when the perfect LAD fit occurs.

2.4. Convergence Analysis

The EM algorithm is iterated until some convergence criterion is met. Let tol be a small tolerance constant and M be the maximum iterations for the proposed algorithm, we believe the algorithm has converged to an ideal state when

$\left|F_{n\mathrm{,}\epsilon }\left({\Psi }^{\left(r+1\right)}\right)-F_{n\mathrm{,}\epsilon }\left({\Psi }^{\left(r\right)}\right)\right|<tol\mathrm{,}$ (28)

or the iterations over the maximum iterations M. See [17] for details regarding the relative merits of convergence criteria.

According to Dempster et al. [2], each iteration of the E step and M step of EM algorithm monotonically non-decreases the objective function (8), i.e. $F_{n\mathrm{,}\epsilon }\left({\Psi }^{\left(r+1\right)}\right)-F_{n\mathrm{,}\epsilon }\left({\Psi }^{\left(r\right)}\right)\ge 0$, for all $r\ge 0$. Moreover, Wu [18] proved that if the EM sequence $\left\{{\Psi }^{\left(r\right)}\right\}$ coverges to some point ${\Psi }^{*}$, ${\Psi }^{*}$ is a stationary point of (8) under some general conditions for $F_{n\mathrm{,}\epsilon }\left(\Psi \right)$ and $E\left[F_{n\mathrm{,}\epsilon }\left(\Psi \mathrm{<mo>\; ;\; </mo>}D\right)\mathrm{<mo>\; |\; </mo>}S\mathrm{,}{\Psi }^{\left(r\right)}\right]$. Given the facts above, in this article, we run multiple times from different initializations ${\Psi }^{\left(0\right)}$ in order to obtain an appropriate limit point.

3. Asymptotic Properties

For the regression coefficient vector ${\beta }_j$ in jth component, we can separate it into ${\beta }_j={\left({\beta }_{1j}^{\text{T}}\mathrm{,}{\beta }_{2j}^{\text{T}}\right)}^{\text{T}}$, where ${\beta }_{1j}$ is the set of non-zero effects and ${\beta }_{2j}$ is the set of zero effects. Naturally, we decompose the parameter $\Psi =\left({\Psi }_1\mathrm{,}{\Psi }_2\right)$ such that ${\Psi }_2$ contains all zero effects, namely ${\beta }_{2j}$, $j=1,\cdots ,g$. The true parameter is denoted as ${\Psi }_0$ and the elements of ${\Psi }_0$ are denoted with a subscript, such as ${\beta }_{0jt}$.

For the purpose of easy discussion, we define $a_n=\underset{j\mathrm{,}t}{\max }\left\{{\lambda }_{jt1}\mathrm{<mo>\; :\; </mo>}{\beta }_{0jt}\ne 0\right\}$, $a_n^{\mathrm{<mo>\; *\; </mo>}}=\underset{j\mathrm{,}k}{\max }\left\{{\lambda }_{jk2}\mathrm{<mo>\; :\; </mo>}{\beta }_{0j}^{\left[k\right]}\ne 0\right\}$, $b_n=\underset{j\mathrm{,}t}{\min }\left\{{\lambda }_{jt1}\mathrm{<mo>\; :\; </mo>}{\beta }_{0jt}=0\right\}$, $b_n^{\mathrm{<mo>\; *\; </mo>}}=\underset{j\mathrm{,}k}{\min }\left\{{\lambda }_{jk2}\mathrm{<mo>\; :\; </mo>}{\beta }_{0j}^{\left[k\right]}=0\right\}$. Furthermore, we let $f\left(z\mathrm{<mo>\; ;\; </mo>}\Psi \right)$ be the joint density function of $z=\left(x\mathrm{,}y\right)$ and $\Omega$ be an open parameter space. In order to prove the asymptotic properties of the proposed algorithm, some regularity conditions on the joint distribution of $z$ are also required.

A1. The density $f\left(z\mathrm{<mo>\; ;\; </mo>}\Psi \right)$ has common support in $z$ for all $\Psi \in \Omega$ and $f\left(z\mathrm{<mo>\; ;\; </mo>}\Psi \right)$ is identifiable in $\Psi$ up to a permutation of the components of the mixture.

A2. For each $\Psi \in \Omega$, the density $f\left(z\mathrm{<mo>\; ;\; </mo>}\Psi \right)$ admits third partial derivatives with respect to $\Psi$ for almost all $z$.

A3. For each ${\Psi }_0\in \Omega$, there are functions $M_1\left(z\right)$ and $M_2\left(z\right)$ (possibly depending on ${\Psi }_0$ ) such that for $\Psi$ in a neighborhood of $N\left({\Psi }_0\right)$,

$\left|\frac{\partial f\left(z\mathrm{<mo>\; ;\; </mo>}\Psi \right)}{\partial {\Psi }_u}\right|\le M_1\left(z\right)\mathrm{,}\left|\frac{{\partial }^{2}f\left(z\mathrm{<mo>\; ;\; </mo>}\Psi \right)}{\partial {\Psi }_u\partial {\Psi }_v}\right|\le M_1\left(z\right)\mathrm{,}\left|\frac{{\partial }^3\log f\left(z\mathrm{<mo>\; ;\; </mo>}\Psi \right)}{\partial {\Psi }_u\partial {\Psi }_v\partial {\Psi }_w}\right|\le M_2(z)$

such that ${\displaystyle \int M_1\left(z\right)\text{d}z}<\infty$ and ${\displaystyle \int M_2\left(z\right)f\left(z\mathrm{<mo>\; ;\; </mo>}\Psi \right)\text{d}z}<\infty$.

A4. The Fisher information matrix

$I\left(\Psi \right)=E\left\{\left[\frac{\partial }{\partial \Psi }\log f\left(z\mathrm{<mo>\; ;\; </mo>}\Psi \right)\right]{\left[\frac{\partial }{\partial \Psi }\log f\left(z\mathrm{<mo>\; ;\; </mo>}\Psi \right)\right]}^{\text{T}}\right\}$

is finite and positive definite for each $\Psi \in \Omega$.

Theorem 1. Let $z_i=\left(x_i\mathrm{,}{y}_i\right)$, $i=1,\cdots ,n$, be a random sample from the joint density function $f\left(z\mathrm{<mo>\; ;\; </mo>}\Psi \right)$ that satisfies the regularity conditions A1-A4. Suppose that $\sqrt{n}{a}_n{\to }^{p}0$ and $\sqrt{n}{a}_n^{\mathrm{<mo>\; *\; </mo>}}{\to }^{p}0$, as $n\to 0$, then there is a local maximizer ${\hat{\Psi }}_n$ of the model (5) for which

$\Vert {\hat{\Psi }}_n-{\Psi }_0\Vert =O_p\left\{n^{-1/2}\right\}\mathrm{,}$

where ${\to }^p$ represents convergence in probability.

Proof. Let $r_n=n^{-1/2}$. It suffices that for any given $\epsilon >0$, there is a constant $M_{\epsilon }$ such that

$\underset{n\to \infty }{\lim }P\left\{\underset{\Vert \mu \Vert =M_{\epsilon }}{\sup }{F}_n\left({\Psi }_0+r_n\mu \right)<F_n\left({\Psi }_0\right)\right\}\ge 1-\epsilon \mathrm{.}$ (29)

Now, there is a local maximum in $\left\{{\Psi }_0+r_n\mu \mathrm{<mo>\; :\; </mo>}\Vert \mu \Vert \le M_{\epsilon }\right\}$ with large probability, and this local maximizer ${\hat{\Psi }}_n$ satisfies $\Vert {\hat{\Psi }}_n-{\Psi }_0\Vert =O_p\left(r_n\right)$. Then we let

$\begin{array}{c}{\Delta }_n\left(\mu \right)=F_n\left({\Psi }_0+r_n\mu \right)-F_n\left({\Psi }_0\right)\\ =\left[L_n\left({\Psi }_0+r_n\mu \right)-L_n\left({\Psi }_0\right)\right]-\left[P_n\left({\Psi }_0+r_n\mu \right)-P_n\left({\Psi }_0\right)\right]\mathrm{.}\end{array}$ (30)

Without loss of generality, we assume that the first $d_j$ coefficients of ${\beta }_{0j}$ are non-zero and the first $K_j$ groups contain all non-zero effects of ${\beta }_{0j}$, where ${\beta }_{0j}$ is the true regression coefficient vector in the jth component of the mixture regression model. Hence, we have

$P_n\left({\Psi }_0\right)=n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{t=1}{\overset{d_j}{\sum }}}{\lambda }_{jt1}\left|{\beta }_{0jt}\right|+n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{k=1}{\overset{K_j}{\sum }}}{\lambda }_{jk2}\Vert {\beta }_{0j}^{\left[k\right]}\Vert .$ (31)

Since $P_n\left({\Psi }_0+r_n\mu \right)$ is a sum of non-negative terms, removing terms corresponding to zero effects makes it smaller,

$\begin{array}{c}{\Delta }_n\left(\mu \right)\le L_n\left({\Psi }_0+r_n\mu \right)-L_n\left({\Psi }_0\right)\\ -n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{t=1}{\overset{d_j}{\sum }}}{\lambda }_{jt1}\left[\left|{\beta }_{0jt}+r_n{\mu }_{jt}\right|-\left|{\beta }_{0jt}\right|\right]\\ -n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{k=1}{\overset{K_j}{\sum }}}{\lambda }_{jk2}\left[\Vert {\beta }_{0j}^{\left[k\right]}+r_n{\mu }_j^{\left[k\right]}\Vert -\Vert {\beta }_{0j}^{\left[k\right]}\Vert \right]\mathrm{.}\end{array}$ (32)

By Taylor’s expansion, triangular inequality and arithmetic-geometric mean inequality,

$\begin{array}{l}{L}_n\left({\Psi }_0+r_n\mu \right)-L_n\left({\Psi }_0\right)=n^{-1/2}{L^{\prime }}_n{\left({\Psi }_0\right)}^{\text{T}}\mu -\frac{1}{2}\left[{\mu }^{\text{T}}I\left({\Psi }_0\right)\mu \right]\left\{1+o_p\left(1\right)\right\}\mathrm{,}\\ \left|n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{t=1}{\overset{d_j}{\sum }}}{\lambda }_{jt1}\left[\left|{\beta }_{0jt}+r_n{\mu }_{jt}\right|-\left|{\beta }_{0jt}\right|\right]\right|\le \sqrt{n}{a}_n\sqrt{{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{d}_j}\Vert \mu \Vert \mathrm{,}\\ \left|n{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{k=1}{\overset{K_j}{\sum }}}{\lambda }_{jk2}\left[\Vert {\beta }_{0j}^{\left[k\right]}+r_n{\mu }_j^{\left[k\right]}\Vert -\Vert {\beta }_{0j}^{\left[k\right]}\Vert \right]\right|\le \sqrt{n}{a}_n^{\mathrm{<mo>\; *\; </mo>}}\sqrt{{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{K}_j}\Vert \mu \Vert \mathrm{.}\end{array}$ (33)

Regularity conditions indicate that ${L^{\prime }}_n\left({\Psi }_0\right)=O_p\left(n^{1/2}\right)$ and $I\left({\Psi }_0\right)$ is positive definite, and it is not difficult to find that the sign of ${\Delta }_n\left(\mu \right)$ is completely determined by $-\frac{1}{2}\left[{\mu }^{\text{T}}I\left({\Psi }_0\right)\mu \right]\left\{1+o_p\left(1\right)\right\}$. Therefore, for any given $\epsilon >0$, there is a sufficiently large $M_{\epsilon }$ such that

$\underset{n\to \infty }{\lim }P\left\{\underset{\Vert \mu \Vert =M_{\epsilon }}{\sup }{\Delta }_n\left(\mu \right)<0\right\}>1-\epsilon \mathrm{,}$ (34)

which implies (29), this completes the proof.

Theorem 2. Suppose that the conditions given in Theorem 1 and g is known, $\sqrt{n}{a}_n{\to }^{p}0$, $\sqrt{n}{a}_n^{\mathrm{<mo>\; *\; </mo>}}{\to }^{p}0$, $\sqrt{n}{b}_n{\to }^p\infty$ and $\sqrt{n}{b}_n^{\mathrm{<mo>\; *\; </mo>}}{\to }^p\infty$, as $n\to 0$. Then, for any $\sqrt{n}$ -consistent maximum penalized log-likelihood estimator ${\hat{\Psi }}_n$, we have the following:

1) Sparsity: As $n\to \infty$, $P\left({\hat{\beta }}_{2j}=0\right)\to 0$, $j=1,\cdots ,g$.

2) Asymptotic normality:

$\sqrt{n}\left\{\left[\frac{{P^{″}}_n\left({\Psi }_{01}\right)}{n}+I_1\left({\Psi }_{01}\right)\right]\left({\hat{\Psi }}_1-{\Psi }_{01}\right)+\frac{{P^{\prime }}_n\left({\Psi }_{01}\right)}{n}\right\}{\to }^{d}N\left(0\mathrm{,}{I}_1\left({\Psi }_{01}\right)\right)\mathrm{,}$

where ${\to }^d$ denotes convergence in distribution and $I_1\left({\Psi }_{01}\right)$ is a Fisher information when all zero effects are removed.

Proof. In order to prove the sparsity of Theorem 2, we consider the partition $\Psi =\left({\Psi }_1\mathrm{,}{\Psi }_2\right)$ and let $\left({\hat{\Psi }}_1\mathrm{,}0\right)$ is the maximizer of the penalized log-likelihood function $F_n\left({\Psi }_1\mathrm{,}0\right)$, which is regarded as a function of ${\Psi }_1$. It suffices to show that in the neighborhood $\Vert \Psi -{\Psi }_0\Vert =O\left(n^{-1/2}\right)$, there is the probability $P\left(F_n\left({\Psi }_1\mathrm{,}{\Psi }_2\right)-F_n\left({\hat{\Psi }}_1\mathrm{,}0\right)<0\right)\to 1$ as $n\to \infty$. Then we have

$\begin{array}{c}{F}_n\left({\Psi }_1\mathrm{,}{\Psi }_2\right)-F_n\left({\hat{\Psi }}_1\mathrm{,}0\right)=\left[F_n\left({\Psi }_1\mathrm{,}{\Psi }_2\right)-F_n\left({\Psi }_1\mathrm{,}0\right)\right]+\left[F_n\left({\Psi }_1\mathrm{,}0\right)-F_n\left({\hat{\Psi }}_1\mathrm{,}0\right)\right]\\ \le \left[F_n\left({\Psi }_1\mathrm{,}{\Psi }_2\right)-F_n\left({\Psi }_1\mathrm{,}0\right)\right]\mathrm{.}\end{array}$ (35)

On the other hand,

$F_n\left({\Psi }_1\mathrm{,}{\Psi }_2\right)-F_n\left({\Psi }_1\mathrm{,}0\right)=\left[L_n\left({\Psi }_1\mathrm{,}{\Psi }_2\right)-L_n\left({\Psi }_1\mathrm{,}0\right)\right]-\left[P_n\left({\Psi }_1\mathrm{,}{\Psi }_2\right)-P_n\left({\Psi }_1\mathrm{,}0\right)\right]\mathrm{.}$ (36)

By the mean value theorem,

$L_n\left({\Psi }_1\mathrm{,}{\Psi }_2\right)-L_n\left({\Psi }_1\mathrm{,}0\right)={\left[\frac{\partial L_n\left({\Psi }_1\mathrm{,}\xi \right)}{\partial {\Psi }_2}\right]}^{\text{T}}{\Psi }_2$ (37)

for some $\Vert \xi \Vert \le \Vert {\Psi }_2\Vert =O\left(n^{-1/2}\right)$. By the mean value theorem and regularity condition A3, we can get

$\begin{array}{l}\Vert \frac{\partial L_n\left({\Psi }_1\mathrm{,}\xi \right)}{\partial {\Psi }_2}-\frac{\partial L_n\left({\Psi }_{01}\mathrm{,}0\right)}{\partial {\Psi }_2}\Vert \\ \le \Vert \frac{\partial L_n\left({\Psi }_1\mathrm{,}\xi \right)}{\partial {\Psi }_2}-\frac{\partial L_n\left({\Psi }_1\mathrm{,}0\right)}{\partial {\Psi }_2}\Vert +\Vert \frac{\partial L_n\left({\Psi }_1\mathrm{,}0\right)}{\partial {\Psi }_2}-\frac{\partial L_n\left({\Psi }_{01}\mathrm{,}0\right)}{\partial {\Psi }_2}\Vert \\ \le \left[{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{M}_1\left(z_i\right)\right]\Vert \xi \Vert +\left[{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{M}_1\left(z_i\right)\right]\Vert {\Psi }_1-{\Psi }_{01}\Vert =\left\{\Vert \xi \Vert +\Vert {\Psi }_1-{\Psi }_{01}\Vert \right\}{O}_p\left(n\right)=O_p\left(n^{1/2}\right)\mathrm{.}\end{array}$ (38)

Here ${\Psi }_{01}$ is a subvector of ${\Psi }_0$ with all zero regression coefficients removed. Regularity conditions imply that $\partial L_n\left({\Psi }_{01}\mathrm{,}0\right)/\partial {\Psi }_2=O_p\left(n^{1/2}\right)$, therefore $\partial L_n\left({\Psi }_1\mathrm{,}\xi \right)/\partial {\Psi }_2=O_p\left(n^{1/2}\right)$. In this case, we have

$L_n\left({\Psi }_1\mathrm{,}{\Psi }_2\right)-L_n\left({\Psi }_1\mathrm{,}0\right)=O_p\left(n^{1/2}\right){\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\displaystyle \underset{t=d_j+1}{\overset{p}{\sum }}}\left|{\beta }_{jt}\right|$ (39)

for large n. And for the penalized function $P_n\left(\Psi \right)$,

$\begin{array}{l}{P}_n\left({\Psi }_1\mathrm{,}{\Psi }_2\right)-P_n\left({\Psi }_1\mathrm{,}0\right)\\ \ge {\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{t=d_j+1}{\overset{p}{\sum }}}\sqrt{n}{b}_n\sqrt{n}\left|{\beta }_{jt}\right|+{\displaystyle \underset{j=1}{\overset{g}{\sum }}}{\pi }_j{\displaystyle \underset{k=K_j+1}{\overset{K}{\sum }}}\sqrt{n}{b}_n^{*}\sqrt{n}\Vert {\beta }_j^{\left[k\right]}\Vert \mathrm{.}\end{array}$ (40)

Since $\sqrt{n}{b}_n\to \infty$ and $\sqrt{n}{b}_n^{\mathrm{<mo>\; *\; </mo>}}\to \infty$ as $n\to \infty$, we have

$\left[L_n\left({\Psi }_1\mathrm{,}{\Psi }_2\right)-L_n\left({\Psi }_1\mathrm{,}0\right)\right]-\left[P_n\left({\Psi }_1\mathrm{,}{\Psi }_2\right)-P_n\left({\Psi }_1\mathrm{,}0\right)\right]<0$ (41)

with probability to one as $n\to \infty$. This completes the proof of the sparsity.

For the asymptotic normality of Theorem 2, we still use the same argument as in Theorem 1 and consider $F_n\left({\Psi }_1\mathrm{,}0\right)$ is a function of ${\Psi }_1$, there is a -consistent local maximizer of this function, say, that satisfies

(42)

By the Taylor’s expansion,

(43)

Substituting into (42), we have

(44)

In addition, regularity conditions imply that

(45)

Finally, we can get

(46)

by the Slutsky’s theorem. This completes the proof of the asymptotic normality.

Now, we know that, as long as the conditions, , and are satisfied when, the conclusions of theorem 1 and theorem 2 are tenable. Since the estimator based on the Lasso penalty, it can be -consistent. Then, for any, we have for and for. Based on the fact, we can take which satisfies the and. Similarly, for with for and for, we also take to satisfy the and.

4. Tuning Parameters and Components Selection

In this section, we need to solve two problems. One concerns the number of components g and the other problem is the selection of the tuning parameters. Until now, there is little theoretical support for the selection of these hyper parameters. In former literatures, the cross validation [19] and the generalized cross validation [20] provided some effective guidances for these problems. Grün and Leisch [21] and Nguyen and McLachlan [22] indicated that the Bayesian information criterion (BIC) has a good performance in solving these problems. In this paper, we still use the BIC,

(47)

where the number of non-zero regression coefficients in the jth regression model.

Suppose that there is a set of parameter combinations. For each parameter combination, , we can obtain the parameter estimate by the proposed algorithm, and there is a which depends on corresponding. Finally, we set the best parameter combination for our robust mixture model, where.

5. Numerical Simulation

To quantify the performance of our proposed robust mixture regression model based on adaptive sparse group Lasso (adaSGL-RMR), we design a numerical simulation and generate sample data from the following mixture regression model

(48)

where Z is a component indicator. There are groups and each group consists of 5 covariates, covariates within the same group are correlated, whereas those in different groups are uncorrelated. In order to generate the covariates, we first generate 30 random variables, , independently from, then obtain, , from a multivariate normal distribution with mean zero and. The generation of the covariates as following:

(49)

The model parameters include, , , , , , , ,. The random error and are independent and we consider the follow four cases: 1); 2) Laplace distribution with mean 0 and variance 1; 3), t-distribution with 3 degrees of freedom; 4) a mixture normal distribution.

We use three methods for comparing. The Gaussian mixture model (GMM) based on Lasso penalty (Lasso-GMM), the GMM model based on adaptive Lasso penalty (adaL-GMM) and the adaSGL-RMR model. The fmr package of the R programming language is used to compute the parameter estimates of Lasso-GMM and adaL-GMM.

In this article, the algorithm is terminated when the change in the penalized complete log-likelihood function is less than 10⁻⁵ or meets the maximum iterations 10⁵. Furthermore, we adopt a threshold value for in the consideration of practical situation. To be specific, will be assigned a value of 0 if, for some j and t. To evaluate the performances of variable selection and data fitting, we introduce the average number of selected non-zero variables (nvars) without intercepts, average number of selected non-zero groups (ngroups), frequency of correct identification of group sparsity structures (cgroups), false negative rate (FNR) of missing important predictors, false positive rate (FPR) of selecting unimportant predictors and average value of root mean square errors (RMSE). Here,

(50)

and

(51)

where TP is the number of true positives, TN is the number of true negatives, FP is the number of false positives and FN is the number of false negatives for each fitted model.

As shown in Table 1, in case (1) through case (4), Lasso-GMM fails to identify non-zero groups and selects too many unimportant predictors, adaL-GMM is inclined to select less unimportant predictors and achieves higher cgroups at the cost of ignoring some important predictors. As a contrast, adaSGL-RMR has a better performance in variable selection and RMSE indicates that the parameter estimates of our algorithm are closer to the true parameters of mixture regression model. Moreover, the simulation results clearly show that adaSGL-RMR still maintain its superiority in identifying the group sparsity structures when the distribution of random errors has a heavy tail. Therefore, no matter from the model complexity or the goodness of fit to the data, our proposed method is more competitive than other methods.

Table 1. Results of Lasso-GMM, adaL-GMM and adaSGL-RMR on 100 replicates,.

6. Real Data Analysis

In this section, we will analyze how NBA players’ performances of regular season affect their salaries. We gather salaries for all players in the NBA from the website https://hoopshype.com/salaries/players/2018-2019, during the period from 2018 to 2019. Performance measures for individuals are gathered from the website https://www.foxsports.com/nba/stats?season=2018 in the 2018-2019 regular season, which include scoring, rebounding, assists and defense statistics. By liminating missing data, we obtain a complete dataset, which contains salaries for 248 NBA players and 27 measures of performance.

These performance measures are divided into five groups and covariates in the same group are correlated. Score: points per game (PPG), points per 48 minutes (PTS/48), field goals made per game (FGM/G), field goal attempts per game (FGA/G), 3 point FG made per game (3FGM/G), 3 point FG attempts per game (3FGA/G), free throws made per game (FTM/G), free throw attempts per game (FTA/G), high point total in a single game (HIGH), points per shot (PPS). Rebound: offensive rebounds per game (ORPG), defensive rebounds per game (DRPG), rebounds per 48 minutes (RPG/48), offensive rebound% (OFF REB%), defensive rebound% (DEF REB%), rebound% (REB%). Assist: assists per game (APG), assists per 48 minutes (AST/48), assist% (AST%), turnovers per game (TPG), turnover% (TO%). Steal: steals per game (SPG), steals per 48 minutes (STL/48), steal% (STL%). Block: blocks per game (BPG), blocks per 48 minutes (BLK/48), block% (BLK%).

The matrix should be column standardized to have mean 0 and variance 1 for avoiding a poor fitting result. Then we use the stepAIC function from R package MASS to realize variable selection of the standard linear model via the BIC, the predicted logged salaries from the stepwise-BIC linear model shows a mean square error (MSE) of 0.60 and adjusted R² of 0.42, these terrible results motivate us to conduct further research for this problem. The logged salaries histogram shows multi-modality from Figure 1, it is acceptable to use the mixture regression model for predicting the logged salaries.

For comparison, we run multiple analyses that include three sets of starting parameters for each of models. The predicted results from the adaSGL-RMR model (BIC = 625) have a MSE of 0.11 and adjusted R² of 0.90. The predicted results from the adaSGL-RMR model (BIC = 598) have a MSE of 0.05 and adjusted R² of 0.95. The predicted results from the adaSGL-RMR model (BIC = 517) have a MSE of 0.04 and adjusted R² of 0.96. See Table 2 for more details. These results suggest that the adaSGL-RMR model has the smallest MSE and explains the largest proportion of variance for the logged salaries from the 2018/19 NBA regular season. Moreover, from Figure 2, the predicted densities show a good characterization of the multi-modality in the logged salaries for the adaSGL-RMR models, with the stepwise-BIC linear model not being able to model this.

Figure 1. Histogram and density estimate for logged salaries

Table 2. Parameter estimates for NBA salary data.

Figure 2. Summary of densities for predicted and observed logged salaries.

7. Conclusion

In this paper, we propose a robust mixture regression model based on a Laplace distribution and consider the adaptive sparse group Lasso for variable selection. Its oracle properties are proved completely in Section 3. In addition, the numerical simulation and real data application show that our method has better performance in parameter estimation and variable selection than other methods. A limitation of this study is that we only consider the mixture regression model with K no-overlapping groups and ignore the case when there are some overlaps between different groups. In our future work, we will pay more attention to this problem.

Conflicts of Interest

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

References

[1]	Quandt, R.E. (1972) A New Approach to Estimating Switching Regressions. Journal of the American Statistical Association, 67, 306-310. https://doi.org/10.1080/01621459.1972.10482378
[2]	Dempster, A.P., Laird, N.M. and Rubin, D.B. (1977) Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society: Series B, 39, 1-22. https://doi.org/10.1111/j.2517-6161.1977.tb01600.x
[3]	Jiang, W.X. and Tanner, M.A. (1999) Hierarchical Mixtures-of-Experts for Exponential Family Regression Models: Approximation and Maximum Likelihood Estimation. The Annals of Statistics, 27, 987-1011. https://doi.org/10.1214/aos/1018031265
[4]	McLachlan, G.J. and Peel, D. (2000) Finite Mixture Models. Wiley, New York. https://doi.org/10.1002/0471721182
[5]	Frühwirth-Schnatter, S. (2006) Finite Mixture and Markov Switching Models. Springer, New York.
[6]	Yao, W.X., Wei, Y. and Yu, C. (2014) Robust Mixture Regression Using the t-Distribution. Computational Statistics and Data Analysis, 71, 116-127. https://doi.org/10.1016/j.csda.2013.07.019
[7]	Song, W.X., Yao, W.X. and Xing, Y.R. (2014) Robust Mixture Regression Model Fitting by Laplace Distribution. Computational Statistics and Data Analysis, 71, 128-137. https://doi.org/10.1016/j.csda.2013.06.022
[8]	Wu, Q. and Yao, W.X. (2016) Mixtures of Quantile Regressions. Computational Statistics and Data Analysis, 93, 162-176. https://doi.org/10.1016/j.csda.2014.04.014
[9]	Khalili, A. and Chen, J.H. (2007) Variable Selection in Finite Mixture of Regression Models. Journal of the American Statistical Association, 102, 1025-1038. https://doi.org/10.1198/016214507000000590
[10]	Städler, N., Bühlmann, P. and van de Geer, S. (2010) L1-Penalization for Mixture Regression Models. TEST, 19, 209-256. https://doi.org/10.1007/s11749-010-0197-z
[11]	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
[12]	Fang, K.N., Wang, X.Y., Zhang, S.W., Zhu, J.P. and Ma, S.G. (2015) Bi-Level Variable Selection via Adaptive Sparse Group Lasso. Journal of Statistical Computation and Simulation, 85, 2750-2760. https://doi.org/10.1080/00949655.2014.938241
[13]	Hunter, D.R. and Li, R.Z. (2005) Variable Selection Using MM Algorithms. The Annals of Statistics, 33, 1617-1642. https://doi.org/10.1214/009053605000000200
[14]	Andrews, D.F. and Mallows, C.L. (1974) Scale Mixtures of Normal Distributions. Journal of the Royal Statistical: Series B, 36, 99-102. https://doi.org/10.1111/j.2517-6161.1974.tb00989.x
[15]	Phillips, R.F. (2002) Least Absolute Deviations Estimation via the EM Algorithm. Statistics and Computing, 12, 281-285. https://doi.org/10.1023/A:1020759012226
[16]	Fan, J.Q. and Li, R.Z. (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
[17]	Lange, K. (2013) Optimization. Springer, New York. https://doi.org/10.1007/978-1-4614-5838-8
[18]	Wu, C.F.J. (1983) On the Convergence Properties of the EM Algorithm. The Annals of Statistics, 11, 95-103. https://doi.org/10.1214/aos/1176346060
[19]	Stone, M. (1974) Cross-Validatory Choice and Assessment of Statistical Predictions. Journal of the Royal Statistical Society: Series B, 36, 111-133. https://doi.org/10.1111/j.2517-6161.1974.tb00994.x
[20]	Craven, P. and Wahba, G. (1978) Smoothing Noisy Data with Spline Functions. Numerische Mathematik, 31, 377-403. https://doi.org/10.1007/BF01404567
[21]	Grün, B. and Leisch, F. (2007) Fitting Finite Mixtures of Generalized Linear Regressions in R. Computational Statistics and Data Analysis, 51, 5247-5252. https://doi.org/10.1016/j.csda.2006.08.014
[22]	Nguyen, H.D. and McLachlan, G.J. (2016) Laplace Mixture of Linear Experts. Computational Statistics and Data Analysis, 93, 177-191. https://doi.org/10.1016/j.csda.2014.10.016