#  Selecting a suitable Parallel Label-propagation based algorithm for
Disjoint Community Detection

Subhajit Sahu
Center for Security, Theory, and Algorithmic Research (CSTAR)
IIIT Hyderabad, India - 500 032
<EMAIL_ADDRESS>
###### Abstract

Community detection is an essential task in network analysis as it helps
identify groups and patterns within a network. High-speed community detection
algorithms are necessary to analyze large-scale networks in a reasonable
amount of time. Researchers have made significant contributions in the
development of high-speed community detection algorithms, particularly in the
area of label-propagation based disjoint community detection. These algorithms
have been proven to be highly effective in analyzing large-scale networks in a
reasonable amount of time. However, it is important to evaluate the
performance and accuracy of these existing methods to determine which
algorithm is best suited for a particular type of network and specific
research problem. In this report, we investigate the RAK, COPRA, and SLPA,
three label-propagation-based static community discovery techniques. We pay
close attention to each algorithm’s minute details as we implement both its
single-threaded and multi-threaded OpenMP-based variants, making any necessary
adjustments or optimizations and obtaining the right parameter values. The RAK
algorithm is found to perform well with a tolerance of 0.05 and OpenMP-based
strict RAK with 12 threads was 6.75x faster than the sequential non-strict
RAK. The COPRA algorithm works well with a single label for road networks and
max labels of 4-16 for other classes of graphs. The SLPA algorithm performs
well with increasing memory size, but overall doesn’t offer a favourable
return on investment. The RAK algorithm is recommended for label-propagation
based disjoint community detection.

##  1 Introduction

There has been an unprecedented growth of real-world interconnected data [1]
that is naturally represented as graphs. These graphs emerging from real-world
applications such as social and biological networks are massive in size
leading to the adoption of parallelism as a way to address the scale issue. In
addition to being massive, most real-world graphs are dynamic in nature with
new edges being deleted and inserted all the time [2]. Therefore, parallel
algorithms for updating graph analytics over dynamic graphs are gaining
significant research attention in recent years. Examples include the dynamic
computation of centrality scores [3, 4, 5, 6, 7], dynamic maintenance of
biconnected components [8, 9, 10], dynamic computation of shortest paths [11,
12, 13].

Community detection is a popular graph problem with practical applications in
fields such as e-commerce, communication networks, and healthcare. It involves
identifying groups of vertices in a graph, known as communities, that
correspond to functional units in the system and are densely connected within
the community but sparsely connected to the rest of the network. When we can
identify these structures based on the topology of the network alone, they are
said to be intrinsic (studied here). On the other hand, extrinsic communities
are defined based on external information or attributes of the nodes, such as
membership in a particular organization or geographic location. There are
various types of communities that can be identified in a graph. Disjoint
communities (studied here) are those in which each vertex belongs to exactly
one community. Alternatively, overlapping communities allow each vertex to
belong to more than one community [14], while hierarchical communities have a
multi-level membership structure.

The problem of community detection is NP-hard, and unlike graph partitioning,
the number and size distribution of communities is not known beforehand. There
are various methods for community detection, including label propagation [15,
14, 16], random walk [17], diffusion [18], spin dynamics [19], fitness metric
optimization [20, 21], statistical inference [22, 23], simulated annealing
[24, 19], clique percolation [25, 26, 27], and more. These methods can be
divided into two main groups: divisive and agglomerative. Divisive methods,
also known as top-down methods, start with all the vertices in a graph being
part of a single community and iteratively identify and remove bridges to
split into more well-connected communities [28, 29]. Agglomerative methods, or
bottom-up methods, merge two or more communities together such that a certain
score is maximized [30, 2]. In certain cases, we have a set of relevant seed
vertices of interest which are expanded to form communities surrounding them.
This is known as seed set expansion [31].

Different community detection methods can return different communities, which
can be evaluated using a quality function. Modularity is a popular fitness
function introduced by Newman and Girvan which compares the number of intra-
community edges to the expected number in a random-null model. It lies between
-0.5 (non-modular clustering) and 1.0 (fully modular clustering). Optimizing
this function theoretically results in the best possible grouping. Conductance
is another fitness score that measures the community cut or inter-community
edges. Constant Potts Model (CPM) is another quality function that overcomes
some limitations of modularity [32].

##  2 Literature survey

There are several techniques that have been developed for detecting
communities in networks. A number of them are based on modularity-
optimization, hierarchical clustering, label propagation, region density, core
clustering [33], game theory, information theory (infomap) [34, 35, 36, 37,
38], and biological evolution (genetics) [39, 40, 41]. Metrics such as the
modularity score [20, 42, 40], Normalized Mutual Information index (NMI) [43,
44], and Jaccard Index [43] are used to compare the quality of communities
obtained with different approaches.

The Raghavan Albert Kumara (RAK), also commonly known as the Label Propagation
Algorithm (LPA), is a method used for identifying communities or groups within
a network by initializing each vertex with a unique label and diffusing these
labels across the graph. LPA is faster and more scalable than the Louvain
algorithm, as it can be parallelized easily and does not require repeated
optimization steps. However, it may not be as effective at identifying high-
quality community structure in some types of networks and tends to produce
communities with lower modularity than the Louvain algorithm. Additionally,
the resulting solution produced by LPA may not be unique and may contain
monster communities [45, 15].

Improvements upon the LPA include the Node Influence Based Label Propagation
Algorithm (NIBLPA) while provides a stable (non-random) mechanism of label
choosing in the case of multiple best labels [46]; Label influence Policy for
Label Propagation Algorithm (LP-LPA) and Spreading Activation Label
Propagation Algorithm (SALPA) [47, 48] which address the issue of monster
communities; the Three Stage (TS) algorithm which first identifies central
nodes and combines communities for improved modularity (no information has
been provided on its execution time) [49]; and Lie et al.’s approach which
uses frontiers with alternating push-pull to reduce the number of edges
visited and improve solution quality [50]. Existing approaches for finding
overlapping communities, such as the Community OVerlap PRopagation Algorithm
(COPRA), or the Speaker-Listener Label Propagation Algorithm (SLPA), may also
be used for disjoint communities by choosing the best label associated with
each vertex [14, 16]. While these modifications to LPA achieve communities of
better quality, they can be (in some cases significantly) slower than the
original LPA. A GPU-accelerated parallel implementation of the original LPA is
available that is able to deal with large-scale datasets that do not fit into
GPU memory [51].

##  3 Evaluation

In this section, we first describe our experimental setup, such as the system
we use and our dataset. We then investigate various static community detection
algorithms, namely, the RAK algorithm, the COPRA algorithm, and the SLPA
algorithm. Here we take note of the subtle details and explore various
optimizations of each of the algorithms, while implementing their single-
threaded and multi-threaded OpenMP-based versions.

###  3.1 Experimental setup

####  3.1.1 System used

For our experiments, we use a system consisting of two Intel Xeon Silver 4116
64-bit @ 2.10 GHz Server CPUs, and 128GB DDR4 Synchronous Registered DRAM @
2666 MHz. Each CPU has 12 x86 cores (with 2 hyper-threads per core) and 16.5M
L3 cache. Our server is running CentOS version 7.9. We use GCC version 9.3 and
OpenMP version 5.0 to compile with optimization level 3 (-O3). For all
experiments we keep simultaneous multi-threading (SMT) enabled.

####  3.1.2 Dataset

