Clustering with STARRS

Daphné Giorgi, Antoine Godichon-Baggioni, Stéphane Robin

1 Introduction

Grouping observations into homogeneous groups is a recurrent task in statistical data analysis. Unsupervised classification methods (also called clustering methods) are typically divided into two main classes: distance-based methods, and model-based methods. The most representative examples in each category are, respectively, the \(k\)-means (McQueen (1967)) algorithm and the Gaussian mixture model (McLachlan and Peel (2000)). The STARRS package propose robust counterparts of each of them, namely the \(k\)-medians algorithm (Cardot et al. (2012)) and a robust version of the EM algorithm for Gaussian mixtures (Godichon-Baggioni and Robin (2022)).

2 Model-based clustering: A first example

2.1 Synthetic data set

A sample of \(n=1500\) drawn from a multivariate Gaussian mixture model with \(K=3\) components (with equal proportions) contaminated with outliers with student distribution can be obtained with the function , as follows.

# Dimensions
K <- 3; nk <- rep(500, K); d <- 5
# Gaussian parameters
mu <- rbind(mu1=rep(0, d), mu2=rep(3, d), mu3=rep(-3, d))
sigma <- array(dim=c(K, d, d))
sigma[1, , ] <- 2*diag(d); sigma[2, , ] <- diag(1:d); sigma[3, , ] <- diag(1/(1:d))
# Student contamination
fraction <- 0.3
dfT <- 1; dfT <- 1; muT <- mu; sigmaT <- sigma
# Simulation
GMMcontStudent <- rGMMcontStudent(nk=nk, muG=mu, sigmaG=sigma,
                                  delta=fraction,
                                  dfT=dfT, muT=muT, sigmaT=sigmaT)
print(head(STARRS::GMMcontStudent$X))
#>             [,1]      [,2]       [,3]      [,4]       [,5]
#> [1,]  3.32812044 24.633765 18.2377841 12.764422 -0.2771513
#> [2,]  0.04142344 -1.227728  0.8033976 -1.251689  0.3896963
#> [3,] -3.06277283 -2.454305 -3.7387710 -3.882074 -3.1689136
#> [4,] -2.92069244 -3.220423 -2.2600471 -3.353260 -2.8802425
#> [5,]  1.06043152  1.875461 -0.1896524 -1.978673  3.7184266
#> [6,]  4.41036880  3.274980  5.5915113  4.667790  2.3246754
print(head(STARRS::GMMcontStudent$Z))
#> [1] 1 1 3 3 1 2
print(head(STARRS::GMMcontStudent$C))
#> [1] 1 0 0 0 0 0

2.2 Robust clustering with STARRS

The robust version of the EM algorithm proposed by Godichon-Baggioni and Robin (2024) for the inference of Gaussian mixture models is implemented in the function .

robustGMMOutput <- multipleRobustMM(GMMcontStudent$X, nclust=1:5)
#> Small eigen values are being replaced by  1e-10 
#> var_k has negative eigen values : l= 1 
#> logPhi has infinity values : l= 1
print(robustGMMOutput$bestresult$centers)
#>             [,1]        [,2]        [,3]       [,4]        [,5]
#> [1,] -0.05728465 -0.02208653 -0.09333779 -0.2188448 -0.01100707
#> [2,] -2.99665734 -3.01476676 -3.00741739 -2.9800422 -3.01183761
#> [3,]  3.02756214  2.95296194  2.96684021  3.2089921  3.34233395

The number of components can be chosen with the BIC or ICL criterion.

plot(robustGMMOutput, what='BIC')

plot(robustGMMOutput, what='ICL')

The classification of the observations (in the first two dimensions of PCA) can be visualized with the function of the package, together with confidence of the classification uncertainty.

plot(robustGMMOutput, what='Two_Dim_Uncertainty')

3 Robustness study

3.1 Distance-based methods

The first element of the list of datasets consists of an uncontaminated sample of \(n=1500\) observations drawn from a multivariate Gaussian mixture with \(K=3\) components, with equal proportions and variance matrices all equal to twice the identity. The next datasets from this list consists of the same data set contaminated with an increasing fraction \(f\) of outliers: \(f = 0, 10\%, \dots 50\%\) drawn from a multivariate student mixture.

