# VarSelLCM

#### Variable Selection for Model-Based Clustering of Mixed-Type Data Set with Missing Values

#### 2025-09-19

References:

- Marbac, M. and Sedki, M. (2017), Variable selection for model-based clustering using the integrated complete-data likelihood, Statistics and Computing, Volume 27, Issue 4, pp 1049–1063.
- Marbac, M., Patin, E. and Sedki, M. (2018), Variable selection for mixed data clustering: Application in human population genomics, Journal of Classification, to appear.

# Introduction

VarSelLCM permits a full model selection (detection of the relevant features for clustering and selection of the number of clusters) in model-based clustering, according to classical information criteria (BIC, MICL or AIC).

Data to analyzed can be composed of continuous, integer and/or categorical features. Moreover, missing values are managed, without any pre-processing, by the model used to cluster with the assumption that values are missing completely at random. Thus, VarSelLCM can also be used for data imputation via mixture models.

An R-Shiny application is implemented to easily interpret the clustering results

# Mixed-type data analysis

## Clustering

This section performs the whole analysis of the Heart data set. Warning the univariate margin distribution are defined by class of the features: numeric columns imply Gaussian distributions, integer columns imply Poisson distribution while factor (or ordered) columns imply multinomial distribution

    library(VarSelLCM)

    Attachement du package : 'VarSelLCM'

    L'objet suivant est masqué depuis 'package:stats':

        predict

    # Data loading:
    # x contains the observed variables
    # z the known status (i.e. 1: absence and 2: presence of heart disease)
    data(heart)
    ztrue <- heart[,"Class"]
    x <- heart[,-13]
    # Add a missing value artificially (just to show that it works!)
    x[1,1] <- NA

Clustering is performed with variable selection. Model selection is done with BIC because the number of observations is large (compared to the number of features). The number of components is between 1 and 3. Do not hesitate to use parallelization (here only two cores are used).

    # Cluster analysis without variable selection
    res_without <- VarSelCluster(x, gvals = 1:3, vbleSelec = FALSE, crit.varsel = "BIC")

    # Cluster analysis with variable selection (with parallelisation)
    res_with <- VarSelCluster(x, gvals = 1:3, nbcores = 2, crit.varsel = "BIC")

Comparison of the BIC for both models: variable selection permits to improve the BIC

    BIC(res_without)

    [1] -6516.216

    BIC(res_with)

    [1] -6509.506

Comparison of the partition accuracy. ARI is computed between the true partition (ztrue) and its estimators. ARI is an index between 0 (partitions are independent) and 1 (partitions are equals). Variable selection permits to improve the ARI. Note that ARI cannot be used for model selection in clustering, because there is no true partition.

    ARI(ztrue, fitted(res_without))

    [1] 0.2218655

    ARI(ztrue, fitted(res_with))

    [1] 0.2661321

To obtained the partition and the probabilities of classification

    # Estimated partition
    fitted(res_with)

      [1] 1 1 2 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 2 1 2 2 2 2 2 1 2 1 1 1 1 2 1 1 1 1
     [38] 1 2 2 2 2 2 2 1 2 1 2 1 1 1 2 2 2 2 2 1 1 1 1 2 1 2 2 1 1 2 2 2 2 1 2 2 1
     [75] 1 1 1 2 2 2 1 1 1 2 1 2 2 1 2 1 2 2 1 1 2 1 1 1 1 2 2 1 2 1 1 1 1 1 1 2 1
    [112] 2 2 2 2 2 2 1 1 1 1 1 1 2 2 2 1 2 2 1 1 1 2 1 2 2 2 1 2 1 1 2 1 1 2 1 2 1
    [149] 2 2 2 2 2 1 2 2 1 2 1 1 1 1 2 1 2 1 2 2 1 1 1 1 1 1 2 1 1 2 1 2 2 1 2 1 2
    [186] 2 1 1 2 1 2 1 2 2 2 2 1 2 1 1 1 1 1 1 1 2 2 1 1 2 1 2 2 1 2 2 2 1 1 2 1 1
    [223] 2 1 2 1 1 1 2 1 1 1 2 1 1 1 2 1 2 2 1 2 1 1 2 1 1 2 2 1 1 2 2 1 2 2 1 1 2
    [260] 2 2 1 2 2 2 2 2 2 1 1

    # Estimated probabilities of classification
    head(fitted(res_with, type="probability"))

           class-1      class-2
    [1,] 0.9999917 8.261532e-06
    [2,] 0.6334688 3.665312e-01
    [3,] 0.1755334 8.244666e-01
    [4,] 1.0000000 4.443137e-08
    [5,] 0.9961153 3.884740e-03
    [6,] 0.9547833 4.521665e-02

To get a summary of the selected model.

    # Summary of the best model
    summary(res_with)

    Model:
       Number of components: 2
       Model selection has been performed according to the BIC  criterion
       Variable selection has been performed, 8  ( 66.67 % ) of the variables are relevant for clustering