The graphs we use in our experiments are shown in Table 1. All of them are
obtained from the SuiteSparse Matrix Collection [52]. The total number of
vertices in the graphs vary from $74.9$ thousand to $12.0$ million, and the
total number of edges vary from $811$ thousand to $304$ million. We consider
all edges to be undirected and weighted with a default weight of one, and add
self-loops to each vertex in all the graphs.

Table 1: List of 17 graphs used in our experiments. All edges of each graph were duplicated in the reverse direction, thus making them undirected, and the weight of each edge was set to 1. Here, $|V|$ is the total number of vertices, $|E|$ is the total number of edges (after making the graph undirected), and $D_{avg}$ is the average degree of vertices. The number of vertices and edges are rounded to the nearest thousand (K) or million (M), as appropriate. Graph | $|V|$ | $|E|$ | $D_{avg}$
---|---|---|---
Web Graphs
web-Stanford | 282K | 3.99M | 14.1
web-BerkStan | 685K | 13.3M | 19.4
web-Google | 916K | 8.64M | 9.43
web-NotreDame | 326K | 2.21M | 6.78
indochina-2004 | 7.41M | 304M | 41.0
Social Networks
soc-Slashdot0811 | 77.4K | 1.02M | 13.2
soc-Slashdot0902 | 82.2K | 1.09M | 13.3
soc-Epinions1 | 75.9K | 811K | 10.7
soc-LiveJournal1 | 4.85M | 86.2M | 17.8
Collaboration Networks
coAuthorsDBLP | 299K | 1.96M | 6.56
coAuthorsCiteseer | 227K | 1.63M | 7.18
coPapersCiteseer | 434K | 32.1M | 74.0
coPapersDBLP | 540K | 30.5M | 56.5
Road Networks
italy_osm | 6.69M | 14.0M | 2.09
great-britain_osm | 7.73M | 16.3M | 2.11
germany_osm | 11.5M | 24.7M | 2.15
asia_osm | 12.0M | 25.4M | 2.12

###  3.2 RAK algorithm

Raghavan Albert Kumara (RAK), also known as LPA, is a popular label-
propagation based community detection algorithm. Here, every vertex is
initialized with a unique label and at every step each vertex adopts the label
that most of its neighbors currently have. In the weighted version of this
algorithm (which we use), each vertex adopts the label with the most
interconnecting weight. Through this iterative process densely connected
groups of vertices form a consensus on a unique label to form communities.

When there exist multiple communities (labels) connected to a vertex with
maximum weight, we randomly pick one of them (non-strict approach), or pick
only the first of them (strict approach). The algorithm converges when at
least $100-n\%$ of vertices don’t change their community membership
(`tolerance` parameter). We conduct experiments to obtain a `tolerance`
parameter value that would be suitable on average for most graphs. In
addition, we are interested in comparing the performance of both the non-
strict and the strict approaches, in terms of modularity (quality) and
computation time.

For the first experiment we implement a single-threaded CPU-based version of
the RAK algorithm 111Single-threaded RAK, https://github.com/ionicf/rak-
communities-seq. We adjust `tolerance` for both the non-strict and strict
approaches from $0.1$ to $0.0001$. From the results, we make the following
observations. Non-strict approach generally achieves better modularity, with
certain exceptions (`coAuthors*`, road networks). We think this has to do with
graph structure. It seems to us that a `tolerance` of $0.05$ would be a good
choice.

For the second experiment, we implement a multi-threaded OpenMP-based version
of the RAK algorithm 222Multi-threaded OpenMP-based RAK,
https://github.com/ionicf/rak-communities-openmp. Each thread is given a
separate hashtable, which it can use for choosing the most popular label among
its neighbors (by weight). We initially packed the hashtables (for each
thread) contiguously on a vector. Nevertheless, we observe that allocating
them separately gives us an almost 2x performance improvement over the
contiguous allocation version (by using pointers to vectors). We use an OpenMP
schedule of “`auto`” and a total of $12$ threads for the time being. For
optimal performance, a better schedule can be chosen through
experimentation/profiling, and up to $N-1$ threads put in use (where $N$ is
the number of cores/threads). Similar to the previous experiment, we adjust
`tolerance` from $0.1$ to $0.0001$ and compare the sequential and OpenMP
implementations (non-strict, strict approaches).

From the results of the second experiment, we observe that OpenMP-based strict
RAK performs better on average, both in terms of time and modularity. In
addition, we generally see it achieve better modularity than sequential
approaches (this could be due to better randomization). For a `tolerance` of
$0.05$, OpenMP-based strict RAK (using 12 threads) is 6.75x faster than
sequential non-strict RAK. However, OpenMP implementations do not achieve
better modularity on `coAuthor*` graphs. On social networks, OpenMP-based non-
strict RAK achieves better modularity than the strict version. Similar to
sequential RAK, we find a `tolerance` of $0.05$ to be a good choice. You can
find code and results for these experiments at our repository.

###  3.3 COPRA algorithm

Community Overlap PRopagation Algorithm (COPRA) is another label-propagation
based community detection algorithm. Unlike RAK, this algorithm uses multiple
labels per vertex, with each label having an associated belonging coefficient
(which sums to $1$). The algorithm is as follows. Each vertex initializes as
its own community (with belonging equal to $1$). At the start of each
iteration, each vertex collects labels from its neighborhood. Here, we exclude
each vertex’s own labels, although not explicitly mentioned in the original
paper. The collected labels are then scaled by interconnecting edge weights
for weighted graphs. Labels above a certain threshold are only picked. This
threshold is inversely proportional to the maximum number of labels allowed
(`max_labels` parameter). If all labels are below threshold, a random best
label is picked. We make a vertex join its own community if it has no labels.
Selected labels are then normalized such that the belonging coefficient sums
to $1$. This process is repeated until convergence.

The authors of this paper suggest using a minimum vertices per community count
to detect convergence. Alas, we do not find it to be helpful. Instead, we use
a convergence condition similar to that used in RAK, i.e., to count the number
of vertices that change their best label. Once this count falls below a
certain fraction (`tolerance` parameter), we consider the algorithm to have
converged. The authors also mention using a synchronous version of the
algorithm, where labels of each vertex are dependent only upon labels in
previous iteration. Nonetheless, we find an asynchronous approach to converge
faster (labels of each vertex can be dependent upon labels of its neighbors in
the current iteration). Since our focus is on finding disjoint communities, we
consider the best label of each vertex as the final result.

We conduct experiments to obtain a `tolerance` parameter value that would be
suitable on average for most graphs, both in terms of modularity and time. In
addition, we adjust the `max_labels` per vertex which in turn controls the
threshold weight above which labels are picked.

For the first experiment we implement a single-threaded CPU-based version of
the COPRA algorithm 333Single-threaded COPRA, https://github.com/ionicf/copra-
communities-seq. We adjust `tolerance` from $0.1$ to $0.0001$, and
`max_labels` from $1$ to $32$. From the results, we make the following
observations. On average, it seems that using a single label is best in terms
of time as well as modularity. This is indeed true in case of road networks,
but does not apply to the other classes of graphs (e.g. web graphs, social
networks, and collaboration networks). For example, web graphs such as `web-
Stanford` and `web-BerkStan` achieve best modularity with `max_labels` of $8$,
`web-Google` attains best modularity with `max_labels` of $32$, and `web-
NotreDame` does best with `max_labels` of $16$. `max_labels` of $4$-$16$ would
be a good choice for such graphs. In addition it seems that on average, making
the `tolerance` tighter than $0.01$ has no beneficial effect on modularity.
Having said that, tightening `tolerance` (to any value lower than $0.1$) does
not help with graphs such as `web-NotreDame`, `coAuthorsDBLP`, and social
networks. In any case, a `tolerance` of $0.01$ would be a good choice on
average.

