Quick tour of OTclust
Lixiang Zhang and Beomseok Seo
2023-10-05
1. Introduction
OTclust is an R package for computing a mean partition
of an ensemble of clustering results by optimal transport alignment
(OTA) and for assessing uncertainty at the levels of both partition and
individual clusters. To measure uncertainty, set relationships between
clusters in multiple clustering results are revealed. Functions are
provided to compute the Covering Point Set (CPS), Cluster Alignment and
Points based (CAP) separability, and Wasserstein distance between
partitions.
2. Mean partition as an ensemble clustering.
library(OTclust)
data(sim1)
Here, we illustrate the usage of OTclust for an ensemble
clustering based on a simulated toy example, sim1, which
has 5000 samples, 2 features, and 4 clusters. ensemble( )
generates nbs number of perturbed partitions based on a
specified clustering method. For the clustering method, user specified
functions or example methods included in package (“kmeans”,
“Mclust”, “hclust”, “dbscan”, “PCAreduce”, “HMM-VB”) can be
used.
# the number of clusters.
C = 4
# generate an ensemble of perturbed partitions.
# if perturb_method is 1 then perturbed by bootstrap resampling, it it is 0, then perturbed by adding Gaussian noise.
ens.data = ensemble(sim1$X, nbs=100, clust_param=C, clustering="kmeans", perturb_method=1)
To find a consensus partition, the function otclust( )
searches mean partition by optimal transport alignment (OTA) between the
ensemble of partitions. As a return, otclust( ) generates
mean partition and its partition-wise and cluster-wise uncertainty
statistics. For the detail of return values, refer to help of
otclust( ).
# find mean partition and uncertainty statistics.
ota = otclust(ens.data)
# calculate baseline method for comparison.
kcl = kmeans(sim1$X,C)
# align clustering results for convenience of comparison.
compar = align(cbind(sim1$z,kcl$cluster,ota$meanpart))
lab.match = lapply(compar$weight,function(x) apply(x,2,which.max))
kcl.algnd = match(kcl$cluster,lab.match[[1]])
ota.algnd = match(ota$meanpart,lab.match[[2]])
# plot the result on two dimensional space.
otplot(sim1$X,sim1$z,con=F,title='Truth') # ground truth
#> Warning: replacing previous import 'lifecycle::last_warnings' by
#> 'rlang::last_warnings' when loading 'tibble'
#> Warning: replacing previous import 'lifecycle::last_warnings' by
#> 'rlang::last_warnings' when loading 'pillar'
otplot(sim1$X,kcl.algnd,con=F,title='Kmeans') # baseline method
otplot(sim1$X,ota.algnd,con=F,title='Mean partition') # mean partition by OTclust
3. Uncertainty assessment of clustering results
Here, as cluster-wise uncertainty measures, we briefly introduce the
usage of topological relationship statistics of mean partitions, cluster
alignment and points based (CAP) separability, and covering point sets
(CPS). The detailed definition of the above statistics can be found in
[1]. Moreover, if you want to carry out CPS Analysis, please next two
sections.
# distance between ground truth and each partition
wassDist(sim1$z,kmeans(sim1$X,C)$cluster) # baseline method
#> [1] 0.2506715
wassDist(sim1$z,ota$meanpart) # mean partition by OTclust
#> [1] 0.2520358
# Topological relationships between mean partition and ensemble clusters
t(ota$match)
#> C1 C2 C3 C4
#> match 99 90 89 90
#> split 0 0 0 0
#> merge 0 0 0 0
#> l.c. 1 10 11 10
# Cluster Alignment and Points based (CAP) separability
ota$cap
#> C1 C2 C3 C4
#> C1 0.0000000 1.0000000 0.9129447 0.9961823
#> C2 1.0000000 0.0000000 0.9992928 0.9524010
#> C3 0.9129447 0.9992928 0.0000000 1.0000000
#> C4 0.9961823 0.9524010 1.0000000 0.0000000
# Covering Point Set(CPS)
otplot(sim1$X,ota$cps[lab.match[[2]][1],],legend.labels=c('','CPS'),add.text=F,title='CPS for C1')
#> Warning: Removed 2 rows containing missing values (geom_text).
otplot(sim1$X,ota$cps[lab.match[[2]][2],],legend.labels=c('','CPS'),add.text=F,title='CPS for C2')
#> Warning: Removed 2 rows containing missing values (geom_text).
otplot(sim1$X,ota$cps[lab.match[[2]][3],],legend.labels=c('','CPS'),add.text=F,title='CPS for C3')
#> Warning: Removed 2 rows containing missing values (geom_text).
otplot(sim1$X,ota$cps[lab.match[[2]][4],],legend.labels=c('','CPS'),add.text=F,title='CPS for C4')
#> Warning: Removed 2 rows containing missing values (geom_text).
The red area of the above plots indicates covering point set (CPS)
for each cluster. The detail of the CPS analysis is addressed in the
next section.
4. CPS Analysis on selection of visualization methods
The functions that are going to be used in this section are
visCPS( ), mplot( ) and cplot(
). First, the function visCPS( ) is used for the
main computation of the CPS Analysis. The input should include: 1.
vlab, which is the visualization coordinates generated by
the visualization method that you are going to assess. 2.
ref, the true cluster labels of the samples. 3.
nEXP, which is optional, the number of perturbed results
for CPS Analysis, default 100. Larger the nEXP is, the
longer time it will take to compute.
# CPS analysis on selection of visualization methods
data(vis_pollen)
c=visCPS(vis_pollen$vis, vis_pollen$ref)
After the computation, we have the return list c, which would be the
input of function mplot( ) or cplot( ).
The mplot( ) will provide the membership heat map of the
required cluster, and the input should be c and the
cluster number. The cplot( ) will provide the covering
point set plot of the required cluster, and the input should be
c and the cluster number.
# visualization of the result
mplot(c,2)
cplot(c,2)
Furthermore, if you want to see the statistics, you can simply view
the return of visCPS( ):
# overall tightness
c$tight_all
#> [1] 0.5107356
# cluster-wise tightness
c$tight
#> 1 2 3 4 5 6
#> Tightness of each cluster 0.2024297 0.7028985 1 0.6016052 0.9138116 0.4588648
#> 7 8 9 10 11
#> Tightness of each cluster 0.4339253 0.2081705 0.1199151 0.4464706 0.5300001
5. CPS Analysis on validation of clustering result
In this section, the relevant functions are clustCPS(
), preprocess( ), perturb( ),
CPS( ), mplot( ), cplot( )
and pplot( ). For most of the users, you just need to use
the clustCPS( ) for the CPS Analysis. It will provide you
a lot of choice: For visualization method, you can choose between
tsne and umap. You can decide to add the
noise before or after the dimension reduction by parameter
noi. Also, you can choose to use Kmeans or
Mclust as the clustering method. Here is the example of a
single cell dataset, choosing to use the log transformation and
preprocessing based on the variance, which can reduce the initial
dimension of the data set. If you want to use other dimension reduction
technique or you need to carry out other preprocessing than we provide,
you just need to set l=FALSE, pre=FALSE, dimr=“None”, and
then input your processed result as parameter data.
# CPS Analysis on validation of clustering result
data(YAN)
y=clustCPS(YAN, k=7, l=FALSE, pre=FALSE, noi="after", cmethod="kmeans", dimr="PCA", vis="tsne")
#> Warning in min(ref): no non-missing arguments to min; returning Inf
#> Warning in if (class(X) == "dist") {: the condition has length > 1 and only the
#> first element will be used
#> sigma summary: Min. : 0.323162264525782 |1st Qu. : 0.686532727791371 |Median : 0.840637685950217 |Mean : 0.832540338898672 |3rd Qu. : 0.996223616580691 |Max. : 1.26695806934483 |
#> Epoch: Iteration #100 error is: 14.2048223360467
#> Epoch: Iteration #200 error is: 0.478304552870317
#> Epoch: Iteration #300 error is: 0.468634033027863
#> Epoch: Iteration #400 error is: 0.431538890961456
#> Epoch: Iteration #500 error is: 0.431229961683245
#> Epoch: Iteration #600 error is: 0.43122985521032
#> Epoch: Iteration #700 error is: 0.431229854930854
#> Epoch: Iteration #800 error is: 0.431229854929765
#> Epoch: Iteration #900 error is: 0.431229854929761
#> Epoch: Iteration #1000 error is: 0.431229854929761
# visualization of the results
mplot(y,4)
cplot(y,4)
# point-wise stability assessment
p=pplot(y)
p$v
If you want to try other clustering method rather than
Kmeans or Mclust, you will need to use the
function CPS( ). For this function, you need to input
several things. First, the reference clustering result, which might be
generated by your own clustering method. Second, the 2-dimension
visualization coordinates of your samples, which will be further used by
mplot( ) or cplot( ). Third, a collection
of clustering results in a matrix format, each column represents one
clustering result. To get this matrix, you might also want to use the
function perturb( ). Suppose the dataset you are going to
use for clustering is X, then perturb(X)
will give you a perturbed version of it. You can use this perturbed
version for clustering to get one clustering result. Repeat this for
several times, you will get a collection of clustering results.
References
[1] Jia Li, Beomseok Seo, and Lin Lin. “Optimal transport, mean
partition, and uncertainty assessment in cluster analysis.” Statistical
Analysis and Data Mining: The ASA Data Science Journal 12.5 (2019):
359-377.
[2] Lixiang Zhang, Lin Lin, and Jia Li. “CPS analysis: self-contained
validation of biomedical data clustering.” Bioinformatics 36.11 (2020):
3516-3521.