Robust and sparse multigroup classification by the optimal scoring approach

Abstract

We propose a robust and sparse classification method based on the optimal scoring approach. It is also applicable if the number of variables exceeds the number of observations. The data are first projected into a low dimensional subspace according to an optimal scoring criterion. The projection only includes a subset of the original variables (sparse modeling) and is not distorted by outliers (robust modeling). In this low dimensional subspace classification is performed by minimizing a robust Mahalanobis distance to the group centers. The low dimensional representation of the data is also useful for visualization purposes. We discuss the algorithm for the proposed method in detail. A simulation study illustrates the properties of robust and sparse classification by optimal scoring compared to the non-robust and/or non-sparse alternative methods. Three real data applications are given.

Appendix

1.1 Derivation of expression (4) for the score vector estimates

Let $\omega _1,\ldots ,\omega _n$ be case weights for each observation. $\varvec{\varOmega }$ is a diagonal matrix with these case weights in the diagonal. Then the weighted data matrices are $\varvec{\tilde{Y}}=\varvec{\varOmega }^{1/2}\varvec{Y}$ and $\varvec{\tilde{X}}=\varvec{\varOmega }^{1/2}\varvec{X}$. The diagonal matrix with weighted class proportions is $\varvec{\tilde{D}}=\frac{1}{\sum \omega _i}\varvec{\tilde{Y}}^T \varvec{\tilde{Y}}$. The optimization problem (2) in step h for a given $\varvec{\hat{\beta }}$ can be rewritten as

$$\begin{aligned} \min _{\varvec{\theta }}\Vert \varvec{\tilde{X}}\varvec{\hat{\beta }} - \varvec{\tilde{Y}}\varvec{\theta } \Vert ^2 \quad \text {s.t.} \quad \varvec{\theta }^T\varvec{\tilde{D}}\varvec{\theta }=1 \text { and } \varvec{C}\varvec{\theta }=\varvec{0} \in \mathbb {R}^h \end{aligned}$$

(6)

with $\varvec{C}=[\hat{\varvec{\theta }}_{1},\ldots , \hat{\varvec{\theta }}_{h-1}]^T\varvec{\tilde{D}}$, and we drop the depending on the index h for ease of notation.

We use the method of Lagrange multipliers. The Lagrangian associated to Eq. (6) is given by

$$\begin{aligned} L=(\varvec{\tilde{X}}\varvec{\hat{\beta }} - \varvec{\tilde{Y}}\varvec{\theta })^T (\varvec{\tilde{X}}\varvec{\hat{\beta }} - \varvec{\tilde{Y}}\varvec{\theta }) - \eta (\varvec{\theta }^T\varvec{\tilde{D}}\varvec{\theta } - 1) - 2\varvec{\gamma }^T\varvec{C}\varvec{\theta }. \end{aligned}$$

The partial derivative set to zero gives

$$\begin{aligned} \frac{\partial L}{\partial \varvec{\theta }}=-2\varvec{\tilde{Y}}^T(\varvec{\tilde{X}} \varvec{\hat{\beta }}-\varvec{\tilde{Y}}\varvec{\theta }) -2\eta \varvec{\tilde{D}}\varvec{\theta } -2\varvec{C}^T\varvec{\gamma }=\varvec{0}. \end{aligned}$$

Hence,

$$\begin{aligned} \varvec{\theta }=(\varvec{\tilde{Y}}^T\varvec{\tilde{Y}} -\eta \varvec{\tilde{D}})^{-1}(\varvec{\tilde{Y}}^T \varvec{\tilde{X}}\varvec{\hat{\beta }} +\varvec{C}^T\varvec{\gamma }). \end{aligned}$$

To solve for the Lagrange multipliers $\eta$ and $\varvec{\gamma }$, the side constraints are used.

$$\begin{aligned} 0=\varvec{C}\varvec{\theta }= \varvec{C}(\varvec{\tilde{Y}}^T\varvec{\tilde{Y}}-\eta \varvec{\tilde{D}})^{-1} \varvec{\tilde{Y}}^T\varvec{\tilde{X}}\varvec{\hat{\beta }}+ \varvec{C}(\varvec{\tilde{Y}}^T\varvec{\tilde{Y}}-\eta \varvec{\tilde{D}})^{-1} \varvec{C}^T\varvec{\gamma } \end{aligned}$$

So

$$\begin{aligned} \varvec{\gamma }=-\left( \varvec{C} (\varvec{\tilde{Y}}^T \varvec{\tilde{Y}}-\eta \varvec{\tilde{D}})^{-1} \varvec{C}^T \right) ^{-1} \varvec{C} (\varvec{\tilde{Y}}^T\varvec{\tilde{Y}}- \eta \varvec{\tilde{D}})^{-1} \varvec{\tilde{Y}}^T \varvec{\tilde{X}}\varvec{\hat{\beta }}. \end{aligned}$$

We conclude