For the second experiment, we implement a multi-threaded OpenMP-based version
of the COPRA algorithm 444Multi-threaded OpenMP-based COPRA,
https://github.com/ionicf/copra-communities-openmp. Similar to RAK, each
thread is given a separate hashtable, which it can use for choosing a set of
labels from its neighbors above a certain threshold (by weight). If no label
is above threshold, we pick only the best weighted label, and if not we simply
join our own community. Hashtables are allocated separately for better
performance (instead of storing them contiguously on a vector). Again, we use
an OpenMP schedule of “`auto`” and a total of $12$ threads for the time being
(not an optimal choice). Similar to the previous experiment, we adjust
`tolerance` from $0.1$ to $0.0001$, and `max_labels` from $1$ to $32$. From
the results, we observe that OpenMP-based approaches are faster than
sequential, and also seem to achieve better modularity. The rest of our
observations are similar to that of the first experiment.

###  3.4 SLPA algorithm

Speaker-listener Label Propagation Algorithm (SLPA) is yet another label-
propagation based community detection algorithm. Here, each vertex has a fixed
label memory where it remembers all the popular labels that it has listened
to. The size of this memory is fixed to the total number of iterations to be
performed plus one. The algorithm is as follows. Each vertex is first
initialized such that it remembers itself as popular. At the start of each
iteration, each neighbor (of a vertex) speaks one of the random labels in its
memory to the vertex. The vertex (listener) adds the most popular label to its
memory. This is repeated until a fixed number of iterations is performed,
equal to `memory_size` minus one. In order to accelerate the algorithm, we
allow early convergence if at least $100-n\%$ of vertices memorize their
previous label (`tolerance` parameter). We use a `tolerance` of $0.05$. As we
are interested in finding disjoint communities, we pick the most popular label
in the memory of each vertex as its community.

Unlike RAK and COPRA, this is a randomized historical-label based technique. A
random number generator (RNG) is used in order for each neighbor of a given
vertex to pick a random label to speak from its memory. We originally used
C++’s default random engine (probably Mersenne twister) which seemed slow. We
updated it to use a xorshift32 random engine. The listener vertex listens to
all of its neighbors, and based on edge weight, picks the most popular label.
In strict mode, the listener picks the first most popular label to record in
its memory, and in non-strict mode, the listener picks randomly one of the
most popular labels to record. It records this in the next index in its
memory.

We conduct an experiment, where we adjust the `memory_size` per vertex which
in turn controls the maximum number of iterations needed to be performed
555Single-threaded SLPA, https://github.com/ionicf/slpa-communities-seq. This
is done in order to obtain a `memory_size` parameter value that would be
suitable on average for most graphs, both in terms of modularity and time.

From the results, we observe that increasing the `memory_size` per vertex
increases the time required for completion, while also increasing the final
modularity. Anyhow, on graphs `soc-Slashdot*` modularity seems to decrease
with increasing `memory_size`. We also observe that the strict variant is
slightly faster than the non-strict one, while achieving similar modularity.
On road networks, the strict approach seems to obtain lower modularity, while
on other classes of graphs, the strict approach achieves higher modularity.
However, it seems that the SLPA algorithm does not provide a good return on
investment (it takes longer but does not achieve great modularity). It should
be noted that we do not perform any post-processing of communities (such as
splitting disconnected communities).

##  4 Conclusion

We examine a number of label-propagation-based static community detection
algorithms, including the RAK, COPRA, and SLPA. While implementing their
single-threaded and multi-threaded OpenMP-based versions, we pay attention to
the minute details of each algorithm, make any useful tweaks/optimizations,
and obtain suitable parameter values. The RAK algorithm performs well with a
`tolerance` of $0.05$ (which determines the maximum percentage of vertices
that change their community membership), with strict RAK performing better
than non-strict version in parallel implementation (except on social
networks), and OpenMP-based strict RAK with 12 threads was 6.75x faster than
the sequential non-strict RAK. The COPRA algorithm works well with a single
label for road networks, and `max_labels` of $4$-$16$ for other classes of
graphs. We use an asynchronous version of the algorithm (which converges
faster), and use the same convergence condition as the RAK algorithm. A
`tolerance` of $0.01$ is a good choice on average, and OpenMP-based
implementation is faster and achieves better modularity than sequential. For
the SLPA algorithm, we use a xorshift32 random engine for generation of random
numbers, and again use the same convergence condition as the RAK algorithm.
The strict variant of SLPA performs better than the non-strict one (except on
road networks), performs well with increasing `memory_size`, but overall
doesn’t offer a favourable return on investment (its slow). We therefore
recommend the reader to stick to the RAK algorithm for label-propagation based
disjoint community detection.

