title: reval: a Python package to determine best clustering solutions with stability-based relative clustering validation
authors: Isotta Landi Veronica Mandelli Michael V. Lombardo
text: 
  H-space and loop space structures for intermediate curvatures
      David J. Wraith
      Received December 30, 2023; accepted December 30, 2023
  =============================================================
  
  
  
  
  Determining the best partition for a dataset can be a challenging task because of 1) the lack of a priori information within an unsupervised learning framework; and 2) the absence of a unique clustering validation approach to evaluate clustering solutions. Here we present reval: a Python package that leverages stability-based relative clustering validation methods to determine best clustering solutions as the ones that best generalize to unseen data. 
      
  Statistical software, both in R and Python, usually rely on internal validation metrics, such as silhouette, to select the number of clusters that best fits the data. Meanwhile, open-source software solutions that easily implement relative clustering techniques are lacking. Internal validation methods exploit characteristics of the data itself to produce a result, whereas relative approaches attempt to leverage the unknown underlying distribution of data points looking for generalizable and replicable results.
  
  The implementation of relative validation methods can further the theory of clustering by enriching the already available methods that can be used to investigate clustering results in different situations and for different data distributions. This work aims at contributing to this effort by developing a stability-based method that selects the best clustering solution as the one that replicates, via supervised learning, on unseen subsets of data. The package works with multiple clustering and classification algorithms, hence allowing both the automatization of the labeling process and the assessment of the stability of different clustering mechanisms.
  
  
  
  
  
  § INTRODUCTION
   
  
  
  Clustering algorithms identify intrinsic subgroups in a dataset by arranging together elements that show smaller pairwise dissimilarities relative to other subgroups <cit.>. They are one of a number of machine learning methods that do unsupervised learning with the aim of identifying patterns in the data in the absence of supervised/external knowledge. While their usage is relatively widespread, the lack of a priori information complicates the evaluation of clustering solutions. Attempts to address this challenge usually rely on internal validation measures, that focus on quantities and features inherent to a grouping solution <cit.>. Many tools are available to compute internal validation measures that help in determining the best number of clusters. For example, the elbow method <cit.> selects the number of clusters for which the within-cluster variability decrement is minimal. Another popular method using internal criteria is the silhouette-based approach <cit.>. This method maximizes cluster cohesion and separation - that is, how similar an object is to other elements of the same cluster compared to elements of other clusters. 
  
  In contrast to internal validation, other approaches, such as relative validation methods, have the potential to transform cluster analysis into a model selection problem and help evaluate the best clustering solution (i.e., best number of clusters). The way these methods are conceived also offer up the possibility to determine the extent to which a clustering solution generalizes to unseen data and hence to enable the replication of the data partition chosen. While a variety of software packages contain internal cluster validation methods and measures, open-source software to easily implement the full potential of relative validation techniques are lacking. Here we present the reval Python package (pronounced “reh-val”, like the word “revel”) to help fill this important missing gap. reval implements an approach for stability-based validation of clustering solutions described by Lange and colleagues <cit.> which allows for the identification and evaluation of partitions that best generalize to unseen data and the automation of the labeling process.
  
  Stability-based methods return the number of clusters that minimizes the expected distance between clustering solutions obtained for different datasets. Several options are available <cit.> to 1) generate the datasets (e.g., random subsampling of the original dataset <cit.>, or adding random noise <cit.>); 2) compare clustering solutions (e.g., overlapping subsamples <cit.>); and 3) compute clustering distances (e.g., the consensus index by Vinh and Epps <cit.>). The method proposed by Lange et al. <cit.> has the advantage of transforming the unsupervised setting into a classification problem and guides selection through the minimization of prediction error. First, a dataset is split into training and validation sets and then independently partitioned into clusters. Second, training set labels are used within supervised classification methods to learn how to best predict the labels. Applying the classification model to the validation set, the model’s predicted labels are then compared to the actual clustering labels derived from the validation set. This procedure is repeated using cross-validation and the optimal number of clusters is identified corresponding to the maximum number of clusters that minimizes prediction error. Prediction performance is defined by the authors as the 0-1 loss in supervised classification <cit.>, namely, the normalized Hamming distance. Nevertheless, other choices are possible, for example, Tibshirani and colleagues used prediction strength - that is, the proportion of observation pairs in the validation set that are assigned to the same cluster by both the clustering algorithm and the partition based on the training set centers <cit.>. 
  
  Libraries and methods for the automated selection of the best number of clusters are available in both Python and R. The yellowbrick Python visual analysis and diagnostic tool suite <cit.> includes the implementation of the elbow method to determine the best number of clusters. In R, NbClust <cit.> is a popular library that compiles 30 different internal metrics and allows for users to compute all or a subset of these metrics for use in a majority vote rule to select the optimal number of clusters. For relative validation approaches there are the clValid <cit.> and cstab <cit.> libraries that apply stability-based relative validation models. clValid was designed to work with highly correlated high-throughput genomic data and computes stability measures comparing clustering solutions based on full data and data with a single column removed. cstab implements the selection of the best number of clusters via model-based and model-free clustering instability <cit.> using a bootstrap approach.
  
  The reval package contributes to this landscape by implementing a stability-based approach that can be easily applied to different datasets using multiple clustering and classification algorithms. Built on top of the stability-based algorithm <cit.>, reval applies a classifier trained on the best clustering solution to a test set, returning classification metrics that help interpret the generalization performance, guide the clustering process, and enable labeling replication. Such a tool can be used in concert with internal measures to assess the underlying structure of a dataset to help avoid the risk of overfitting. With respect to clustering errors internal and relative indices can exhibit similar behavior with the advantage of the former being less computationally expensive <cit.>. However, in the case of complex models and clusters, an approach based on the minimization of prediction error may be particularly advantageous because internal indices tend to fail to correlate well with errors <cit.>.
  
  
  
  
  § THEORETICAL DETAILS
   
  
  
  
  
   §.§ Stability measure
  
  
  The notion of stability by Lange et al. <cit.> is used to assess solutions of clustering algorithms based on the rationale that true clusters are those that can always be identified by a clustering algorithm when applied to different datasets from the same generating process. Formally, let 𝒜_k be a clustering algorithm with k the number of clusters, ϕ a classifier, and (𝐗, 𝐘) the training set and clustering labels, i.e., 𝒜_k(𝐗)=𝐘. After training ϕ on (𝐗, 𝐘), both the clustering algorithm and trained classifier are applied to a separate dataset 𝐗'. The distance between the two solutions is the normalized Hamming distance:
  
  
      d_S_k(ϕ(𝐗'),𝐘') = min_σ∈ S_k1/n∑_i=1^n1_{σ(ϕ(X_i'))≠ Y_i'}
  
  
  with S_k the set of all possible permutations of k elements. Supervised labels are permuted to overcome the non-uniqueness of clustering labeling and σ is the permutation that minimizes the solutions dissimilarity. Averaging out the distance between any pair of partitions 𝐗, 𝐗' from Eq. <ref> we define the stability index of the clustering algorithm as:
  
  
      𝒮(𝒜_k) = _𝐗,𝐗'[d_S_k(ϕ(𝐗'), 𝐘')]
    
  The stability index ranges in [0, 1], with lower values indicating more stable and reproducible solutions <cit.>. Because this measure scales with the number of clusters, the measure suggested by the authors is the normalized stability 𝒮̅_k, i.e., the stability from Eq. <ref> normalized for the stability of random labeling ℛ_k.
  
  
  
   §.§ reval algorithm
  
  The algorithm implemented in reval allows the user to: 1) automatically select the number of clusters for a dataset by minimizing validation stability, via repeated cross validation (see Algorithm <ref>); and 2) compute classification performance obtained when generalizing the solution to an held-out dataset (see Algorithm <ref>). An overview of the framework is reported in Figure <ref>.
  
  A dataset 𝐗 is first split  into training 𝐗_tr and test 𝐗_ts sets and a clustering 𝒜 and classifier ϕ are selected. Let n_fold be the number of folds for cross validation, n_r the number of repetitions, n_rnd the number of random labeling iterations, and k the number of clusters in set K. In Algorithm <ref> we indicate with 𝐗^i_j_itr and 𝐗^i_j_val the internal training and validation splits of training set 𝐗_tr, respectively, for cross-validation ith fold split at the jth shuffled repetition. These correspond to 𝐗 and 𝐗' sets introduced in the “Stability measure” section. With (K× n_fold× n_r) we indicate the Cartesian product of the sets of number of clusters and repeated cross-validation splits. The fitted model becomes the one trained on 𝐗_itr that returns the maximum number of clusters with minimum stability. That model can then be used within Algorithm <ref> for generalization on the test set.
  
  Among clustering methods that work within reval, density-based clustering HDBSCAN <cit.> does not need any assumption on the number of clusters. Hence, we do not need to iterate over different number of clusters to select the best solution. Instead, normalized stability is computed within the repeated cross-validation loops that return the same number of clusters.
  
  
  
  
  
  
  
  
  
  
  § TECHNICAL DETAILS
  
  
  
   §.§ Implementation
  
  The reval package has four core modules:
  
  
  
    * : This module includes training and test methods that return misclassification errors obtained by comparing classification labels and clustering labels. It also includes the random labeling method which allows users to compute the asymptotic misclassification rate.
  
  
    * : This module implements repeated cross validation and returns the best clustering solution together with normalized stability scores, obtained from the average of the misclassification scores divided by the asymptotic misclassification rate. Repeated cross validation leads to unbiased stability estimates and it can also be performed stratifying the repeated randomized splits according to a desired variable. To control for the size imbalance that derives from cross validation we initialized the repeated cross-validation loop to a 1× 2 schema as default. Users can change this configuration according to dataset size and available stratifiers that can be useful to overcome imbalance issues. The evaluation method applies the fitted model with the returned number of clusters to the held-out dataset and returns accuracy (ACC). Other metrics, such as Matthews correlation coefficient (MCC) <cit.>, F1 score, precision and recall scores can also be computed (see utils file).
  
  
    * : This module enables hyperparameter tuning to select the best configuration of classifier/clustering ( class) and the parameters within clustering and classifier themselves ( class). Best parameters are those that report minimum normalized stability. If the number of true classes is available this module also returns the best solution that correctly identifies the true number of clusters, if it exists.
  
  
    * : This module includes the function to plot cross-validation performance metrics with 95% confidence intervals for a varying number of clustering solutions. The threshold of random labeling stability can be displayed to visually investigate model performance.
  
  
  As suggested by Lange and colleagues <cit.>, we used the Kuhn-Munkres algorithm <cit.> to obtain the label permutation that minimizes misclassification error. However, differently from the work of Lange and colleagues <cit.>, reval permutes the clustering labels instead of the classification labels, i.e., the normalized Hamming distance in Eq. <ref> becomes:
  
  
      d_S_k(ϕ(𝐗'),𝐘') = min_σ∈S_k1/n∑_i=1^n1_{ϕ(X_i')≠σ(Y_i')}
  
  
  This approach allows the test set to preserve the partition structure of the training set, to better investigate results replicability and to aid visual comparison. Figure <ref> shows the rationale behind the need to permute clustering labels when training a classifier within reval. We simulated three Gaussian blobs and divided them into training (N=20) and test (N=10) sets. Figure <ref> shows clustering labels for the training set and Figure <ref> the clustering labels for the test set. Figure <ref> shows what might happen when training a classifier on the training set labels and then predict labels for the test set. Because two out of three classes show label discordance the trained classifier fails to correctly predict the classes. Nevertheless, if we permute class labels 0 and 1 in the example, the trained classifier will correctly identified all three classes on the test set returning minimum prediction error and consistent label ordering is preserved. The kuhn_munkres_algorithm function can be found in utils file.
  
  A more thorough description of the code and its usage can be found at <https://reval.readthedocs.io/en/latest/code_description.html>. reval mainly works with the scikit-learn Python library for machine learning <cit.>. In particular, among clustering methods, users can select those with number of clusters parameter, i.e., k-means, hierarchical clustering, and spectral clustering, but also density-based clustering HDBSCAN <cit.> from hdbscan library. Moreover, any classifier from scikit-learn can be selected.
  
  
  
   §.§ Algorithm complexity
  
  We report here the complexity analysis of the two core methods included in  module. In particular, we focus on best_nclustcv, which enables the selection of the best number of clusters via repeated cross validation, and evaluate, which implements the testing of the best solution on held-out datasets. The best_nclustcv method includes sequential calls to train (train), test (test), and random labeling (rndlabels_traineval) methods from the  module, whereas evaluate sequentially calls the train and test functions.
  
  The overall complexity to perform cross validation and evaluation depends on the number of calls to the relative validation module functions and their intrinsic cost. In particular, the complexity of the training method is led by the sum of the costs of the chosen classification (O(classifier)) and clustering (O(clustering)) algorithms, which usually depend on the number of data samples and features. The complexity of the test method mainly depends on the complexity of the clustering algorithm in that only prediction is performed for the classifier. The random labeling algorithm increases the classifier complexity by a factor of n_rnd (i.e. the number of random labeling iterations). To perform cross validation and compute normalized stability we need to set the following parameters: n_C = |K| (i.e. the length of the list with the different number of clusters to try); n_fold and n_r, which corresponds to the number of cross-validation folds and repetitions respectively; and n_rnd. In conclusion, the cost of best_nclustcv is equal to n_C· n_fold· n_r (O(classifier) + O(clustering) + n_rnd O(classifier)) and the complexity of evaluate is O(classifier) + 2 · O(clustering).
  
  We run the methods considered on simulated datasets to empirically investigate the execution times. Simulations were run on a MacBook Pro 2020 with a 2.3 GHz Quad-Core Intel Core i7 processor and 32 GB of RAM. The complexity of clustering algorithms and state-of-the-art classifiers that can be used with reval can be found in Table <ref>. Figure <ref> and <ref> show the execution times of best_nclustcv and evaluate, respectively, on a two-blob dataset of 100 samples and 10 features with different combinations of classifiers and clustering algorithms. n_C is set at 5, n_fold and n_r are 2 and 10, respectively, and n_rnd is equal to 10. Default parameters were used for all classifiers and clustering algorithms. It is straightforward to observe that execution times largely depend on the chosen algorithms, with HDBSCAN the least expensive and spectral clustering the most expensive choice among clustering techniques, irrespective of classifier, and random forest the most expensive choice among classifiers. Overall the execution time of the evaluate method is reduced when compared to repeated cross validation. The package implementation also includes a multiprocessing feature which speeds up computations for repeated cross validation but not for the evaluation method, see Supplementary Figure <ref>, where 7 jobs are simultaneously run. Figure <ref> shows the execution times with varying number of samples and features when sequentially running best_nclustcv and evaluate with k-nearest neighbors (KNN) classifier and k-means clustering. We can observe a moderate increase for low number of samples (<10^3) with varying number of features and a steep growth for 10^3 samples. Corresponding execution times with multiprocessing can be observed in Supplementary Figure <ref>. 
  
  
  
  § TECHNICAL VALIDATION
   
  
  To investigate reval performance, first we investigate group detection in simulated non-overlapping Gaussian blobs, and the handwritten digits dataset from the MNIST digit recognition database (Modified National Institute of Standards and Technology database - <http://yann.lecun.com/exdb/mnist/>). Then, we leverage the  class to determine the combination of classifier and clustering algorithms that best identifies the number of classes for 19 different datasets from the UCI Machine Learning repository <cit.>.
  
  
  
   §.§ Blobs dataset
  
  To provide a simple example of how to use reval, we generated 5 isotropic Gaussian blobs for a total of 1000 samples with 2 features. To run the example refer to the code in the “Blobs dataset” section of the supplementary material. The dataset was first split into training and test sets (70/30%). Then we selected KNN with the number of neighbors equal to 15 and k-means with Euclidean distance, as suggested by the experiment reported in the “Algorithm selection” section. We then run the stability-based algorithm with a 10× 2 repeated cross validation framework, with 10 random labeling repetitions, and the number of clusters varying from 2 to 6. We evaluated the solution found on the test set and report ACC and MCC scores as metrics. Last, we report the AMI score to compare true labels and predicted labels on the test set. We also determined the best number of clusters maximizing and minimizing internal measures such as silhouette and Davies-Bouldin. We computed the metrics independently on the training and test sets and return the best number of clusters found with k-means. 
  
  The normalized stability for varying number of clusters in Figure <ref> shows that the model perfectly identified 5 clusters with 0.0 normalized stability. The comparison between true and clustering labels can be observed in Supplementary Figure <ref>. The ACC and MCC scores on the test set are equal to 1.0 and the AMI is equal to 1.0. The maximum silhouette score is 0.83 on both training and test sets with the correct number of clusters and AMI scores equal to 1.0. The Davies-Bouldin index results do not replicate between training and test sets. The index is equal to 0.23 with 4 clusters on the training set and to 0.23 with 5 clusters on the test set. AMI scores are 0.91 for training and 1.0 for testing.
  
  
  
   §.§ MNIST dataset
  
  For the real-world dataset example, we considered the handwritten digits MNIST dataset (<http://yann.lecun.com/exdb/mnist/>) that includes 70,000 samples corresponding to 28× 28 images of digits from 0 to 9. When flattened, each sample has 784 features. The dataset has 10 classes, with ∼7000 samples for each class.
  
  First, we split the dataset into training and test sets of 60,000 and 10,000 samples, respectively. Then, we preprocessed the dataset with UMAP <cit.> for dimensionality reduction with 2 components. We run HDBSCAN clustering (as suggested in <https://umap-learn.readthedocs.io/en/latest/clustering.html>) with different classifiers and selected the best configuration. In particular, we considered KNN with the number of neighbors equal to 30, and SVM, RF, LR with default parameters from scikit-learn library. The HDBSCAN algorithm was initialized with minimum samples equal to 10 and minimum cluster size equal to 200. The relative clustering validation procedure was run with 10 repetitions of 2-fold cross validation, and number of random labeling iterations equal to 10. Because HDBSCAN does not need the number of clusters specified a priori, the normalized stability is computed averaging the misclassification error over the solutions that return the same number of clusters. We selected the one that minimizes the normalized stability and the trained classifier applied to test set is the one trained on the corresponding cross-validation fold. Because we are not preselecting the number of clusters results might differ between training and test sets.
  
  The best classifier-clustering combination is RF and HDBSCAN, which correctly identifies 10 clusters with a misclassification error equal to 0.06 (0.06) on the training set (see Figure <ref>), and replicates on the test set, identifying 9 clusters with ACC equal to 0.86 and MCC equal to 0.86. The AMI score on the test set is equal to 0.87. When the experiment was run with internal validation measures, we obtained 10 clusters in training and 9 clusters during testing, with silhouette scores equal to 0.60 and 0.64, respectively, and Davies-Bouldin scores equal to 0.96 and 0.82, respectively. In Figure <ref>, we observe that internal measures fail to detect the actual digit classes during training, whereas reval with HDBSCAN and RF successfully identifies them (see Figures <ref>, <ref>). On test set the result is the same among all methods as demonstrated by equal AMI scores, see Supplementary Figure <ref>. Based on these results, we can speculate that clustering on a subset of data (N=35,000) better detects the digits classes.
  
  It is worth noting that the classifier performance highly depends on the clustering solution in that it tends to overfit to the training dataset. To guide the choice of a classifier, users should first select an appropriate clustering algorithm, then select the classifier that minimizes stability and that reports the best performance on the held-out dataset. In case of equal or very similar outcomes, algorithm complexity should guide the classifier selection.
  
  
  
   §.§ Algorithm selection
  
  We considered 19 datasets from the UCI Machine Learning repository <cit.> including the test set of the handwritten digits dataset that can be found in scikit-learn “toy datasets” - <https://scikit-learn.org/stable/datasets/index.html#toy-datasets>. In Table <ref> we report the dataset names along with the number of samples, features, and inherent classes. We applied the relative validation algorithm with different combinations of classifier and clustering algorithms. In particular, for clustering we selected: hierarchical clustering (HC) with Ward's method and Euclidean distance; k-means clustering with Euclidean distance; and HDBSCAN with minimum cluster size equal 5 and Euclidean distance. Among classifiers we opted for: KNN with 1 neighbor and Euclidean distance; random forest (RF) classifier with 100 estimators; support vector machines (SVM) with C=1 and γ=1/N samples; and logistic regression (LR). We did not consider spectral clustering because of its computational cost (see Table <ref>). To improve the performance, in addition to raw datasets, we repeatedly run the experiments after preprocessing with: 1) standard scaler, which removes the mean and scales to unit variance; 2) uniform manifold approximation and projection (UMAP) algorithm <cit.> for dimensionality reduction; and 3) standard scaler and UMAP. UMAP parameters are chosen according to those suggested in the documentation - <https://umap-learn.readthedocs.io/en/latest/clustering.html>. We ran the algorithm with 10 replications of 2-fold cross validation and 10 random labeling iterations. The number of clusters ranges from 2 to n_(classes)+2, where n_classes is the number of classes for each dataset.
  
  Table <ref> reports the best solutions selected as those reporting minimum stability along with the correct number of clusters. If no experiment identified the correct number of classes, we chose, among the solutions with minimum stability for each preprocessed dataset, the one with maximum adjusted mutual information (AMI), ACC, and MCC. AMI external score measures the similarity of two labelings of the same data, irrespective of label order and ranges from [0, 1], with a perfect match equal to 1. For comparison, we computed the silhouette score and Davies-Bouldin index <cit.> internal measures independently on both training and test sets with the clustering algorithm selected by the stability-based approach and with the same varying number of clusters. We report in Table <ref> the best solutions as those that maximize and minimize those measures respectively. Davies-Bouldin index <cit.> measures clusters separation and “tightness” with Euclidean distance. Lower indices correspond to better solutions. AMI scores are reported to evaluate the similarity of true and cluster labels on test sets.
  
  We observe that the stability-based approach identified the correct number of classes in 68% of cases (N=13). Moreover, 6 out 13 experiments selected k-means clustering and KNN classifier as the best choice. Of these, only 1 utilized raw data with no preprocessing, whereas the rest required either UMAP-preprocessed and/or scaled datasets. Because k-means works well with center-based spherical clusters and usually cannot find a good representation if clusters are very elongated or have complicated shapes, data preprocessing can relax this issue. From this experiment, it emerged that UMAP can be used as a preprocessing tool in this sense, in that it unwraps manifolds to find manifold boundaries. Nevertheless, because each dataset has its own intrinsic characteristics, preprocessing steps must be chosen with care to avoid breaking clusters into several erroneous small spherical clusters (see example with k-means at <https://umap-learn.readthedocs.io/en/latest/clustering.html>). It has been observed <cit.> that if the number of clusters k is too large with respect to the true clusters, the k-means algorithm tends to be unstable. On the contrary, if k is smaller or equal to the true number of clusters, the algorithm tends to be stable. Therefore, it is argued that k-means stability depends on the number of true clusters in the dataset, which should be on the order of 10 to provide stable solutions. This seems to hold when we look at the UCI datasets with ≥ 10 classes for which k-means is often selected as the best clustering algorithm although it returns smaller number of clusters respect to true classes. For this reason, despite KNN/k-means provides the best algorithm configuration in more than half the experiments, the choice of a classifier/clustering should be done carefully taking into account the dataset dimension and the computational cost of the algorithms.
  
  The number of clusters selected with the silhouette score and Davies-Bouldin index (see Table <ref>) is equal to that reported in Table <ref>. The exception here was the iris dataset <cit.>, whereby reval selects 3 clusters during training and 2 during testing with RF and HDBSCAN, whereas internal measures result in 4 and 2 cluster solutions. Nevertheless, the validation stability of the relative approach is equal to 0.33 (0.19), suggesting that the partition does not generalize well because the solution is not stable. If we could only rely on internal measures we would have failed to acknowledge the quality of the solution found. Generally, because 2 out of 3 classes are not linearly separable and the data is displayed into two separate groups <cit.> the iris dataset is not a good candidate for clustering and only the relative validation approach can clearly show that.
  
  More generally, when comparing the AMI scores for the internal measure results, we have no ability to say how well the results are actually doing. On the contrary, we have a sense of how good the clustering is with relative validation through the generalization process. See the example of the climate dataset, which has a validation normalized stability of 0.20 (0.04), indicating poor generalization, and a silhouette score on test set of 0.49, although both methods provide the same clustering solution on test set given that AMI score is equal to -0.005 in both cases. 
  
  In conclusion, grid-search for classifiers and clustering methods is easily and effectively implemented with reval and can be handy to avoid a priori selection of a classifier. Furthermore, the stability-based approach helps evaluate the goodness of a clustering solution by means of its generalization process. This information about generalization is absent with internal measures.
  
  
  
  § STABILITY REGIME TO GUIDE CLUSTER SELECTION
  
  Ben-Hur and colleagues presented a stability-based method with data sub-sampling and reported on the risk of underestimating/overestimating the true number of clusters <cit.>. Introducing prediction strength computed via repeated cross validation, Tibshirani et al. <cit.> linked unsupervised to supervised learning in the attempt to overcome the best cluster estimation issue. The reval implementation moves forward adding the evaluation of the best solution on unseen data, an approach that is particularly suitable for datasets with large sample size. While large sample sizes were less common in the past, dataset size has increased substantially over time and reval can be particularly important and well suited for modern data science contexts. Solution generalizability can be leveraged to select the best number of clusters not only based on validation metrics (e.g., prediction strength greater than 0.8/0.9 <cit.>), but also on the performance of the solution applied to a new set of data. In this way, we are able to compare test set distribution to the one on which the result was based, hence reinforcing the decision. Selecting the number of clusters that best generalizes to new data (i.e., investigation of the stability regime) holds promise for overcoming the underestimation issue. To give a better sense of how this works in practice we present a case study in the context of autism research, where clustering solutions within a stability regime could be further investigated. 
  
  Autism spectrum conditions (ASCs) are characterized by difficulties in social-communication alongside heightened restricted and repetitive behaviors <cit.>. The spectrum of affected individuals with ASC is highly heterogeneous and this heterogeneity is present at multiple levels, from genome to phenome, and can co-exist with differing levels of severity and comorbidities <cit.>. Given the high level of heterogeneity, data-driven clustering could be a promising approach to isolating different types of autisms. To split ASC into data-driven subtypes, we applied reval to clinical data obtained from the National Database for Autism Research (NDAR - <https://nda.nih.gov/>). NDAR is a database that includes a heterogeneous collection of de-identified human subjects data for autism research. We focus here on clinical behavioral data from the Vineland Adaptive Behavior Scales (VABS) <cit.>. Within the VABS, there are three domain total scores for communication, living skills, and socialization skills. Using these 3 domain total scores with UMAP preprocessing, we trained the stability-based model on 420 subjects (mean age 41.65 (17.91) months, female/male counts 109/311) . As for analysis choices, we ran this analysis with 100× 3 repeated cross validation and 100 iterations of random labeling, with number of clusters ranging from 2 to 10, using k-means clustering, and a KNN classifier with number of neighbors equal to 15. 
  
  From the performance plot in Figure <ref> we find that both a 2-cluster and 3-cluster solution results in small stabilities [2-cluster solution (error): 0.027 (0.004); 3-cluster solution (error): 0.036 (0.003)] and thus form a stability regime whereby either solution might be a promising solution to follow-up with future work. Based on the minimization of the normalized stability measure, the default behavior when using reval would be to select 2 as the best number of clusters. However, upon evaluation of these solutions on the unseen test set (n = 344 subjects; mean age 43.12 (17.04) months, female/male counts 90/254), we find that both solutions reach 94% accuracy, further confirming the presence of a stability regime whereby more than one solution might be a plausibly good model for follow-up work. The generalization performance alongside visual inspection (see Supplementary Figure <ref>) indicates that the default selection of a 2-cluster solution would possibly underestimate the true number of clusters <cit.>. Thus, based on the stability regime present here, we could select 3 as the true number of clusters since the stability differences are likely neglibgle and because both 2 and 3 cluster solutions generalize equally well. In practice, we would ultimately follow-up with examination of both 2 and 3 cluster solutions and utilize other datasets to better understand which solution might be most illuminating for decomposing clinical and biological heterogeneity of importance for autism research. 
  
  To contrast the reval stability regime results here with internal measures, we would have obtained 2 cluster solutions for both silhouette (scores of 0.41 and 0.42 on training and test respectively)  and Davies-Bouldin measures (0.86 on both training and test sets). If we force the number of clusters to 3 we obtain lower silhouette scores (i.e., 0.35 in both training and test sets), but higher Davies-Bouldin indices (i.e., 0.93 and 0.95, respectively). In this example, internal measures do not reveal a regime of possible cluster solutions. Given the additional lack of information about generalization from internal measures, such a regime might be easily missed. This example illustrates a real-world example in data science for how relative validation implemented with reval may reveal insights regarding regimes of clustering solutions that may be missed with internal validation approaches.
  
  
  
  § DISCUSSION
  
  In this work, we introduce the reval package for relative clustering validation and describe how it can be utilized as well as provide examples for how it performs in simulations and several real datasets. In many cases reval successfully identifies the correct number of clusters and confers several other advantages over and above other internal validation approaches. First, when compared to internal validation measures, reval is able to report the extent different clustering solutions fit to the data at hand and how well those solutions may generalize or replicate on unseen data. Moreover, because reval works with multiple clustering algorithms, it can facilitate a more thorough investigation of clustering mechanisms. In fact, although a thorough theoretical and experimental analysis of stability-based model selection with k-means clustering have been done <cit.>, the effectiveness of such an approach needs to be further investigated with different clustering algorithms. Furthermore, reval can be included in ensemble learning pipelines <cit.> or integrated in ensemble clustering frameworks for the selection of the best clustering solution <cit.>. Last, the ability to identify stability regimes and evaluate such regimes based on generalization to unseen data can inform the underestimation issue of the best number of clusters identified in real-world datasets. 
  
  A primary caveat or limitation to the approach reval takes is primarily one of data size. reval identifies the best clustering solution within a cross-validation framework, and hence needs large sample size to preserve cluster distribution between training and validation sets. Moreover, a separate held-out dataset is also needed to generalize the solution found. Smaller sized datasets may not allow for sufficient splitting within a cross-validation framework to allow for robust clustering solutions to generalize in unseen datasets. Finally, reval does not address the possibility of finding unrealistic data partitions. Because classifiers can overfit to their training set, a stable solution does not necessarily imply the true presence of subgroups in the data.  Future work will focus on the implementation of other relative validation methods, e.g., based on prediction strength <cit.> with the aim to create a comprehensive library.
  
  
  
  § COMPUTATIONAL DETAILS AND CODE AVAILABILITY
  
  Code, package installation instructions, and documentation with working examples, can be found at <https://github.com/IIT-LAND/reval_clustering>. The code to replicate the simulations presented in the “Technical validation” section is in . Code was written in Python 3.8 and simulations were run on a MacBook Pro 2020 with a 2.3 GHz Quad-Core Intel Core i7 processor and 32 GB RAM.
  
  
  
  § ACKNOWLEDGMENTS
  
  This project was supported by funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No 755816 (ERC Starting Grant to MVL).
  
  
  
  § AUTHOR CONTRIBUTIONS
  
  Conceptualization, I.L. and M.V.L; Methodology, I.L.; Software, I.L. and V.M.; Investigation,
  I.L. and V.M.; Writing – Original Draft, I.L. and M.V.L.; Writing –
  Review & Editing, I.L., V.M., and M.V.L.; Funding Acquisition, M.V.L.; Supervision, M.V.L.
  
  
  
  § DECLARATION OF INTERESTS
  
  The authors declare no competing interests.
  
  
  [title=References]
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  [pages=-]supplementary.pdf
  