Discriminative power of the variables (here, the most discriminative variable is MaxHeartRate). The greater this index, the more the variable distinguishes the clusters.

    plot(res_with)

Distribution of the most discriminative variable per clusters

    # Boxplot for the continuous variable MaxHeartRate
    plot(x=res_with, y="MaxHeartRate")

Empirical and theoretical distributions of the most discriminative variable (to check that the distribution is well-fitted)

    # Empirical and theoretical distributions (to check that the distribution is well-fitted)
    plot(res_with, y="MaxHeartRate", type="cdf")

    Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
    ℹ Please use `linewidth` instead.
    ℹ The deprecated feature was likely used in the VarSelLCM package.
      Please report the issue to the authors.
    This warning is displayed once every 8 hours.
    Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
    generated.

Distribution of a categorical variable per clusters

    # Summary of categorical variable
    plot(res_with, y="Sex")

To have details about the selected model

    # More detailed output
    print(res_with)

    Data set:
       Number of individuals: 270
       Number of continuous variables: 3
       Number of count variables: 1
       Percentile of missing values for the integer variables: 0.37
       Number of categorical variables: 8

    Model:
       Number of components: 2
       Model selection has been performed according to the BIC  criterion
       Variable selection has been performed, 8  ( 66.67 % ) of the variables are relevant for clustering

    Information Criteria:
       loglike: -6403.136
       AIC:     -6441.136
       BIC:     -6509.506
       ICL:     -6638.116

To print the parameters

    # Print model parameter
    coef(res_with)

    An object of class "VSLCMparam"
    Slot "pi":
      class-1   class-2
    0.5221134 0.4778866

    Slot "paramContinuous":
    An object of class "VSLCMparamContinuous"
    Slot "pi":
    numeric(0)

    Slot "mu":
                       class-1  class-2
    RestBloodPressure 131.3444 131.3444
    SerumCholestoral  249.6593 249.6593
    MaxHeartRate      135.4167 165.2587

    Slot "sd":
                       class-1  class-2
    RestBloodPressure 17.82850 17.82850
    SerumCholestoral  51.59043 51.59043
    MaxHeartRate      20.98139 13.14845

    Slot "paramInteger":
    An object of class "VSLCMparamInteger"
    Slot "pi":
    numeric(0)

    Slot "lambda":
         class-1 class-2
    Age 58.11337 50.3206

    Slot "paramCategorical":
    An object of class "VSLCMparamCategorical"
    Slot "pi":
    numeric(0)

    Slot "alpha":
    $Sex
                    0         1
    class-1 0.2358079 0.7641921
    class-2 0.4166339 0.5833661

    $ChestPainType
                     1          2         3         4
    class-1 0.08922364 0.03291641 0.1738646 0.7039954
    class-2 0.05752247 0.28954453 0.4223085 0.2306245

    $FastingBloodSugar
                    0         1
    class-1 0.8518519 0.1481481
    class-2 0.8518519 0.1481481

    $ResElectrocardiographic
                    0           1         2
    class-1 0.4851852 0.007407407 0.5074074
    class-2 0.4851852 0.007407407 0.5074074

    $ExerciseInduced
                    0          1
    class-1 0.4484672 0.55153282
    class-2 0.9128099 0.08719007

    $Slope
                    1         2          3
    class-1 0.2266442 0.6884266 0.08492919
    class-2 0.7599031 0.1933829 0.04671401

    $MajorVessels
                    0         1          2           3
    class-1 0.4104432 0.2830467 0.17928534 0.127224825
    class-2 0.7915993 0.1402683 0.05987811 0.008254233

    $Thal
                    3            6         7
    class-1 0.3183111 9.931148e-02 0.5823774
    class-2 0.8302565 4.578551e-12 0.1697435

Probabilities of classification for new observations

    # Probabilities of classification for new observations
    predict(res_with, newdata = x[1:3,])

           class-1      class-2
    [1,] 0.9999914 8.635608e-06
    [2,] 0.6231271 3.768729e-01
    [3,] 0.1692163 8.307837e-01

The model can be used for imputation (of the clustered data or of a new observation)

    # Imputation by posterior mean for the first observation
    not.imputed <- x[1,]
    imputed <- VarSelImputation(res_with, x[1,], method = "sampling")
    rbind(not.imputed, imputed)

      Age Sex ChestPainType RestBloodPressure SerumCholestoral FastingBloodSugar
    1  NA   1             4               130              322                 0
    2  74   1             4               130              322                 0
      ResElectrocardiographic MaxHeartRate ExerciseInduced Slope MajorVessels Thal
    1                       2          109               0     2            3    3
    2                       2          109               0     2            3    3

## Shiny application

All the results can be analyzed by the Shiny application…

    # Start the shiny application
    VarSelShiny(res_with)

… but this analysis can also be done on R.