## References

  * [1] S. Sakr, A. Bonifati, H. Voigt, A. Iosup, K. Ammar, R. Angles, W. Aref, M. Arenas, M. Besta, P. A. Boncz, and et al., “The future is big graphs,” _Communications of the ACM_ , vol. 64, no. 9, pp. 62–71, Sep 2021.
  * [2] N. Zarayeneh and A. Kalyanaraman, “Delta-Screening: A Fast and Efficient Technique to Update Communities in Dynamic Graphs,” _IEEE transactions on network science and engineering_ , vol. 8, no. 2, pp. 1614–1629, Apr 2021.
  * [3] K. Lerman, R. Ghosh, and J. H. Kang, “Centrality metric for dynamic networks,” in _Proceedings of the Eighth Workshop on Mining and Learning with Graphs_ , 2010, pp. 70–77.
  * [4] C. Li, J. Han, G. He, X. Jin, Y. Sun, Y. Yu, and T. Wu, “Fast computation of simrank for static and dynamic information networks,” in _Proceedings of the 13th International Conference on Extending Database Technology_ , 2010, pp. 465–476.
  * [5] E. Bergamini and H. Meyerhenke, “Approximating betweenness centrality in fully dynamic networks,” _Internet Mathematics_ , vol. 12, no. 5, pp. 281–314, 2016.
  * [6] E. Nathan and D. A. Bader, “A dynamic algorithm for updating katz centrality in graphs,” in _Proceedings of the 2017 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining 2017_ , 2017, pp. 149–154.
  * [7] S. C. Regunta, S. H. Tondomker, K. Shukla, and K. Kothapalli, “Efficient parallel algorithms for dynamic closeness-and betweenness centrality,” _Concurrency and Computation: Practice and Experience_ , p. e6650, 2021.
  * [8] Z. Galil and G. F. Italiano, “Maintaining biconnected components of dynamic planar graphs,” in _International Colloquium on Automata, Languages, and Programming_. Springer, 1991, pp. 339–350.
  * [9] W. Liang, R. Brent, and H. Shen, “Fully dynamic maintenance of k-connectivity in parallel,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 12, no. 8, pp. 846–864, 2001.
  * [10] C. A. Haryan, G. Ramakrishna, K. Kothapalli, and D. S. Banerjee, “Shared-memory parallel algorithms for fully dynamic maintenance of 2-connected components,” in _2022 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 2022, pp. 1195–1205.
  * [11] P. Narvaez, K.-Y. Siu, and H.-Y. Tzeng, “New dynamic algorithms for shortest path tree computation,” _IEEE/ACM Transactions On Networking_ , vol. 8, no. 6, pp. 734–746, 2000.
  * [12] E. P. Chan and Y. Yang, “Shortest path tree computation in dynamic graphs,” _IEEE Transactions on Computers_ , vol. 58, no. 4, pp. 541–557, 2008.
  * [13] A. Khanda, S. Srinivasan, S. Bhowmick, B. Norris, and S. K. Das, “A parallel algorithm template for updating single-source shortest paths in large-scale dynamic networks,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 33, no. 4, pp. 929–940, 2021.
  * [14] S. Gregory, “Finding overlapping communities in networks by label propagation,” _New Journal of Physics_ , vol. 12, p. 103018, 10 2010. [Online]. Available: https://dx.doi.org/10.1088/1367-2630/12/10/103018
  * [15] U. N. Raghavan, R. Albert, and S. Kumara, “Near linear time algorithm to detect community structures in large-scale networks,” _Physical Review E_ , vol. 76, no. 3, Sep 2007.
  * [16] J. Xie, B. K. Szymanski, and X. Liu, “SLPA: Uncovering overlapping communities in social networks via a speaker-listener interaction dynamic process,” in _2011 ieee 11th international conference on data mining workshops_. IEEE, 2011, pp. 344–349.
  * [17] M. Rosvall and C. T. Bergstrom, “Maps of random walks on complex networks reveal community structure,” _Proceedings of the national academy of sciences_ , vol. 105, no. 4, pp. 1118–1123, 2008.
  * [18] K. Kloster and D. F. Gleich, “Heat kernel based community detection,” in _Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining_ , 2014, pp. 1386–1395.
  * [19] J. Reichardt and S. Bornholdt, “Statistical mechanics of community detection,” _Physical review E_ , vol. 74, no. 1, p. 016110, 2006.
  * [20] M. E. Newman, “Finding community structure in networks using the eigenvectors of matrices,” _Physical review E_ , vol. 74, no. 3, p. 036104, 2006.
  * [21] S. Fortunato, “Community detection in graphs,” _Physics reports_ , vol. 486, no. 3-5, pp. 75–174, 2010.
  * [22] E. Côme and P. Latouche, “Model selection and clustering in stochastic block models based on the exact integrated complete data likelihood,” _Statistical Modelling_ , vol. 15, pp. 564–589, 3 2015, doi: 10.1177/1471082X15577017. [Online]. Available: https://doi.org/10.1177/1471082X15577017
  * [23] M. E. Newman and G. Reinert, “Estimating the number of communities in a network,” _Physical review letters_ , vol. 117, no. 7, p. 078301, 2016.
  * [24] R. Guimera and L. A. Nunes Amaral, “Functional cartography of complex metabolic networks,” _nature_ , vol. 433, no. 7028, pp. 895–900, 2005.
  * [25] I. Derényi, G. Palla, and T. Vicsek, “Clique percolation in random networks,” _Physical review letters_ , vol. 94, no. 16, p. 160202, 2005\.
  * [26] S. Maity and S. K. Rath, “Extended Clique percolation method to detect overlapping community structure,” in _2014 International Conference on Advances in Computing, Communications and Informatics (ICACCI)_. IEEE, 2014, pp. 31–37.
  * [27] S. K. Gupta, D. P. Singh, and J. Choudhary, “A review of clique-based overlapping community detection algorithms,” _Knowledge and Information Systems_ , pp. 1–36, 2022.
  * [28] M. Girvan and M. E. Newman, “Community structure in social and biological networks,” _Proceedings of the national academy of sciences_ , vol. 99, no. 12, pp. 7821–7826, 2002.
  * [29] S. Souravlas, A. Sifaleras, M. Tsintogianni, and S. Katsavounis, “A classification of community detection methods in social networks: a survey,” _International journal of general systems_ , vol. 50, no. 1, pp. 63–91, Jan 2021.
  * [30] A. Zakrzewska and D. A. Bader, “A dynamic algorithm for local community detection in graphs,” in _2015 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining (ASONAM)_. IEEE, 2015, pp. 559–564.
  * [31] J. J. Whang, D. F. Gleich, and I. S. Dhillon, “Overlapping community detection using seed set expansion,” in _Proceedings of the 22nd ACM international conference on Information & Knowledge Management_, 2013, pp. 2099–2108.
  * [32] V. A. Traag, P. Van Dooren, and Y. Nesterov, “Narrow scope for resolution-limit-free community detection,” _Physical Review E_ , vol. 84, no. 1, p. 016114, 2011.
  * [33] Y. Ruan, D. Fuhry, J. Liang, Y. Wang, and S. Parthasarathy, _Community discovery: simple and scalable approaches_ , ser. Human–Computer Interaction Series. Springer International Publishing, 2015, pp. 23–54.
  * [34] J. Zeng and H. Yu, “Effectively Unified optimization for Large-scale Graph Community Detection,” in _2019 IEEE International Conference on Big Data (Big Data)_. IEEE, Dec 2019, pp. 475–482.
  * [35] ——, “A Distributed Infomap Algorithm for Scalable and High-Quality Community Detection,” in _Proceedings of the 47th International Conference on Parallel Processing - ICPP 2018_. New York, New York, USA: ACM Press, Aug 2018, pp. 1–11.
  * [36] M. A. M. Faysal and S. Arifuzzaman, “Distributed Community Detection in Large Networks using An Information-Theoretic Approach,” in _2019 IEEE International Conference on Big Data (Big Data)_. IEEE, Dec 2019, pp. 4773–4782.
  * [37] M. Rosvall, D. Axelsson, and C. T. Bergstrom, “The map equation,” _The European Physical Journal Special Topics_ , vol. 178, no. 1, pp. 13–23, Nov 2009\.
  * [38] L. Rita, “Infomap Algorithm. An algorithm for community finding,” Apr 2020. [Online]. Available: https://towardsdatascience.com/infomap-algorithm-9b68b7e8b86
  * [39] T. B. D. P. Aditama and A. Sn, “Determining Community Structure and Modularity in Social Network using Genetic Algorithm,” _Indonesian Journal of Computing and Cybernetics Systems_ , vol. 14, pp. 219–230, 2020.
  * [40] A. K. Ghoshal, N. Das, S. Bhattacharjee, and G. Chakraborty, “A fast parallel genetic algorithm based approach for community detection in large networks,” in _2019 11th International Conference on Communication Systems & Networks (COMSNETS)_. IEEE, Jan 2019, pp. 95–101.
  * [41] Y.-C. Lu and G. Chakraborty, “Improving Efficiency of Graph Clustering by Genetic Algorithm Using Multi-Objective Optimization,” _International Journal of Applied Science and Engineering_ , vol. 17, no. 2, pp. 157–173, Jun 2020.
  * [42] V. D. Blondel, J.-L. Guillaume, R. Lambiotte, and E. Lefebvre, “Fast unfolding of communities in large networks,” _Journal of Statistical Mechanics: Theory and Experiment_ , vol. 2008, no. 10, p. P10008, Oct 2008.
  * [43] A. Jain, R. Nasre, and B. Ravindran, “DCEIL: Distributed Community Detection with the CEIL Score,” in _2017 IEEE 19th International Conference on High Performance Computing and Communications; IEEE 15th International Conference on Smart City; IEEE 3rd International Conference on Data Science and Systems (HPCC/SmartCity/DSS)_. IEEE, Dec 2017, pp. 146–153.
  * [44] P. Chopade and J. Zhan, “A Framework for Community Detection in Large Networks Using Game-Theoretic Modeling,” _IEEE Transactions on Big Data_ , vol. 3, no. 3, pp. 276–288, Sep 2017.
  * [45] M. E. J. Newman, “Detecting community structure in networks,” _The European Physical Journal B - Condensed Matter_ , vol. 38, no. 2, pp. 321–330, Mar 2004.
  * [46] Y. Xing, F. Meng, Y. Zhou, M. Zhu, M. Shi, and G. Sun, “A node influence based label propagation algorithm for community detection in networks,” _The Scientific World Journal_ , vol. 2014, 2014.
  * [47] K. Berahmand and A. Bouyer, “LP-LPA: A link influence-based label propagation algorithm for discovering community structures in networks,” _International Journal of Modern Physics B_ , vol. 32, no. 06, p. 1850062, mar 2018. [Online]. Available: http://www.worldscientific.com/doi/abs/10.1142/S0217979218500625
  * [48] M. Sattari and K. Zamanifar, “A spreading activation-based label propagation algorithm for overlapping community detection in dynamic social networks,” _Data & knowledge engineering_, vol. 113, pp. 155–170, jan 2018. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/{S0169023X16303500}
  * [49] X. You, Y. Ma, and Z. Liu, “A three-stage algorithm on community detection in social networks,” _Knowledge-Based Systems_ , vol. 187, p. 104822, 2020\.
  * [50] X. T. Liu, M. Halappanavar, K. J. Barker, A. Lumsdaine, and A. H. Gebremedhin, “Direction-optimizing label propagation and its application to community detection,” in _Proceedings of the 17th ACM International Conference on Computing Frontiers_. New York, NY, USA: ACM, May 2020, pp. 192–201.
  * [51] Y. Kozawa, T. Amagasa, and H. Kitagawa, “GPU-Accelerated Graph Clustering via Parallel Label Propagation,” in _Proceedings of the 2017 ACM on Conference on Information and Knowledge Management - CIKM ’17_. New York, New York, USA: ACM Press, Nov 2017, pp. 567–576.
  * [52] S. Kolodziej, M. Aznaveh, M. Bullock, J. David, T. Davis, M. Henderson, Y. Hu, and R. Sandstrom, “The SuiteSparse matrix collection website interface,” _The Journal of Open Source Software_ , vol. 4, no. 35, p. 1244, Mar 2019\.