Both the regular \(k\)-means algorithm and the \(k\)-medians algorithm can be applied to each element of this list.

K <- 3; 
# k-means
kmeans <- lapply(1:length(homoGMMStudentCont_K3_d5_n1500), function(f){
  kmeans(homoGMMStudentCont_K3_d5_n1500[[f]]$X, centers=K, nstart=20, iter.max=100)})
# k-medians
kmedians <- lapply(1:length(homoGMMStudentCont_K3_d5_n1500), function(f){
  kmedians(homoGMMStudentCont_K3_d5_n1500[[f]]$X)})

The next plot illustrate the effect of contamination on the classification accuracy as measure by the adjusted Rand index (ARI).

# Classification accuracy
fractionList <- (0:5)/10; library(mclust)
#> Package 'mclust' version 6.1.1
#> Type 'citation("mclust")' for citing this R package in publications.
plot(fractionList, sapply(1:length(fractionList), function(ff){
  adjustedRandIndex(heteroGMMStudentCont_K3_d5_n1500[[ff]]$Z, kmeans[[ff]]$cluster)
  }), type='b', pch=20, xlab='fraction of outliers', ylab='ARI')
points(fractionList, sapply(1:length(fractionList), function(ff){
  adjustedRandIndex(heteroGMMStudentCont_K3_d5_n1500[[ff]]$Z, kmedians[[ff]]$bestresult$cluster)
  }), type='b', pch=20, col='red')

3.2 Distance-based methods

The first element of the list of datasets has been contructed is the same way as the except that the components of the multivariate Gaussian mixture have different variance matrices.

Again, both the regular EM algorithm (as implemente in the package: Scrucca et al. (2016)) and its robust version can be applied to each element of this list.

# Regular EM
gmm <- lapply(1:length(heteroGMMStudentCont_K3_d5_n1500), function(f){
  Mclust(heteroGMMStudentCont_K3_d5_n1500[[f]]$X, G=K)})
# Robust EM
rgmm <- lapply(1:length(heteroGMMStudentCont_K3_d5_n1500), function(f){
  multipleRobustMM(heteroGMMStudentCont_K3_d5_n1500[[f]]$X, nclust=K)})

Again, the effect of contamination on the classification accuracy (in terms of ARI) can be described as follows.

fractionList <- (0:5)/10
plot(fractionList, sapply(1:length(fractionList), function(ff){
  adjustedRandIndex(heteroGMMStudentCont_K3_d5_n1500[[ff]]$Z, gmm[[ff]]$classification)
  }), type='b', pch=20, xlab='fraction of outliers', ylab='ARI')
points(fractionList, sapply(1:length(fractionList), function(ff){
  adjustedRandIndex(heteroGMMStudentCont_K3_d5_n1500[[ff]]$Z, rgmm[[ff]]$bestresult$clusters)
  }), type='b', pch=20, col='red')

4 References

Cardot, H., P. Cénac, and J-M. Monnez. 2012. “A Fast and Recursive Algorithm for Clustering Large Datasets with k-Medians.” Computational Statistics and Data Analysis 56: 1434–49.
Godichon-Baggioni, A., and S. Robin. 2022. “A Robust Model-Based Clustering Based on the Geometric Median and the Median Covariation Matrix.” Statistics and Computing 34 (1): 55.
Godichon-Baggioni, A., and S. Robin. 2024. “A Robust Model-Based Clustering Based on the Geometric Median and the Median Covariation Matrix.” Statistics and Computing, 24–55.
McLachlan, G., and D. Peel. 2000. Finite Mixture Models. Wiley.
McQueen, J. B. 1967. “Some Methods of Classification and Analysis of Multivariate Observations.” Proc. Of 5th Berkeley Symposium on Math. Stat. And Prob., 281–97.
Scrucca, L., M. Fop, B. T. Murphy, and A. E Raftery. 2016. “mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models.” The R Journal 8 (1): 289–317. https://doi.org/10.32614/RJ-2016-021.