$$\begin{aligned} \varvec{\theta }= & {} (\varvec{\tilde{Y}}^T \varvec{\tilde{Y}} -\eta \varvec{\tilde{D}})^{-1}\nonumber \\&\left\{ \varvec{I} - \varvec{C}^T(\varvec{C}(\varvec{\tilde{Y}}^T \varvec{\tilde{Y}}- \eta \varvec{\tilde{D}})^{-1}\varvec{C}^T)^{-1} \varvec{C} (\varvec{\tilde{Y}}^T\varvec{\tilde{Y}} - \eta \varvec{\tilde{D}})^{-1} \right\} (\varvec{\tilde{Y}}^T \varvec{\tilde{X}}\varvec{\hat{\beta }}). \end{aligned}$$

(7)

Since $\varvec{\tilde{Y}}^T\varvec{\tilde{Y}}$ is proportional to $\varvec{\tilde{D}}$, there exists a scalar c such that

$$\begin{aligned} (\varvec{\tilde{Y}}^T\varvec{\tilde{Y}} - \eta \varvec{\tilde{D}})^{-1}=c\varvec{\tilde{D}}^{-1}. \end{aligned}$$

Formula (7) can be simplified to

$$\begin{aligned} \varvec{\theta }=c\left\{ \varvec{I} - \varvec{\tilde{D}}^{-1} \varvec{C}^T (\varvec{C} \varvec{\tilde{D}}^{-1} \varvec{C}^T)^{-1} \varvec{C}\right\} \varvec{\tilde{D}}^{-1} \varvec{\tilde{Y}}^T \varvec{\tilde{X}} \varvec{\hat{\beta }}. \end{aligned}$$

Due to the symmetry of $\varvec{\tilde{D}}$ and with the definition of $\varvec{C}=\varvec{Q}^T\varvec{\tilde{D}}$ we obtain

$$\begin{aligned} \varvec{\theta }=c\left\{ \varvec{I} - \varvec{Q} (\varvec{Q}^T \varvec{\tilde{D}} \varvec{Q})^{-1} \varvec{Q}^T \varvec{\tilde{D}} \right\} \varvec{\tilde{D}}^{-1} \varvec{\tilde{Y}}^T \varvec{\tilde{X}} \varvec{\hat{\beta }}. \end{aligned}$$

The scalar c can then be scaled so that the side constraint $\varvec{\theta }^T\varvec{\tilde{D}}\varvec{\theta }=1$ is fulfilled.

1.2 Algorithm for the computation of the initial estimates for $\varvec{\beta }_h$ and $\varvec{\theta }_h$

Input: $h, \varvec{Q}_h, \varvec{X}, \varvec{Y}$, $\lambda$

• (i) Compute $\varvec{D}=\frac{1}{n}\varvec{Y}^T\varvec{Y}$.

• (ii) Generate $\varvec{\theta }_{*}$, a random vector from N(0, 1) of length K.

• (iii) Compute $\hat{\varvec{\theta }}_h=c\left\{ \varvec{I} - \varvec{Q}_h (\varvec{Q}_h^T \varvec{D} \varvec{Q}_h)^{-1} \varvec{Q}_h^T \varvec{D} \right\} \varvec{\theta }_{*}$, with c so that $\hat{\varvec{\theta }}_h^T\varvec{D}\hat{\varvec{\theta }}_h=1$.

Apply twice the following steps:

• 1. For fixed $\hat{\varvec{\theta }}_h$ apply sparse least trimmed squares (sparse LTS) regression (Alfons et al. 2013) to the response $\varvec{Y}\hat{\varvec{\theta }}_h$ and predictors $\varvec{X}$.

  Let $a=0.5n$ and $\Vert \varvec{r}\Vert ^2_{1:a}=\sum _{i=1}^a r_{(i)}^2$ denote the sum of the a smallest squared elements of the vector $\varvec{r}$. The sparse LTS estimator is a robust version of the Lasso and defined as

  $$\begin{aligned} \min _{\varvec{\beta }}\frac{1}{a}\Vert \varvec{Y} \varvec{\theta }_h-\varvec{X}\varvec{\beta }\Vert ^2_{(1):(a)} + \lambda \Vert \varvec{\beta }\Vert _1. \end{aligned}$$

  As in Alfons et al. (2013), a re-weighting step is carried out afterwards yielding $\hat{\varvec{\beta }}_h$.

• 2. For fixed $\hat{\varvec{\beta }}_h$ apply least absolute deviation (LAD) regression with response $\varvec{X}\hat{\varvec{\beta }}_h$ and predictor matrix $\varvec{Y}$:

  $$\begin{aligned} \varvec{\theta }^*=\mathop {\text{ argmin }}_{\varvec{\theta }}\Vert \varvec{Y}\varvec{\theta }-\varvec{X} \hat{\varvec{\beta }}_h\Vert _1. \end{aligned}$$

  The LDA estimator is robust to outliers in the dependent variable, but not to leverage points (i.e. outliers in the covariate space). Since the covariates are dummy variables here, leverage points cannot occur. Then we apply the transformation for satisfying the side constraints:

  $$\begin{aligned} \hat{\varvec{\theta }}_h=c\left\{ \varvec{I} - \varvec{Q}_h (\varvec{Q}_h^T \varvec{D} \varvec{Q}_h)^{-1} \varvec{Q}_h^T \varvec{D} \right\} \varvec{\theta }^*. \end{aligned}$$

Output: Initial estimators $\hat{\varvec{\beta }}_h$ and $\hat{\varvec{\theta }}_h$.