# A Proposed Method for Assessing Cluster Heterogeneity

John W. G. Addy  Malaria Immunology Laboratory, Francis Crick Institute,
London, NW1 1AT, UK  Jean Langhorne  Malaria Immunology Laboratory, Francis
Crick Institute, London, NW1 1AT, UK

###### Abstract

Assessing how adequate clusters fit a dataset and finding an optimum number of
clusters is a difficult process. A membership matrix and the degree of
membership matrix is suggested to determine the homogeneity of a cluster fit.
Maximisation of the ratio of the overall degree of membership at cluster
number lag 1 is also suggested as a method to optimise the number of clusters
in a dataset. A threshold factor upon the degree of membership is also
suggested for homogeneous clusters. Cluster simulations were given to compare
how well the proposed method compares against established methods. This method
may be applied to the output of both hierarchical and k-means clustering.

##  1 Introduction

Optimising the number of clusters through the partitioning of a dataset is a
difficult process, especially with a dataset of heterogeneous or nested
clusters. Methods of optimising the number of clusters in multivariate
analyses typically consider how well the cluster optimisation fits given the
within-cluster-sum-of-squares ($W_{k}$), by observing an elbow in $W_{k}$, as
cluster number ($k$) increases. However, the $W_{k}$, or manipulations of it,
may not consider the proximity of an observation to another cluster it has not
been assigned to. Identifying this on the individual basis has been given and
suggested through the use of fuzzy clustering (Bezdek, 1974; Dunn, 1973),
calculating the degree of closeness of an observation $i$ belonging to cluster
$k$. However, the membership of observation $i$ belonging to cluster $k$ may
only inform about a single observation and not all observations within a
cluster. Proposed is a method of determining the homogeneity of clusters given
their membership. This method does not consider which type of clustering that
has been used, but only the final cluster list and the associated distance
matrix of all individuals.

An empirical estimation of where an elbow occurs in the $W_{k}$ as cluster
number increases is given by the GAP statistic and the 1-standard-error rule
(Tibshirani et al., 2001). Other methods of optimising cluster number include
the comparison of the between-cluster-sum-of-squares along with the $W_{k}$
(Caliński and Harabasz, 1974) and the average Silhouette method (Rousseeuw,
1987). However, by ignoring the heterogeneity of the overall cluster fit,
through $W_{k}$, it may not be understood how well represented each cluster is
given the overall membership of observations. The Silhouette method
(Rousseeuw, 1987) provides some information about cluster membership, by
averaging over all calculated Silhouette values from the nearest neighbouring
cluster. Once the homogeneity of clusters has been considered, this allows the
identification of the degree of membership of a cluster, given the membership
of all clusters. The maximising of the $\phi$ ratio of such degrees of
membership at cluster lag 1 is suggested as an estimation of the optimum
number of clusters.

Demonstrated is the application of the proposed method on two datasets,
followed by a simulation study, similar to Tibshirani et al. (2001), showing
the utility of the method compared to other established methods. A short
review of other methods of optimising cluster classification has been given in
an effort to show how this method may be used in conjunction with others
methods. Hierarchical clustering was considered in this paper, although the
methods proposed may be considered for other clustering methods.

##  2 Distance Matrix and Membership

Consider a dataset ($x_{ip}$) over $i$ observations ($i=1,2,...,n$) and $p$
variates. Where the distance between each observation can be expressed as a
squared distance matrix ($d^{2}_{ij}$), derived by the squared Euclidean
distance between two observations

$d^{2}_{ij}=\sum_{p}(x_{ip}-x_{jp})^{2}.$

The aim of cluster analysis is to classify $x_{ip}$ into appropriate $k$
clusters, $C_{1}$, $C_{2}$, …, $C_{k}$, by optimising a given a metric, such
as the within-cluster-sum-of-squares ($W_{k}$).

The method proposed intends to identify the heterogeneity of a performed
cluster analysis by considering the allocation of $x_{ip}$ to a cluster,
$C_{k}$, given the observation membership ($M_{1},M_{2},...,M_{m}$).

The degree of closeness ($\delta_{ik}$) may be considered as the closeness of
point $x_{i}$ being within cluster $C_{k}$ given by,

$\delta_{ik}=\frac{1/{\overline{d^{2}}_{\\!\\!ik}}}{\sum_{k}1/{\overline{d^{2}}_{\\!\\!ik}}}.$
(1)

Where, $\delta_{ik}$ (Equation 1) is calculated as the degree of closeness of
observation $x_{i}$ occurring in cluster $C_{k}$ by considering the mean
distance of $x_{i}$ to all other points allocated to every cluster in $C_{k}$
in Euclidean space, $d_{ij}$. Resulting in,

$\overline{d^{2}}_{\\!\\!ik}=\frac{1}{|n_{k}|}\sum_{j\in k}d^{2}_{ij}.$

Other fuzzy clustering algorithms consider a similar conditional membership,
such as FANNY (Kaufman and Rousseeuw, 2009) and Fuzzy C-Means (Bezdek, 1974;
Dunn, 1973). The degree of closeness of an observation to a cluster is
therefore derived directly from the distance matrix.

$\delta_{ik}=\frac{1/\left({\frac{1}{|n_{k}|}\sum_{j\in
k}d^{2}_{ij}}\right)}{\sum_{k}1/{\left({\frac{1}{|n_{k}|}\sum_{j\in
k}d^{2}_{ij}}\right)}}$ (2)

From these, a consideration of a degree of closeness for each observation to a
cluster (Equation 1) may be given in Equation 2.

Although informative about $x_{i}$ to $C_{k}$, $\delta_{ik}$ provides no
insight into how well the data fits the clusters regarding overall membership
or how many clusters is optimum. The membership ($m$) of each cluster is
considered to be all observations ($x_{i}$) within each cluster $k$. A
summation of all distances within a cluster is considered to form matrix
$\gamma_{mk}=\sum_{i\in m}{\overline{d^{2}}_{\\!\\!ik}}$, over $x_{i}$, for
membership $m$ and cluster $k$. Considerations of the distances of $x_{i}$ are
no longer given, but of the related membership distance over all points to
each cluster. In this instance $x_{i}$ may be considered as a parameter to be
integrated over.

The membership matrix, $1/\gamma_{mk}$, has dimensions of $m$-by-$k$. In the
case of $k=1$, $\gamma_{mk}$ is a $1$-by-$1$ matrix, and the degree of
membership of a cluster occurring given its membership would therefore become
1. With $k=1$ it is assumed there are no other clusters to which membership
may be ascribed ($m=1$). Therefore, this method of optimising may only be
considered at $k>1$. In the case where singleton clusters appear (one
observation to a cluster) $\gamma$ was set to 1 to prevent dividing by zero
and being dominated by singletons.

By considering membership, the degree of closeness regarding the overall
membership of a cluster or that of another cluster may be obtained. The degree
of membership may be considered as

$\delta_{mk}=\frac{1/\gamma_{mk}}{\sum_{k}1/\gamma_{mk}}$ (3)

with an overall measure of validity for each cluster $\delta_{\bullet k}$ is
given below,

$\delta_{\bullet
k}=\frac{\sum_{m}1/\gamma_{mk}}{\sum_{m}\sum_{k}1/\gamma_{mk}}$

and a measure of membership validity,

$\delta_{m\bullet}=\frac{\sum_{k}1/\gamma_{mk}}{\sum_{k}\sum_{m}1/\gamma_{mk}}.$

$\delta_{m\bullet}$ in Equation 3, may be generalised to the distance matrix
as,

$\delta_{mk}=\frac{1/\left(\sum_{i\in k}\left({\frac{1}{|n_{k}|}\sum_{j\in
k}d^{2}_{ij}}\right)\right)}{\sum_{k}1/\left(\sum_{i\in
k}\left({\frac{1}{|n_{k}|}\sum_{j\in k}d^{2}_{ij}}\right)\right)}.$ (4)

##  3 Cluster Optimisation

Considering the $\delta_{mk}$ matrix (Equations 3-4), marginalising over the
clusters’ membership within $\delta_{mk}$ where $k\equiv m$ may be used to
obtain the degree of fit of cluster $C_{k}$ occurring given its own
membership, $M_{m}$, over all clusters.

$\delta_{T}=\sum_{k\equiv m}\delta_{mk}\delta_{m\bullet}$

$\delta_{m\bullet}$ is given as above. Therefore, the degree of fit of the
overall clusters being true given their membership, $\delta_{T}$, may be
derived for every cluster number, $k$, as $k\subset\mathbb{N}$. $\delta_{T}$
may be given as a summary of homogeneity over all clusters, given their
overall membership.

$\delta_{T}$ may be considered to be a method of identifying an adequacy of
the estimated clusters to their respective individual membership. If
$\delta_{T}$ was small, then the estimated clusters may be considered
heterogeneous and therefore inadequate. Similarly, as $\delta_{T}\rightarrow
1$, then the estimated clusters may be considered homogeneous and provide an
adequate fit to their membership.

The degree of fit cluster membership may be represented in terms of ratios, in
order to understand if there has been a large discrepancy in $\delta_{T}$ as
the number of clusters increases. Where, the ratio of clustering at $k$
clusters are,

$\phi_{k}=\frac{(\delta_{T})_{k}}{1-(\delta_{T})_{k}}.$

Optimising the number of clusters by identifying an elbow in the $\phi_{k}$
may be given by maximising the $\phi_{k}$ ratio at lag 1.

$\Phi_{1}=\frac{\phi_{k}}{\phi_{k+1}}$ (5)

Maximising the $\Phi_{1}$ in favour for cluster number $k$ against $k+1$
investigates how likely the clusters fit the data as $k$ increases. It should
be noted if two or more maxima with similar estimates of $\Phi_{1}$ occur at
lag 1, for separate cluster numbers, a direct $\Phi$ estimate and modification
of Equation (5) may be considered.

##  4 Thresholding on Degree of Membership

$\delta_{mk}$ is a $m$-by-$k$ matrix, which forces to use all the distances
within $\gamma_{mk}$. Therefore, the degree of membership of a cluster is
still estimated even if one or many clusters are homogeneous
($\delta_{mk^{\neg}}\rightarrow 0$). These small estimates from homogeneous
clusters may influence the overall estimate of $\delta_{T}$ and $\phi$ ratio.
This may be adjusted for by acknowledging a threshold upon
$\delta_{mk^{\neg}}$ and set the level of $\delta_{mk^{\neg}}$ where it is
believed there to be a more homogeneous cluster. For this study a threshold of
$0.1$ is suggested, where any value of $\delta_{mk^{\neg}}<0.1$ is assumed
$0$.

##  5 Other Cluster Number Optimising Methods

The classical approach to determining the number of clusters is to consider
within-cluster-sum-of-squares ($W_{k}$) as cluster number increases and choose
a number of clusters where there appears to be an elbow.

This method has been discussed and built upon by Tibshirani et al. (2001) to
derive the GAP statistic which considers the 1-standard-error rule by
comparing the log of the pooled within-cluster-sum-of-squares against a
reference distribution derived through multiple sampling. The computational
implementation of the GAP statistic is given in Tibshirani et al. (2001), with
a short description given here. Consider the within-cluster-sum-of-squares,
$W_{k}$, a comparison between the logarithm of $W_{k}$ against the expectation
of $log(\overline{W}_{\\!\\!kb})$. Where, $\overline{W}_{\\!\\!kb}$ is the
average of $b$ Bootstrapped samples.

$GAP_{k}=\frac{1}{b}\sum_{b}log(\overline{W}_{\\!\\!kb})-log(W_{k})$ (6)

The optimum cluster number is then derived from a 1-standard-error rule,
$GAP_{k}\geq GAP_{k+1}-s_{k+1}$. Where, $s_{k+1}$ is the standard error of the
Monte Carlo estimates of $E\\{log(\overline{W}_{kb})\\}$. The method of
determining a reference distribution is given as a uniform distribution or
Principal Component rotating.

Another method of cluster number optimisation includes one proposed by
Caliński and Harabasz (1974). This method compares the between-cluster-sum-of-
squares ($B_{k}$) by $W_{k}$.

$CH_{k}=\frac{B_{k}/(n-k)}{W_{k}/(k-1)}$ (7)

Although the GAP statistic and the Caliński and Harabasz (1974) method provide
good ways to optimise the number of clusters, both only consider the adequacy
of the cluster fit and not the membership of observations to their own
respective clusters or to others.

A method which considers the membership of clusters is the Silhouette
statistic (Rousseeuw, 1987). The Silhouette statistic ($Sl_{k}$) considers the
average distance between points within an observations cluster
($\overline{w}_{k}$) and the average distance between points in the nearest
cluster ($\overline{b}_{k}$).

$Sl_{k}=\frac{\overline{b}_{k}-\overline{w}_{k}}{max\\{\overline{w}_{k},\overline{b}_{k}\\}}$
(8)

A comparison of the different methods given above are compared against the
marginalisation of the degree of membership, and the derivation of $\Phi_{1}$,
at lag 1. This was given to compare which method may give the most appropriate
clustering and where any potential biases occur. The GAP statistic was derived
from the Cluster package (Maechler et al., 2018), and the Silhouette and CH
method from the ClusterCrit package (Desgraupes, 2018).

##  6 Examples

###  6.1 Iris Datset

The method of deriving a degree of membership is applied to the petal width
and length from the Iris dataset given by Anderson (1935). Complete linkage
hierarchical clusters was used to cluster in this example.

The distribution of the petal width and length (Figure 1(a)) show overlap
between clusters 2 ($\square$) and 3 ($\blacksquare$). Both GAP statistics
estimate the optimum cluster number to be 4 (Figure 1(b-c)). The Silhouette
method correctly estimates the number of 3 clusters (Figure 1(d) and the CH
method over estimates the optimum cluster number with 6 (Figure 1(e)). The
maximum $\phi$ ratio at lag 1 occurred at cluster number 3 and had a value of
2.61, with a degree of membership at cluster 3 ($\delta_{T}$) of 0.940 (Figure
1(f)). Therefore, suggesting the clusters are not heterogeneous, regarding the
membership fitted each cluster. Table 1 shows how clusters 2 (centres; 1.87,
5.29) and 3 (1.21, 3.96) are not that well defined compared to cluster 1
(0.25, 1.46). The $\delta_{11}$ was 0.985, which provides a good estimate that
the membership of cluster 1 was represented by its own membership.
$\delta_{22}$ and $\delta_{33}$ were 0.714 and 0.872, respectively, with
$\delta_{23}$ estimated at 0.246 and $\delta_{32}$ estimated at 0.092.
Therefore, the membership of cluster 2 did not fit cluster 2 as well as the
membership of cluster 3 to cluster 3. The lack of representation of the
membership of cluster 2 given its cluster is contributing the most to the
overall $\delta_{T}$ at cluster 3.

###  6.2 Transcriptional Dataset

A further example of applying the degree of membership was to 23 mice and 5011
combined orthologues gene recounts for isolates CB and AS of the malaria
parasite given in Lin et al. (2018). These data include both CB and AS
isolates given to serially blood passaged (SBP) and mosquito-transmitted (MT)
infected mice, with treatments SBP-AS ($n=5$), MT-AS ($n=5$), SBP-CB ($n=5$)
and MT-CB ($n=6$). Complete linkage hierarchical clusters was used to cluster
in this example.

A principal components analysis was conducted on all 5011 genes across 23
mice, the first 3 Principal Components explained 98.89% of the variance.
Principal Components 2 and 3 seem to show the most obvious separation between
treatments of infected mice (Figure 2 (a-b)). Both GAP statistics estimated
the optimum cluster number to be 1, suggesting no clusters over 3 Principal
Components (Figure 2 (c-d)). The maximum $\phi$ ratio at lag 1 occurred at 4
clusters, with a value of 2.511 (Figure 2 (e)). The estimated degree of
membership at 4 clusters was 0.647 (Figure 2 (f)), which suggests these
clusters identified within the transcriptome may not well defined over 3
Principal Components and that there is noise between clusters within the
dataset. The treatment group of infected mice SBP-AS was the only group which
came out as their own separate cluster (Cluster 1 $\bullet$; Figure 3). Where,
$\delta_{11}$ had the highest estimated value (0.695) compared to all other
clusters regarding their membership ($\delta_{22}=0.660$, $\delta_{33}=0.596$
and $\delta_{44}=0.662$) (Table 2 (a)). Clusters 3 and 4 were the least
defined and had the most overlap, with $\delta_{34}$ estimated as 0.254. With
thresholding upon $\delta_{mk}$ of 0.01, the optimum number of clusters is
still 4, where the estimated $\phi$ ratio was 3.256 and $\delta_{T}$ increases
to 0.745, which suggests some clusters are well defined. However, clusters 3
and 4 are still not well separated, with $\delta_{34}$ now estimated as 0.299
(Table (2)b).

##  7 Simulation Study

Below is a simulation study, similar to Tibshirani et al. (2001) when
evaluating the GAP statistic. This was applied to the following methods
described in Section 5 and the $\phi$ ratio presented in Section 3.

  1. (1)

Gap statistic using a uniform distribution (Tibshirani et al., 2001) (Equation
6),

  2. (2)

Gap statistic using PCA (Tibshirani et al., 2001),

  3. (3)

Caliński and Harabasz (1974) (Equation 7),

  4. (4)

Silhouette (Rousseeuw, 1987) (Equation 8),

  5. (5)

$\phi$ ratio at lag 1 (Equation 5),

  6. (6)

$\phi$ ratio at lag 1 with threshold factor upon $\delta_{mk}$ at 0.1.

Four simulation scenarios are given and have been sampled 100 times. The
results from these simulations are presented in Table 3. Complete linkage
hierarchical clustering was used.

  1. (a)

A uniform distribution with 100 observations over six dimensions with a
minimum of 0 and a maximum of 1.

  2. (b)

Four clusters four dimensions with centres ($-4,-6$), ($-8,1$), ($4,-5$) and
($-6,9$). The first simulation has a sample size of 25, 25, 10 and 10 for each
cluster, the second simulation has a sample of 5 for all clusters. This was to
test the adequacy of the proposed method of degree of membership with a low
number of samples.

  3. (c)

Four clusters in four dimensions with centres ($11,-8,0,-3$), ($-8,4,4,2$),
($9,$ $-2,$ $-2,$ $3$), ($3$, $7,$ $-4,$ $0$) with 30, 30, 20 and 15 samples.
To evaluate the bias of higher dimensional datasets.

  4. (d)

Nested Clusters with two nesting at each cluster, 6 in total, with centres
($12,-15$), ($15,-18$), ($-16,-15$), ($-16,-18$), ($17,14$) and ($14,11$) with
$25,$ $25,$ $15,$ $15,$ $10$ and $10$ samples. First sampled with a standard
deviation of 0.5, then 1, to illustrate the effect when nested clusters (SD =
0.5) become a single cluster (SD = 1). At SD = 0.5, both 3 and 6 clusters are
thought to be equally likely.

Uniform distribution: Both GAP methods adequately found there were no
clusters. This is due to the Silhouette, CH and $\phi$ ratio unable to be
calculated for 1 cluster. Both the Silhouette and the CH both had a tendency
to over cluster, whereas the $\phi$ ratio optimised around a total cluster
number of 2.

Four clusters two dimensions: All methods found 4 to be the optimum number of
clusters. The $\phi$ ratio and both GAP methods were the worst in estimating
the correct number of clusters compared to all other methods. The Silhouette,
CH $\phi$ ratio with threshold provided the best methods in this simulation.
With fewer observations, the Silhouette and $\phi$ ratio with threshold were
the better methods at estimating the correct number of clusters. The CH method
became worse at estimating the correct number of clusters with fewer
observations. In both cases, estimation on the optimum number of clusters
through the $\phi$ ratio was improved by including a threshold at 0.1.

Four clusters, four dimensions: With higher dimensions, the $\phi$ ratio was
similar at optimising the number of clusters at 2 dimensions. The GAP, unif,
PCA and CH performed better in higher dimensions. The Silhouette method
provided worst estimates in more dimensions. The $\phi$ ratio with threshold
of 0.1 still provided an optimisation of 4 clusters.

Nested clusters: The $\phi$ ratio with threshold was the least biased when
optimising the number of clusters, with a 49-50 split for 3 and 6 clusters,
respectively at SD = 0.5. The Silhouette method had a bias for fewer clusters.
With a more variable dataset (SD = 1), and 6 clusters becomes 3, the GAP
statistics and CH methods over-performed compared to the Silhouette and $\phi$
ratio. By applying a threshold of 0.1 to the method of $\phi$ ratio for nested
clusters, the method accurately identified both 3 clusters and not their
nests. With increased variability (SD = 1), and therefore the disappearance of
nested clusters, the $\phi$ ratio with and without a threshold and the
Silhouette method outperformed the GAP statistics and the CH method.

There was no overall consensus on which method out-performed the others across
all simulations (Table 3), as each method is suited to different types of
data, where the GAP statistics are the only methods which may consider no
clusters. What is consistent is that the $\phi$ ratio estimate was more
consistent across simulations. With the nested simulation suiting the $\phi$
ratio best compared to other methods. Considerations of the definition of a
defined cluster and threshold should also be given for different cluster
scenarios.

##  8 Discussion

The use of the $\phi$ ratio at lag 1 (Equation 5) and other established
clustering statistics (Equations 6-8) provided an adequate estimation of
optimisation of the number of clusters for all scenarios where $k>1$. Although
some have potential biases and some work better with different types of
datasets. The method of $\phi$ ratio was shown to be the least bias when
identifying nested clusters and with more variable datasets. Generally the
$\phi$ ratio was more consistent across given scenarios. The degree of
membership provided an quantitative estimate of the homogeneity, not only for
the whole cluster and membership fit, but also for each cluster. This
statistic makes it possible to identify which clusters within an analysis may
not be representative of their membership, or how well their membership is
represented by other clusters. The Silhouette method considers a similar
method but only of the neighbouring cluster. It may be the case that from a
cluster analysis, more so with nested clusters, that an optimum number of
clusters may be found, however, the membership of closely related clusters may
blur the accuracy of identified clusters. It is proposed the heterogeneity may
be evaluated for each cluster and each membership through $\delta_{mk}$.

The method of optimising cluster number through the $\phi$ ratio was improved
by introducing a threshold upon $\delta_{mk^{\neg}}$ for defined clusters.
Although cluster analysis may not be an application-independent mathematical
problem (Von Luxburg et al., 2012), where prior experimental knowledge may
explain clusters, each cluster and its membership within one analysis is not
application-independent, such that some clusters may be more biologically
relevant and obvious in their membership, while others may be more varied and
so a consideration of a threshold upon $\delta_{mk}$ should be considered.
However, considerations on the threshold should be given as this may vary
depending on the dimensionality and if any nested clusters are present within
the dataset. Overall the use of the same threshold for all cluster analyses
was shown to be inadequate and varied between simulations. Therefore, similar
to the application of cluster analysis (Von Luxburg et al., 2012), the use of
a threshold upon $\delta_{mk}$ is not application-independent and is dependent
on the relevance of the applied problem and leads to more of a problem-centric
perspective.

The application of classifying a dataset into optimal meaningful clusters is a
difficult process dependent on the sampled data, method of clustering and
method of optimising used. This paper presents that considerations of the
heterogeneity of each cluster through cluster membership should be
acknowledged. The method gives consideration of how membership of a cluster
should also be considered along with the cluster fit itself.

##  9 Acknowledgements

A special thanks to JJ Valletta, K Hassall, V Koutra, H Metcalf and G Kelly
for reading an early draft of this manuscript and providing valuable feedback.

The Francis Crick Institute receives its core funding from Wellcome
(FC001101), the UK Medical Research Council (FC001101) and Cancer Research UK
(FC001101). This study was funded by the UK Medical Research Council
(MR/M003906/1). Jean Langhorne is a Wellcome Trust Senior investigator
(102907).

## References

  * (1)
  * Anderson (1935) Anderson, E. (1935), ‘The irises of the gaspe peninsula’, Bulletin of the American Iris Society 59, 2–5.
  * Bezdek (1974) Bezdek, J. C. (1974), ‘Cluster validity with fuzzy sets’, Journal of Cybernetics 3, 58–73.
  * Caliński and Harabasz (1974) Caliński, T. and Harabasz, J. (1974), ‘A dendrite method for cluster analysis’, Communications in Statistics-theory and Methods 3(1), 1–27.
  * Desgraupes (2018) Desgraupes, B. (2018), clusterCrit: Clustering Indices.  R package version 1.2.8.
https://CRAN.R-project.org/package=clusterCrit

  * Dunn (1973) Dunn, J. C. (1973), ‘A fuzzy relative of the isodata process and its use in detecting compact well-separated clusters’, Journal of Cybernetics 3, 32–57.
  * Kaufman and Rousseeuw (2009) Kaufman, L. and Rousseeuw, P. J. (2009), Finding groups in data: an introduction to cluster analysis, Vol. 344, John Wiley & Sons.
  * Lin et al. (2018) Lin, J.-W., Reid, A. J., Cunningham, D., Böhme, U., Tumwine, I., Keller-Mclaughlin, S., Sanders, M., Berriman, M. and Langhorne, J. (2018), ‘Genomic and transcriptomic comparisons of closely related malaria parasites differing in virulence and sequestration pattern’, Wellcome Open Research 3.
  * Maechler et al. (2018) Maechler, M., Rousseeuw, P. J., Struyf, A., Hubert, M. and Hornik, K. (2018), cluster: Cluster Analysis Basics and Extensions.  R package version 2.0.7-1.
  * R Core Team (2018) R Core Team (2018), R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria.
https://www.R-project.org/

  * Rousseeuw (1987) Rousseeuw, P. J. (1987), ‘Silhouettes: a graphical aid to the interpretation and validation of cluster analysis’, Journal of computational and applied mathematics 20, 53–65.
  * Tibshirani et al. (2001) Tibshirani, R., Walther, G. and Hastie, T. (2001), ‘Estimating the number of clusters in a data set via the gap statistic’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63(2), 411–423.
  * Von Luxburg et al. (2012) Von Luxburg, U., Williamson, R. C. and Guyon, I. (2012), Clustering: Science or art?, in ‘Proceedings of ICML Workshop on Unsupervised and Transfer Learning’, pp. 65–79.